Download data

d <- read.table("http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts.txt",
                skip = 7, 
                header = TRUE, 
                nrows = 143, 
                colClasses = "character")

Clean data

d <- d[-grep("Year", d$Year),]
rownames(d) <- d[,1]
d <- d[,-1]
d <- d[,-13:-19]
d[d == "****"] <- NA

# convert to numeric and celsius
for (i in 1:12) {
    d[,i] <- as.numeric(d[,i]) / 100
}

Download seasonal adjustments

s <- read.table("http://data.giss.nasa.gov/gistemp/faq/merra2_seas_anom.txt",
                skip = 3, 
                header = TRUE, 
                colClasses = "character")
sa <- as.numeric(s$seas_anom)

# create colours from blue through to red
colours <- colorRampPalette(c("grey", "blue", "red"))(nrow(d))

Create initial plot

mmin <- -3
mmax <- 3
par(mar = c(2, 2, 2, 0))
plot(1:12, d[1,]+sa, type="l", ylim=c(mmin,mmax), yaxt="n", xaxt="n", bty="n")

# add axes
axis(side = 1, 
     at = 1:12, 
     lwd = 0, 
     lwd.ticks = 1, 
     cex.axis = 0.8,
     labels = c("Jan","Feb","Mar",
                "Apr","May","Jun",
                "Jul","Aug","Sep",
                "Oct","Nov","Dec"))

axis(side = 2, at = -3:3, labels = -3:3, 
     lwd = 0, lwd.ticks = 1, las = 2, 
     cex.axis = 0.8)

# add title
title(main = expression(paste("Temperature ","Anomaly ","(",degree,"C)")),
      adj = 0.5, 
      cex.main = 0.8)
title(main = "(Difference from 1980-2015 annual mean)", 
      adj = 0.5, 
      cex.main = 0.6, 
      line = 0, 
      font.main = 1)

# add horizontal lines
for (j in -3:3) {
    lines(1:12, rep(j, 12), lty = 3)
}

# add yearly temperature lines
for (j in 1:nrow(d)) {
    lines(1:12, as.numeric(d[j,]) + sa, col = colours[j])
}