# joy style plot for gistemp # output into pdf or png by uncommenting the relevant lines pdf()/dev.off. Last frame needs to be done by hand # animation uses "convert -layers OptimizePlus -strip -density 150 -delay 30 joyseas5_1.pdf joyseas5_2.pdf joyseas5_3.pdf joyseas5_4.pdf joyseas5_5.pdf joyseas5_6.pdf joyseas5_7.pdf joyseas5_8.pdf joyseas5_9.pdf joyseas5_10.pdf joyseas5_11.pdf joyseas5_12.pdf joyseas5_13.pdf joyseas5_14.pdf joyseas5_15.pdf joyseas5_16.pdf joyseas5_17.pdf joyseas5_18.pdf joyseas5_19.pdf joyseas5_20.pdf joyseas5_21.pdf joyseas5_22.pdf -delay 200 joyseas5_N.pdf joyseas5.gif" or similar # inputs (should stay up-to-date for the rest of 2017) giss=read.csv("https://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.csv",skip=1,colClasses = "character",header=TRUE) # make a month x year array of anomalies bymon=array(as.numeric(as.character(t(giss[1:137,2:13]))),c(12,137)) # seasonal cycle anomalies merra=read.table("https://data.giss.nasa.gov/gistemp/faq/merra2_seas_anom.txt",skip=3,header=TRUE) merra_seas=merra[,3] baseseas=rowMeans(bymon[,102:131]) # baseline shift for MERRA climatology bymonseas=bymon+merra_seas-baseseas # mon=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec") ########## # joy/slug frames for basic anomalies # set number of years in each step dyr=5. # 2 or 5 or 10 ns=(136 - 30) %/% dyr +1. # number of steps # define a color bar that is smooth in rgb ncol=ns+1 colfr=NA*(0:ncol) for (i in 0:ncol) {ic=i/ncol ; colfr[i+1]=rgb(ic,3*(1-ic)*ic,1-ic)} for (j in 1:ns) { # pdf(paste("joy",dyr,"_",j,".pdf",sep="")) par(xpd=FALSE) plot(c(-1.5,1.5),c(1,13),type="n",xlab="Temperature (w.r.t. 1951-1980)",ylab="",yaxt="n",bty="n",family="Helvetica",main="How have temperature distributions changed since the 19th Century?") for (i in 1:12) { polygon(c(density(bymon[i,2:31])$x,1.5,-1.5),c(0.5*density(bymon[i,2:31])$y+(13-i),13-i,13-i),col=colfr[1],lwd=0.5) if (j > 1) { polygon(c(density(bymon[i,(j-1)*dyr+2:31])$x,1.5,-1.5),c(0.5*density(bymon[i,(j-1)*dyr+2:31])$y+(13-i),13-i,13-i),col=colfr[j],lwd=0.5) legend(0,13.8,box.lwd=NA,legend=paste(1881+(j-1)*dyr,"-",1910+(j-1)*dyr)) } legend(-1.67,13.,legend=mon,box.lwd=NA,y.intersp=2,text.col=1) legend(-0.8,13.8,box.lwd=NA,legend="1881-1910") axis(1,col.ticks=2,col=1,col.axis=2,padj=1.6,labels=c("-2.0","-1.0","0.0","1.0","2.0"),at=c(-1.1111,-0.5555,0,0.5555,1.1111),xlab="(ºF)") par(xpd=TRUE) text(1.64,-0.25,"ºC") text(1.64,-0.8,"ºF",col=2) } #dev.off() readline(prompt="Press [enter] to continue") } # plot 2016 and 2017 separately # for output redo last climatology than this points(as.numeric(giss[137,2:13]),13.12-(1:12),pch=0,col=1,cex=1.4) points(as.numeric(giss[137,2:13]),13.12-(1:12),pch=15,col=2,cex=1.4) points(as.numeric(giss[138,2:13]),13.12-(1:12),pch=15,col=1,cex=1.4) points(as.numeric(giss[138,2:13]),13.12-(1:12),pch=0,col=1,cex=1.4) legend(0.9,13.25,cex=1.,legend=c("2016"),text.col=2,box.lwd=NA) legend(0.65,13.25,cex=1.,legend=c("2017"),,text.col=1,box.lwd=NA) #dev.copy(pdf,paste("joy",dyr,"_N.pdf",sep="")) #dev.off() readline(prompt="Press [enter] to continue") # now with seasonal climatology as well. for (j in 1:ns) { # pdf(paste("joyseas",dyr,"_",j,".pdf",sep="")) par(xpd=FALSE) plot(c(-3.7,3.7),c(1,13),type="n",xlab="Temperature (w.r.t. 1981-2010)",ylab="",yaxt="n",bty="n",family="Helvetica",main="How have seasonal temperatures changed since the 19th Century?",xlim=c(-3.7,3.7),xaxs="i") for (i in 1:12) { polygon(c(density(bymonseas[i,2:31])$x,3.7,-3.7),c(0.5*density(bymonseas[i,2:31])$y+(13-i),13-i,13-i),col=colfr[1],lwd=0.5) if (j > 1) { polygon(c(density(bymonseas[i,(j-1)*dyr+2:31])$x,3.7,-3.7),c(0.5*density(bymonseas[i,(j-1)*dyr+2:31])$y+(13-i),13-i,13-i),col=colfr[j],lwd=0.5) legend(0,13.8,box.lwd=NA,legend=paste(1881+(j-1)*dyr,"-",1910+(j-1)*dyr)) } legend(-4.1,13.,legend=mon,box.lwd=NA,y.intersp=2,text.col=1) legend(-2.5,13.8,box.lwd=NA,legend="1881-1910") axis(1,col.ticks=2,col=1,col.axis=2,padj=1.6,labels=c("-4","-2","0","2","4"),at=c(-2.2222,-1.1111,0,1.1111,2.2222),xlab="(ºF)") par(xpd=TRUE) text(3.2,-0.25,"ºC") text(3.2,-0.8,"ºF",col=2) } #dev.off() readline(prompt="Press [enter] to continue") } # 2016/2017 separately points(as.numeric(giss[137,2:13])+merra_seas-baseseas,13.12-(1:12),pch=0,col=1,cex=1.4) points(as.numeric(giss[137,2:13])+merra_seas-baseseas,13.12-(1:12),pch=15,col=2,cex=1.4) points(as.numeric(giss[138,2:13])+merra_seas-baseseas,13.12-(1:12),pch=15,col=1,cex=1.4) points(as.numeric(giss[138,2:13])+merra_seas-baseseas,13.12-(1:12),pch=0,col=1,cex=1.4) legend(-1.6,13.25,cex=1.,legend=c("2016"),text.col=2,box.lwd=NA) legend(-2.1,13.25,cex=1.,legend=c("2017"),,text.col=1,box.lwd=NA) #dev.copy(pdf,paste("joyseas",dyr,"_N.pdf",sep="")) #dev.off()