#Why R for NOAA/NCEI: Demonstrated by Examples? #By Sam Shen, San Diego State Univeristy #May 29, 2017 #Display the trends of 12 months using 12 panels #Read directly from the NCEI URL for the NOAA GlobalTemp monthly global aevrage data: #https://www1.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/timeseries/aravg.mon.land_ocean.90S.90N.v4.0.1.201704.asc setwd("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch12-13-Rgraphics") aveNCEI <- read.table("https://www1.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/timeseries/aravg.mon.land_ocean.90S.90N.v4.0.1.201704.asc", header=FALSE) dim(aveNCEI) #Jan 1880-Feb 2017 #an extra month to be deleted #[1] 1648 10 par(mar=c(4,5,2,1)) timemo=seq(aveNCEI[1,1],aveNCEI[length(aveNCEI[,1]),1],len=length(aveNCEI[,1]) ) #create matrix of 136 yrs of data matrix # row=year from 1880 to 2016, col=mon 1 to 12 ave=aveNCEI[,3] myear=length(ave)/12 nyear=floor(myear) nmon=nyear*12 avem = matrix(ave[1:nmon], ncol=12, byrow=TRUE) #compute annual average annv=seq(0,length=nyear) for(y in 1:nyear){annv[y]=mean(avem[y,])} #Put the monthly averages and annual ave in a matrix avemy=cbind(avem,annv) #Plot 12 panels on the same figure: Trend for each month dev.off() quartz(width=10,height=16,pointsize = 16) #quartz(display = "name", width = 5, height = 5, pointsize = 12, # family = "Helvetica", antialias = TRUE, autorefresh = TRUE) plot.new() #png(file = 'monthtrend.png') #Automatical saving of a figure timeyr=seq(aveNCEI[1,1], aveNCEI[1,1]+nyear-1) par(mfrow = c(4, 3)) # 4 rows and 3 columns par(mgp=c(2,1,0)) for (i in 1:12) { plot(timeyr, avemy[,i],type="l", ylim=c(-1.0,1.0), xlab="Year",ylab="Temp [oC]", main = paste("Month =", i, split = "")) abline(lm(avemy[,i]~timeyr),col="red") text(1945,0.7, paste("Trend oC/century=", round(digits=3, (100*coefficients(lm(avemy[,i]~timeyr))[2]))), col="red") } #Reproduce the IPCC plot of the NOAAGlobalTemp time series NOAATemp <- read.table("https://www1.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/timeseries/aravg.ann.land_ocean.90S.90N.v4.0.1.201704.asc", header=FALSE) plot.new() dev.off() quartz(width=10,height=7,pointsize = 12) par(mar=c(4,4,3,1)) x<-NOAATemp[,1] y<-NOAATemp[,2] z<-rep(-99,length(x)) for (i in 3:length(x)-2) z[i]=mean(c(y[i-2],y[i-1],y[i],y[i+1],y[i+2])) n1<-which(y>=0) x1<-x[n1] y1<-y[n1] n2<-which(y<0) x2<-x[n2] y2<-y[n2] x3<-x[2:length(x)-2] y3<-z[2:length(x)-2] plot(x1,y1,type="h",xlim=c(1880,2016),lwd=3, tck=0.02, ylim=c(-0.7,0.7), #tck>0 makes ticks inside the plot ylab="Temperature [deg C]", xlab="Time",col="red", main="NOAA Global Average Annual Mean Temperature Anomalies") lines(x2,y2,type="h", lwd=3, tck=-0.02, col="blue") lines(x3,y3,lwd=2) #Plot a temperature map from NOAAGlobalTemp gridded data #Read the data online from an URL #Plot the climate data map for a given time #Install maps package if not done before #install.packages("maps") gpcpst=read.csv("http://shen.sdsu.edu/data/NOAAGlobalTJ1880_J2017.csv") library(maps) Lat= seq(-87.5, 87.5, length=36) Lon=seq(2.5, 357.5, length=72) mapmat=matrix(gpcpst[,1634],nrow=72) #column 1634 corresponding to Dec 2015 #Covert the vector into a lon-lat matrix for R map plotting mapmat=pmax(pmin(mapmat,6),-6) #This command compresses numbers larger than 6 to 6 dev.off() quartz(width=14,height=7,pointsize = 12) plot.new() par(mar=c(4,5,3,0)) int=seq(-6,6,length.out=81) rgb.palette=colorRampPalette(c('black','blue', 'darkgreen','green', 'yellow','pink','red','maroon'),interpolate='spline') mapmat= mapmat[,seq(length(mapmat[1,]),1)] filled.contour(Lon, Lat, mapmat, color.palette=rgb.palette, levels=int, plot.title=title(main="NOAAGlobalTemp Anomalies Dec 2015 [deg C]", xlab="Latitude",ylab="Longitude", cex.lab=1.5), plot.axes={axis(1, cex.axis=1.5); axis(2, cex.axis=1.5);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]"), key.axes={axis(4, cex.axis=1.5)}) #Animmation of the NOAAGlobalTemp gpcpst=read.csv("http://shen.sdsu.edu/data/NOAAGlobalTJ1880_J2017.csv") library(maps) install.packages("animation") library(animation) int=seq(-5,5,length.out=81) rgb.palette=colorRampPalette(c('black','blue','darkgreen','green', 'white','yellow','pink','red','maroon'), interpolate='spline') n1 = length(gpcpst[1,])-15 n2 = length(gpcpst[1,])-3 Lat= seq(-87.5, 87.5, length=36) Lon=seq(2.5, 357.5, length=72) ## set up an empty frame, then add points one by one par(bg = "white") # ensure the background color is white ani.record(reset = TRUE) # clear history before recording dev.off() quartz(width=12,height=6,pointsize = 12) for (i in n1:n2) { precnc=matrix(gpcpst[,i],nrow=72) mapmat=pmax(pmin(precnc,5),-5) filled.contour(Lon, Lat, mapmat, zlim=c(-5,5), color.palette=rgb.palette, levels=int, plot.title=title( main=paste(c("Monthly surface air temperature anomalies [deg C]:"), (colnames(gpcpst))[i]), xlab="Longitude",ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) ani.record() # record the current frame } ## now we can replay it, with an appropriate pause between frames oopts = ani.options(interval = 0.6) ani.replay() ## or export the animation to an HTML page saveHTML(ani.replay(), img.name = "NOAAGlobalT", width=12, height=6) #Animation for El Nino regions #Animmation of the NOAAGlobalTemp gpcpst=read.csv("http://shen.sdsu.edu/data/NOAAGlobalTJ1880_J2017.csv") library(maps) install.packages("animation") library(animation) int=seq(-5,5,length.out=81) rgb.palette=colorRampPalette(c('black','blue','darkgreen','green', 'white','yellow','pink','red','maroon'), interpolate='spline') n1 = length(gpcpst[1,])-15 n2 = length(gpcpst[1,])-3 Lat= seq(-87.5, 87.5, length=36) Lon=seq(2.5, 357.5, length=72) ## set up an empty frame, then add points one by one par(bg = "white") # ensure the background color is white ani.record(reset = TRUE) # clear history before recording dev.off() quartz(width=12,height=6,pointsize = 12) for (i in n1:n2) { precnc=matrix(gpcpst[,i],nrow=72) mapmat=pmax(pmin(precnc,5),-5) filled.contour(Lon, Lat, mapmat, zlim=c(-5,5), xlim=c(100,300),ylim=c(-60,60), color.palette=rgb.palette, levels=int, plot.title=title( main=paste(c("Monthly surface air temperature anomalies [deg C]:"), (colnames(gpcpst))[i]), xlab="Longitude",ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) ani.record() # record the current frame } ## now we can replay it, with an appropriate pause between frames oopts = ani.options(interval = 0.6) ani.replay() ## or export the animation to an HTML page saveHTML(ani.replay(), img.name = "NOAAGlobalT", width=12, height=6)