#R Programming for Climate Data Analysis and Visualization #Codes #A Short Course at NOAA-NCEI, May 30-June 2, 2017 #By Samuel Shen #San Diego State University #Email: sam.shen@sdsu.edu #The same sequence can be generated by different commands, such as 1:8 seq(1,8) seq(8) seq(1,8, by=1) seq(1,8, length=8) seq(1,8, length.out =8) #Define a function samfctn <- function(x) x*x samfctn(4) #[1] 16 fctn2 <- function(x,y,z) x+y-z/2 fctn2(1,2,3) #[1] 1.5 #Simple plots plot(sin, -pi, 2*pi) #plot the curve of y=sin(x) from -pi to 2 pi square <- function(x) x*x #Define a function plot(square, -3,2) # Plot the defined function # Plot a 3D surface x <- seq(-1, 1, length=100) y <- seq(-1, 1, length=100) z <- outer(x, y, function(x, y)(1-x^2-y^2)) #outer (x,y, function) renders z function on the x, y grid persp(x,y,z, theta=330) # yields a 3D surface with perspective angle 330 deg #Contour plot contour(x,y,z) #lined contours filled.contour(x,y,z) #color map of contours #Symbolic calculation D(expression(x^2,'x'), 'x') # Take derivative of x^2 w.r.t. x 2 * x #The answer is 2x fx= expression(x^2,'x') #assign a function D(fx,'x') #differentiate the function w.r.t. x 2 * x #The answer is 2x fx= expression(x^2*sin(x),'x') #Change the expression and use the same derivative command D(fx,'x') 2 * x * sin(x) + x^2 * cos(x) fxy = expression(x^2+y^2, 'x','y') #One can define a function of 2 or more variables fxy #renders an expression of the function in terms of x and y #expression(x^2 + y^2, "x", "y") D(fxy,'x') #yields the partial derivative with respect to x: 2 * x D(fxy,'y') #yields the partial derivative with respect to y: 2 * y square = function(x) x^2 integrate (square, 0,1) #Integrate x^2 from 0 to 1 equals to 1/3 with details below #0.3333333 with absolute error < 3.7e-15 integrate(cos,0,pi/2) #Integrate cos(x) from 0 to pi/2 equals to 1 with details below #1 with absolute error < 1.1e-14 #Vectors, matrices and solve linear equations c(1,6,3,pi,-3) #c() gives a vector and is considered a 4X1 column vector #[1] 1.000000 6.000000 3.000000 3.141593 -3.000000 seq(2,6) #Generate a sequence from 2 to 6 #[1] 2 3 4 5 6 seq(1,10,2) # Generate a sequence from 1 to 10 with 2 increment #[1] 1 3 5 7 9 x=c(1,-1,1,-1) x+1 #1 is added to each element of x #[1] 2 0 2 0 2*x #2 multiplies each element of x #[1] 2 -2 2 -2 x/2 # Each element of x is divided by 2 #[1] 0.5 -0.5 0.5 -0.5 y=seq(1,4) x*y # This multiplication * multiples each pair of elements #[1] 1 -2 3 -4 x%*%y #This is the dot product of two vectors and yields # [,1] #[1,] -2 t(x) # Transforms x into a row 1X4 vector # [,1] [,2] [,3] [,4] #[1,] 1 -1 1 -1 t(x)%*%y #This is equivalent to dot product and forms 1X1 matrix # [,1] #[1,] -2 x%*%t(y) #This column times row yields a 4X4 matrix # [,1] [,2] [,3] [,4] #[1,] 1 2 3 4 #[2,] -1 -2 -3 -4 #[3,] 1 2 3 4 #[4,] -1 -2 -3 -4 my=matrix(y,ncol=2) #Convert a vector into a matrix of the same number of elements #The matrix elements go by column, first column, second, etc #Commands matrix(y,ncol=2, nrow=2) or matrix(y,2) #or matrix(y,2,2) does the same job my # [,1] [,2] #[1,] 1 3 #[2,] 2 4 dim(my) #find dimensions of a matrix #[1] 2 2 as.vector(my) #Convert a matrix to a vector, again via columns #[1] 1 2 3 4 mx <- matrix(c(1,1,-1,-1), byrow=TRUE,nrow=2) mx*my #multiplication between each pair of elements # [,1] [,2] #[1,] 1 3 #[2,] -2 -4 mx/my #division between each pair of elements # [,1] [,2] #[1,] 1.0 0.3333333 #[2,] -0.5 -0.2500000 mx-2*my # [,1] [,2] #[1,] -1 -5 #[2,] -5 -9 mx%*%my #This is the real matrix multiplication in matrix theory # [,1] [,2] #[1,] 3 7 #[2,] -3 -7 det(my) #determinant #[1] -2 myinv = solve(my) #yields the inverse of a matrix myinv # [,1] [,2] #[1,] -2 1.5 #[2,] 1 -0.5 myinv%*%my #verifies the inverse of a matrix # [,1] [,2] #[1,] 1 0 #[2,] 0 1 diag(my) #yields the diagonal vector of a matrix #[1] 1 4 myeig=eigen(my) #yields eigenvalues and unit eigenvectors myeig myeig$values #[1] 5.3722813 -0.3722813 myeig$vectors # [,1] [,2] #[1,] -0.5657675 -0.9093767 #[2,] -0.8245648 0.4159736 mysvd = svd(my) #SVD decomposition of a matrix M=UDV' #SVD can be done for a rectangular matrix of mXn mysvd$d #[1] 5.4649857 0.3659662 mysvd$u # [,1] [,2] #[1,] -0.5760484 -0.8174156 #[2,] -0.8174156 0.5760484 mysvd$v # [,1] [,2] #[1,] -0.4045536 0.9145143 #[2,] -0.9145143 -0.4045536 ysol=solve(my,c(1,3)) #solve linear equations matrix %*% x = b ysol #solve(matrix, b) #[1] 2.5 -0.5 my%*%ysol #verifies the solution # [,1] #[1,] 1 #[2,] 3 #Basic statistics of a sequence x=rnorm(10) #generate 10 normally distributed numbers x #[1] 2.8322260 -1.2187118 0.4690320 -0.2112469 0.1870511 #[6] 0.2275427 -1.2619005 0.2855896 1.7492474 -0.1640900 mean(x) #[1] 0.289474 var(x) #[1] 1.531215 sd(x) #[1] 1.237423 median(x) #[1] 0.2072969 quantile(x) # 0% 25% 50% 75% 100% #-1.2619005 -0.1994577 0.2072969 0.4231714 2.8322260 range(x) #yields the min and max of x #[1] -1.261900 2.832226 max(x) #[1] 2.832226 boxplot(x) #yields the box plot of x w=rnorm(1000) summary(rnorm(12)) #statistical summary of the data sequence # Min. 1st Qu. Median Mean 3rd Qu. Max. #-1.9250 -0.6068 0.3366 0.2309 1.1840 2.5750 hist(w) #yields the histogram of 1000 random numbers with a normal distribution #Read the 5 deg gridded NOAAGlobalTemp .asc data rm(list=ls(all=TRUE)) # Download .asc file and read it into R #setwd("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/NOAAGlobalTemp") #da1=scan("NOAAGlobalTemp.gridded.v4.0.1.201701.asc") #Or read it directly from the data website URL da1=scan("http://shen.sdsu.edu/data/NOAAGlobalTemp.gridded.v4.0.1.201701.asc") length(da1) #[1] 4267130 da1[1:3] #[1] 1.0 1880.0 -999.9 #means mon, year, temp #data in 72 rows (2.5, ..., 357.5) and #data in 36 columns (-87.5, ..., 87.5) tm1=seq(1,4267129, by=2594) tm2=seq(2,4267130, by=2594) length(tm1) length(tm2) mm1=da1[tm1] #Extract months yy1=da1[tm2] #Extract years head(mm1) head(yy1) length(mm1) length(yy1) rw1<-paste(yy1, sep="-", mm1) #Combine YYYY with MM head(tm1) head(tm2) tm3=cbind(tm1,tm2) tm4=as.vector(t(tm3)) head(tm4) #[1] 1 2 2595 2596 5189 5190 da2<-da1[-tm4] #Remote the months and years data from the scanned data length(da2)/(36*72) #[1] 1645 #months, 137 yrs 1 mon: Jan 1880-Jan 2017 da3<-matrix(da2,ncol=1645) #Generate the space-time data #2592 (=36*72) rows and 1645 months (=137 yrs 1 mon) #Convert the data into standard space-time data, #Add the latitude and longitude coordinates for each #grid box as the first two columns, and the time mark #for each month as the first row colnames(da3)<-rw1 lat1=seq(-87.5, 87.5, length=36) lon1=seq(2.5, 357.5, length=72) LAT=rep(lat1, each=72) LON=rep(lon1,36) gpcpst=cbind(LAT, LON, da3) head(gpcpst) dim(gpcpst) #[1] 2592 1647 #The first two columns are Lat and Lon #-87.5 to 87.5 and then 2.5 to 375.5 #The first row for time is header, not counted as data. write.csv(gpcpst,file="NOAAGlobalT.csv") #Output the data as a csv file #Plot the climate data map for a given time #Install maps package if not done before #install.packages("maps") 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 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)}) #Extract the data for a specific region #Keep only the data for the Pacific region n2<-which(gpcpst[,1]>-20&gpcpst[,1]<20&gpcpst[,2]>160&gpcpst[,2]<260) dim(gpcpst) length(n2) #[1] 160 #4 latitude bends and 20 longitude bends pacificdat=gpcpst[n2,855:1454] #from 1951-2000 #(1951-1880)*12 + lat col + lon col =854 #Thus, Jan 1951 data from column 855 #Plot a climate map for a specific region Lat=seq(-17.5,17.5, by=5) Lon=seq(162.5, 257.5, by=5) plot.new() par(mar=c(4,5,3,0)) mapmat=matrix(pacificdat[,564], nrow=20) int=seq(-5,5,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, xlim=c(120,300),ylim=c(-40,40), plot.title=title(main="Tropic Pacific SAT Anomalies [deg C]: Dec 1997", 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)}) #Extract data from only one grid box #Extract data for a specified box with given lat and lon n2 <- which(gpcpst[,1]==32.5&gpcpst[,2]==242.5) SanDiegoData <- gpcpst[n2,855:1454] plot(seq(1880,2017, len=length(SanDiegoData)), SanDiegoData, type="l", xlab="Year", ylab="Temp [oC]", main="San Diego temperature anomalies history") #Fill the few missing data temp=gpcpst areaw=matrix(0,nrow=2592,ncol = 1647) dim(areaw) #[1] 2592 1647 areaw[,1]=temp[,1] areaw[,2]=temp[,2] veca=cos(temp[,1]*pi/180) #Here, don't forget to convert degrees to radians #Area-weight matrix equal to cosine lat of the boxes with data #and to zero for the boxes of missing data -999.9 for(j in 3:1647) { for (i in 1:2592) {if(temp[i,j]> -290.0) {areaw[i,j]=veca[i]} } } #Compute an area-weighted temperature data matrix and its average: #area-weight data matrix’s first two columns as lat-lon tempw=areaw*temp tempw[,1:2]=temp[,1:2] #create monthly global average vector for 1645 months #Jan 1880- Jan 2017 avev=colSums(tempw[,3:1647])/colSums(areaw[,3:1647]) #Plot time series of the global ave temp timemo=seq(1880,2017,length=1645) plot(timemo,avev,type="l", cex.lab=1.4, xlab="Year", ylab="Temperature anomaly [oC]", main="Area-weighted global average of monthly SAT anomalies: Jan 1880-Jan 2017") abline(lm(avev ~ timemo),col="blue",lwd=2) text(1930,0.7, "Linear trend: 0.69 [oC] per century", cex=1.4, col="blue") #Plot the data coverage time series #Plot this time series motime=seq(1880, 2017, length=1645) rcover=100*colSums(areaw[,3:1647])/sum(veca) plot(motime,rcover,type="l",ylim=c(0,100), main="NOAAGlobalTemp Data Coverage: Jan 1880-Jan 2017", xlab="Year",ylab="Percent area covered [%]") #Work with the NOAAGlobalTemp time series data #Download the NCEI spatial average time series of monthly data #https://www1.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/ #timeseries/aravg.mon.land_ocean.90S.90N.v4.0.1.201702.asc #setwd("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch12-13-Rgraphics") #aveNCEI <- read.table("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch12-13-RGraphics/aravg.mon.land_ocean.90S.90N.v4.0.1.201703.asc.txt", # header=FALSE) aveNCEI <- read.table("http://shen.sdsu.edu/data/aravg.mon.land_ocean.90S.90N.v4.0.1.201703.asc.txt", header=FALSE) dim(aveNCEI) #Jan 1880-Feb 2017 #an extra month to be deleted #[1] 1647 10 avediff<-avev-aveNCEI[1:1645,3] par(mar=c(4,5,2,1)) plot(timemo,avediff,type="l", cex.lab=1.4, xlab="Year", ylab="Diffences [oC]", main="Difference of R average minus NCEI average of global temp") #Plot the annual time series of the NOAAGlobalTemp plot.new() avem = matrix(avev[1:1644], ncol=12, byrow=TRUE) #compute annual average annv=rowMeans(avem) #Plot the annual mean global average temp timeyr<-seq(1880, 2016) plot(timeyr,annv,type="s", cex.lab=1.4, lwd=2, xlab="Year", ylab="Temperature anomaly [oC]", main="Area-weighted global average of annual SAT anomalies: 1880-2016") abline(lm(annv ~ timeyr),col="blue",lwd=2) text(1940,0.4, "Linear trend: 0.69 [oC] per century", cex=1.4, col="blue") text(1900,0.07, "Base line",cex=1.4, col="red") lines(timeyr,rep(0,137), type="l",col="red") #One can compare with the NOAAGlobalTemp annual time series #There are some small differences aveannNCEI <- read.table("http://shen.sdsu.edu/data/aravg.ann.land_ocean.90S.90N.v4.0.1.201703.asc.txt", header=FALSE) dim(aveannNCEI) #[1] 138 6 diff2=annv-aveannNCEI[1:137,2] range(diff2) #[1] -0.016094969 0.007508731 #Nonlinear trend from polynomial fitting #Polynomial fitting to the global average annual mean #poly9<-lm(annv ~ poly(timeyr,9, raw= TRUE)) #raw=TRUE means regular polynomial a0+a1x^2+..., non-orthogonal polyor9<-lm(annv ~ poly(timeyr,9, raw= FALSE)) polyor20<-lm(annv ~ poly(timeyr,20, raw= FALSE)) #raw=FALSE means orthongonal polynomial of 9th order #Orthogonal polynomial fitting is usually better plot(timeyr,annv,type="s", cex.lab=1.4, lwd=2, xlab="Year", ylab="Temperature anomaly [oC]", main="Annual SAT time series and its orthogonal polynomial fits: 1880-2016") lines(timeyr,predict(polyor9),col="blue", lwd=3) legend(1880, 0.6, col=c("blue"),lty=1,lwd=2.0, legend=c("9th order orthogonal polynomial fit"), bty="n",text.font=2,cex=1.5) lines(timeyr,predict(polyor20),col="red", lwd=3) legend(1880, 0.7, col=c("red"),lty=1,lwd=2.0, legend=c("20th order orthogonal polynomial fit"), bty="n",text.font=2,cex=1.5) #Non-parametric smooth by LOWESS fitting scatter.smooth(annv ~ timeyr, span=2/18, cex=0.6) #Compute and plot the trend of gridded data #Compute the trend for each box for the 20th century timemo1=seq(1900,2000, len=1200) temp1=temp temp1[temp1 < -490.00] <- NA trendgl=rep(0,2592) for (i in 1:2592){ if(is.na(temp1[i,243])==FALSE & is.na(temp1[i,1442])==FALSE) {trendgl[i]=lm(temp1[i,243: 1442] ~ timemo1, na.action=na.omit)$coefficients[2]} else {trendgl[i]=NA} } library(maps) Lat= seq(-87.5, 87.5, length=36) Lon=seq(2.5, 357.5, length=72) mapmat=matrix(100*trendgl,nrow=72) mapmat=pmax(pmin(mapmat,2),-2) #Matrix flipping is not needed since the data goes from 2.5 to 375.5 plot.new() par(mar=c(4,5,3,0)) int=seq(-2,2,length.out=21) 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="Jan 1900-Dec 1999 temperature trends: [oC/century]", 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)}) #Compute trend for each box for the 20th century #Version 2: Allow 2/3 of data, i.e., 1/3 missing #Compute the trend timemo1=seq(1900,2000, len=1200) temp1=temp[,243:1442] temp1[temp1 < -490.00] <- NA temptf=is.na(temp1) bt=et=rep(0,2592) for (i in 1:2592) { if (length(which(temptf[i,]==FALSE)) !=0) { bt[i]=min(which(temptf[i,]==FALSE)) et[i]=max(which(temptf[i,]==FALSE)) } } trend20c=rep(0,2592) for (i in 1:2592){ if(et[i]-bt[i] > 800) {trend20c[i]=lm(temp1[i,bt[i]:et[i]] ~ seq(bt[i],et[i]), na.action=na.omit)$coefficients[2]} else {trend20c[i]=NA} } #plot the 20C V2 trend map plot.new() #par(mar=c(4,5,3,0)) mapmat=matrix(120*trend20c,nrow=72) mapmat=pmax(pmin(mapmat,0.2),-0.2) int=seq(-0.2,0.2,length.out=41) rgb.palette=colorRampPalette(c('black','blue', 'darkgreen','green', 'yellow','pink','red','maroon'),interpolate='spline') filled.contour(Lon, Lat, mapmat, color.palette=rgb.palette, levels=int, plot.title=title(main="Jan 1900-Dec 1999 temperature trends: [oC/decade]", 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)}) #Trend from 1976-2016 timemo2=seq(1976,2017, len=492) temp1=temp temp1[temp1 < -490.00] <- NA trend7616=rep(0,2592) for (i in 1:2592){ if(is.na(temp1[i,1155])==FALSE & is.na(temp1[i,1646])==FALSE) {trend7616[i]=lm(temp1[i,1155: 1646] ~ timemo2, na.action=na.omit)$coefficients[2]} else {trend7616[i]=NA} } #Plot US temp and prec times series on the same figure plot.new() Time <- 2001:2010 Tmean <- c(12.06, 11.78,11.81,11.72,12.02,12.36,12.03,11.27,11.33,11.66) Prec <- c(737.11,737.87,774.95,844.55,764.03,757.43,741.17,793.50,820.42,796.80) plot(Time,Tmean,type="o",col="red",xlab="Year", ylab="Tmean [dec C]",lwd=1.5, main="Contiguous U.S. Annual Mean Temperature and Total Precipitation") legend(2000.5,12.42, col=c("red"),lty=1,lwd=2.0, legend=c("Tmean"),bty="n",text.font=2,cex=1.0) #Allows a figure to be overlaid on the first plot par(new=TRUE) plot(Time, Prec,type="o",col="blue",lwd=1.5,axes=FALSE,xlab="",ylab="") legend(2000.5,839, col=c("blue"),lty=1,lwd=2.0, legend=c("Prec"),bty="n",text.font=2,cex=1.0) #Suppress the axes and assign the y-axis to side 4 axis(4) mtext("Precipitation [mm]",side=4,line=3) #legend("topleft",col=c("red","blue"),lty=1,legend=c("Tmean","Prec"),cex=0.6) #Plotting two legends at the same time is difficult to adjust the font size #because of different scales #Margins, math symbol, and figure setups plot.new() par(mar=c(6,4,3,3)) x<-0.25*(-30:30) y<-sin(x) x1<-x[which(sin(x) >=0)] y1<-sin(x1) x2<-x[which(sin(x) < 0)] y2<-sin(x2) plot(x1,y1,xaxt="n", xlab="",ylab="",lty=1,type="h", lwd=3, tck=-0.02, ylim=c(-1,1), col="red", col.lab="purple",cex.axis=1.4) lines(x2,y2,xaxt="n", xlab="",ylab="",lty=3,type="h", col="blue",lwd=8, tck=-0.02) axis(1, at=seq(-6,6,2),line=3, cex.axis=1.8) axis(4, at=seq(-1,1,0.5), lab=c("A", "B", "C", "D","E"), cex.axis=2,las=2) text(0,0.7,font=3,cex=6, "Sine waves", col="darkgreen") #Itatlic font size 2 mtext(side=2,line=2, expression(y==sin(theta-hat(phi))),cex=1.5, col="blue") mtext(font=2,"Text outside of the figure on side 3",side=3,line=1, cex=1.5)#Bold font mtext(font=1, side=1,line=1, expression(paste("Angle in radian: ", theta-phi[0])),cex=1.5, col="red") #change the font sizes for axis labels, main title, and sub-title par(mar=c(8,6,3,2)) par(mgp=c(2.5,1,0)) plot(1:200/20, rnorm(200),sub="Sub-title: 200 random values", xlab= "Time", ylab="Random values", main="Normal random values", cex.lab=1.5, cex.axis=2, cex.main=2.5, cex.sub=2.0) #A fancy IPCC plot of the NOAAGlobalTemp time series NOAATemp=read.table("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch12-13-RGraphics/aravg.ann.land_ocean.90S.90N.v4.0.1.201703.asc.txt", header=FALSE) plot.new() 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 US temp and prec times series on the same figure par(mfrow=c(2,1)) par(mar=c(0,5,3,1)) #Zero space between (a) and (b) Time <- 2001:2010 Tmean <- c(12.06, 11.78,11.81,11.72,12.02,12.36,12.03,11.27,11.33,11.66) Prec <- c(737.11,737.87,774.95,844.55,764.03,757.43,741.17,793.50,820.42,796.80) plot(Time,Tmean,type="o",col="red",xaxt="n", xlab="",ylab="Tmean [dec C]") text(2006, 12,font=2,"US Annual Mean Temperature", cex=1.5) text(2001.5,12.25,"(a)") #Plot the panel on row 2 par(mar=c(3,5,0,1)) plot(Time, Prec,type="o",col="blue",xlab="Time",ylab="Prec [mm]") text(2006, 800, font=2, "US Annual Total Precipitation", cex=1.5) text(2001.5,840,"(b)") #Stack multiple panels together as a single figure plot.new() layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE), widths=c(3,3), heights=c(2,2)) plot(sin,xlim=c(0,20)) plot(sin,xlim=c(0,10)) plot(sin,xlim=c(10,20)) #Plot a 5-by-5 grid global map of standard normal random values library(maps) plot.new() #Step 1: Generate a 5-by-5 grid (pole-to-pole, lon 0 to 355) Lat<-seq(-90,90,length=37) #Must increasing Lon<-seq(0,355,length=72) #Must increasing #Generate the random values mapdat<-matrix(rnorm(72*37),nrow=72) #The matrix uses lon as row going and lat as column #Each row includes data from south to north #Define color int=seq(-3,3,length.out=81) rgb.palette=colorRampPalette(c('black','purple','blue','white', 'green', 'yellow','pink','red','maroon'), interpolate='spline') #Plot the values on the world map filled.contour(Lon, Lat, mapdat, color.palette=rgb.palette, levels=int, plot.title=title(xlab="Longitude", ylab="Latitude", main="Standard Normal Random Values on a World Map: 5 Lat-Lon Grid"), plot.axes={ axis(1); axis(2);map('world2', add=TRUE);grid()} ) #filled.contour() is a contour plot on an x-y grid. #Background maps are added later in plot.axes={} #axis(1) means ticks on the lower side #axis(2) means ticks on the left side #Save image with width=800, maintain aspect ratio #Plot a regional map of a random dataset on a map #Plot a 5-by-5 grid regional map to cover USA and Canada Lat3<-seq(10,70,length=13) Lon3<-seq(230,295,length=14) mapdat<-matrix(rnorm(13*14),nrow=14) int=seq(-3,3,length.out=81) rgb.palette=colorRampPalette(c('black','purple','blue','white', 'green', 'yellow','pink','red','maroon'), interpolate='spline') filled.contour(Lon3, Lat3, mapdat, color.palette=rgb.palette, levels=int, plot.title=title(main="Standard Normal Random Values on a World Map: 5-deg Lat-Lon Grid", xlab="Lon", ylab="Lat"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}) #R plot of NCEP/NCAR Reanalysis PSD monthly temp data .nc file #http://www.esrl.noaa.gov/psd/data/gridded/data.ncep. #reanalysis.derived.surface.html rm(list=ls(all=TRUE)) setwd("/Users/sshen/Desktop/Papers/KarlTom/Recon2016/Test-with-Gregori-prec-data") # Download netCDF file # Library install.packages("ncdf") library(ncdf4) # 4 dimensions: lon,lat,level,time nc=ncdf4::nc_open("air.mon.mean.nc") nc nc$dim$lon$vals # output values 0.0->357.5 nc$dim$lat$vals #output values 90->-90 nc$dim$time$vals #nc$dim$time$units #nc$dim$level$vals Lon <- ncvar_get(nc, "lon") Lat1 <- ncvar_get(nc, "lat") Time<- ncvar_get(nc, "time") head(Time) #[1] 65378 65409 65437 65468 65498 65529 #install.packages("ncdf") library(chron) month.day.year(1297320/24,c(month = 1, day = 1, year = 1800)) #1948-01-01 precnc<- ncvar_get(nc, "air") dim(precnc) #[1] 144 73 826, i.e., 826 months=1948-01 to 2016-10, 68 years 10 mons #plot a 90S-90N precip along a meridional line at 160E over Pacific dev.off() plot(seq(-90,90,length=73),precnc[15,,1], type="l", xlab="Lat", ylab="Temp [oC]", main="90S-90N temperature [mm/day] along a meridional line at 35E: Jan 1948", lwd=3) #Compute and plot NCEP Reanalysis climatology #Compute and plot climatology and standard deviation Jan 1948-Dec 2015 library(maps) climmat=matrix(0,nrow=144,ncol=73) sdmat=matrix(0,nrow=144,ncol=73) Jmon<-12*seq(0,67,1) for (i in 1:144){ for (j in 1:73) {climmat[i,j]=mean(precnc[i,j,Jmon]); sdmat[i,j]=sd(precnc[i,j,]) } } mapmat=climmat #Note that R requires coordinates increasing from south to north -90->90 #and from west to east from 0->360. We have to make Lat and Lon this way. #Correpondingly, we have to flip the data matrix left to right according to #the data matrix precnc[i,j,]: 360 (i.e. 180W) lon and from North Pole #and South Pole, then lon 178.75W, 176.75W, ..., 0E. This puts Greenwich #at center, China on the right, and USA on the left. But our map should #have Pacific at center, and USA on teh right. Thus, we make a flip. Lat=-Lat1 mapmat= mapmat[,length(mapmat[1,]):1]#Matrix flip around a column #mapmat= t(apply(t(mapmat),2,rev)) int=seq(-50,50,length.out=81) rgb.palette=colorRampPalette(c('black','blue','darkgreen','green', 'white','yellow','pink','red','maroon'),interpolate='spline') filled.contour(Lon, Lat, mapmat, color.palette=rgb.palette, levels=int, plot.title=title(main="NCEP RA 1948-2015 January climatology [deg C]", xlab="Longitude",ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) #plot standard deviation plot.new() par(mgp=c(2,1,0)) par(mar=c(3,3,2,2)) mapmat= sdmat[,seq(length(sdmat[1,]),1)] int=seq(0,20,length.out=81) rgb.palette=colorRampPalette(c('black','blue', 'green','yellow','pink','red','maroon'), interpolate='spline') filled.contour(Lon, Lat, mapmat, color.palette=rgb.palette, levels=int, plot.title=title(main="NCEP 1948-2015 Jan SAT RA Standard Deviation [deg C]", xlab="Longitude", ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) #Plot the January 1983 temperature using the above setup mapmat83J=precnc[,,421] mapmat83J= mapmat83J[,length(mapmat83J[1,]):1] int=seq(-50,50,length.out=81) rgb.palette=colorRampPalette(c('black','blue','darkgreen', 'green', 'white','yellow','pink','red','maroon'),interpolate='spline') filled.contour(Lon, Lat, mapmat83J, color.palette=rgb.palette, levels=int, plot.title=title(main="January 1983 surface air temperature [deg C]", xlab="Longitude",ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) #Plot the January 1983 temperature anomaly from NCEP data plot.new() anomat=precnc[,,421]-climmat anomat=pmin(anomat,6) anomat=pmax(anomat,-6) anomat= anomat[,seq(length(anomat[1,]),1)] int=seq(-6,6,length.out=81) rgb.palette=colorRampPalette(c('black','blue','darkgreen','green', 'white','yellow','pink','red','maroon'),interpolate='spline') filled.contour(Lon, Lat, anomat, color.palette=rgb.palette, levels=int, plot.title=title(main="January 1983 surface air temperature anomaly [deg C]", xlab="Longitude",ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) #Zoom in to a specific lat-lon region: Pacific int=seq(-6,6,length.out=81) rgb.palette=colorRampPalette(c('black','blue','darkgreen','green', 'white','yellow','pink','red','maroon'), interpolate='spline') filled.contour(Lon, Lat, anomat, xlim=c(100,300),ylim=c(-50,50),zlim=c(-6,6), color.palette=rgb.palette, levels=int, plot.title=title( main="January 1983 surface air temperature anomaly [deg C]", xlab="Longitude",ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) #Wind directions due to the balance between PGF and Coriolis force #using an arrow plot for vector fields on a map library(fields) library(maps) library(mapproj) lat<-rep(seq(-75,75,len=6),12) lon<-rep(seq(-165,165,len=12),each=6) x<-lon y<-lat #u<- rep(c(-1,1,-1,-1,1,-1), each=12) #v<- rep(c(1,-1,1,-1,1,-1), each=12) u<- rep(c(-1,1,-1,-1,1,-1), 12) v<- rep(c(1,-1,1,-1,1,-1), 12) wmap<-map(database="world", boundary=TRUE, interior=TRUE) grid(nx=12,ny=6) #map.grid(wmap,col=3,nx=12,ny=6,label=TRUE,lty=2) points(lon, lat,pch=16,cex=0.8) arrow.plot(lon,lat,u,v, arrow.ex=.08, length=.08, col='blue', lwd=2) box() axis(1, at=seq(-165,135,60), lab=c("165W","105W","45W","15E","75E","135E"), col.axis="black",tck = -0.05, las=1, line=-0.9,lwd=0) axis(1, at=seq(-165,135,60), col.axis="black",tck = 0.05, las=1, labels = NA) axis(2, at=seq(-75,75,30),lab=c("75S","45S","15S","15N","45N","75N"), col.axis="black", tck = -0.05, las=2, line=-0.9,lwd=0) axis(2, at=seq(-75,75,30), col.axis="black", tck = 0.05, las=1, labels = NA) text(30, 0, "Intertropical Convergence Zone (ITCZ)", col="red") text(75, 30, "Subtropical High", col="red") text(75, -30, "Subtropical High", col="red") mtext(side=3, "Polar High", col="red", line=0.0) #ggplot example #ggplot for USA States library(ggplot2) states <- map_data("state") p<- ggplot(data = states) + geom_polygon(aes(x = long, y = lat, fill = region, group = group), color = "white") + coord_fixed(1.3) #guides(fill=TRUE) # This leaves off the color legend on the right p<- p + xlab("Longitude")+ ylab("Latitude") p<- p + ggtitle("Color Map of the 48 Lower States") p #Plot the wind field over the ocean #Ref: https://rpubs.com/alobo/vectorplot #Agustin.Lobo@ictja.csic.es #20140428 library(ncdf4) library(chron) library(RColorBrewer) library(lattice) download.file("ftp://eclipse.ncdc.noaa.gov/pub/seawinds/SI/uv/clm/uvclm95to05.nc", "uvclm95to05.nc", method = "curl") mincwind <- nc_open("uvclm95to05.nc") dim(mincwind) #print.nc(mincwind) u <- ncvar_get(mincwind, "u") class(u) dim(u) v <- ncvar_get(mincwind, "v") class(v) dim(v) u9 <- raster(t(u[, , 9])[ncol(u):1, ]) v9 <- raster(t(v[, , 9])[ncol(v):1, ]) filled.contour(u[, , 9]) filled.contour(u[, , 9],color.palette = heat.colors) filled.contour(u[, , 9],color.palette = colorRampPalette(c("red", "white", "blue"))) contourplot(u[, , 9]) install.packages("raster") library(raster) library(sp) library(rgdal) u9 <- raster(t(u[, , 9])[ncol(u):1, ]) v9 <- raster(t(v[, , 9])[ncol(v):1, ]) w <- brick(u9, v9) wlon <- ncvar_get(mincwind, "lon") wlat <- ncvar_get(mincwind, "lat") range(wlon) range(wlat) projection(w) <- CRS("+init=epsg:4326") extent(w) <- c(min(wlon), max(wlon), min(wlat), max(wlat)) w plot(w[[1]]) plot(w[[2]]) install.packages("rasterVis") install.packages("latticeExtra") library(latticeEtra) library(rasterVis) vectorplot(w * 10, isField = "dXY", region = FALSE, margin = FALSE, narrows = 10000) slope <- sqrt(w[[1]]^2 + w[[2]]^2) aspect <- atan2(w[[1]], w[[2]]) vectorplot(w * 10, isField = "dXY", region = slope, margin = FALSE, par.settings=BuRdTheme, narrows = 10000, at = 0:10) vectorplot(stack(slope * 10, aspect), isField = TRUE, region = FALSE, margin = FALSE) #Generate 3D data x<-seq(0, 2*pi, len=100) y<-seq(0, 2*pi, len=100) #t<-seq(1,20,by=1) mydat<-array(0,dim=c(100,100,10)) for(t in 1:10){z<-function(x,y){z=sin(t)*(1/pi)*sin(x)*sin(y) + exp(-0.3*t)*(1/pi)*sin(8*x)*sin(8*y)} mydat[,,t]=outer(x,y,z) } #Plot the original z(x,y,t) waves for a given t plot.new() int=seq(-0.5,0.5,length.out=61) filled.contour(x,y,mydat[,,1], color.palette =rainbow, levels=int, plot.title=title(main="Orignal field at t=10", xlab="x", ylab="y", cex.lab=1.0), key.title = title(main = "Scale"), plot.axes = {axis(1,seq(0,3*pi, by = 1), cex=1.0) axis(2,seq(0, 2*pi, by = 1), cex=1.0)} ) #3D array to 2D space-time matrix da1<- matrix(0,nrow=length(x)*length(y),ncol=10) for (i in 1:10) {da1[,i]=c(t(mydat[,,i]))} #SVD for EOFs da2<-svd(da1) uda2<-da2$u #EOFs vda2<-da2$v #PCs dda2<-da2$d #Energy for the corresponding EOFs dda2 #[1] 3.589047e+01 1.596154e+01 7.764115e-14 6.081008e-14 dda2^2/sum(dda2^2) # [1] 8.348751e-01 1.651249e-01 1.939476e-29 #The first mode counts for variance 83% variance #The second 17%, and almost zero for higher modes #Plot EOF par(mgp=c(2,1,0)) filled.contour(x,y,matrix(-uda2[,1],nrow=100), color.palette =rainbow, plot.title=title(main="SVD Mode 1: EOF1", xlab="x", ylab="y", cex.lab=1.0), key.title = title(main = "Scale"), plot.axes = {axis(1,seq(0,2*pi, by = 1), cex=1.0) axis(2,seq(0, 2*pi, by = 1), cex=1.0)}) #Accurate spatial patterns from functions that generate data z1 <- function(x,y){(1/pi)*sin(x)*sin(y)} z2 <- function(x,y){(1/pi)*sin(5*x)*sin(5*y)} fcn1<-outer(x,y,z1) fcn2<-outer(x,y,z2) par(mgp=c(2,1,0)) filled.contour(x,y,fcn1, color.palette =rainbow, plot.title=title(main="Accurate Mode 1", xlab="x", ylab="y", cex.lab=1.0), key.title = title(main = "Scale"), plot.axes = {axis(1,seq(0,3*pi, by = 1), cex=1.0) axis(2,seq(0, 2*pi, by = 1), cex=1.0)} ) #Plot PCs and coefficients of the functional patterns t=1:10 plot(1:10, vda2[,1],type="o", ylim=c(-1,1), lwd=2, ylab="PC or Coefficient", xlab="Time", main="SVD PCs vs. Accurate Temporal Coefficients") legend(0.5,1.15, lty=1, legend=c("PC1: 83% variance"), bty="n",col=c("black")) lines(1:10, vda2[,2],type="o", col="red", lwd=2) legend(0.5,1.0, lty=1, legend=c("PC2: 17% varance"), col="red", bty="n",text.col=c("red")) lines(t, -sin(t), col="blue", type="o") legend(0.5,0.85, lty=1, legend=c("Original mode 1 coefficient -sin(t): 91% variance"), col="blue", bty="n",text.col="blue") lines(t, -exp(-0.3*t), type="o",col="purple") legend(0.5,0.70, lty=1, legend=c("Original mode 2 coefficient -exp(-0.3t): 9% variance"), col="purple", bty="n",text.col="purple") #Orthonormal verification t=1:10 t(-sin(t))%*%(-exp(-0.3*t)) v1=var(-sin(t)) v2=var(-exp(-0.3*t)) v1/(v1+v2) #[1] 0.9098719 #Verify orthogonality of PCs t(vda2[,1])%*%vda2[,2] # [1,] -5.551115e-17 t=1:10 t(-sin(t))%*%(-exp(-0.3*t)) #[1,] 0.8625048 #Data reconstruction from 2 modes B<-uda2[,1:2]%*%diag(dda2)[1:2,1:2]%*%t(vda2[,1:2]) B1<-uda2%*%diag(dda2)%*%t(vda2) #Plot the reconstructed data plot.new() filled.contour(x,y,matrix(B[,5],nrow=100), color.palette =rainbow, plot.title=title(main="2-mode SVD reduced field t=5", xlab="x", ylab="y", cex.lab=1.0), key.title = title(main = "Scale"), plot.axes = {axis(1,seq(0,3*pi, by = 1), cex=1.0) axis(2,seq(0, 2*pi, by = 1), cex=1.0)}) #Download and visualize the NCEP temperature data # Download netCDF file # Library install.packages("ncdf") library(ncdf4) # 4 dimensions: lon,lat,level,time nc=ncdf4::nc_open("/Users/sshen/Desktop/Papers/KarlTom/Recon2016/Test-with-Gregori-prec-data/air.mon.mean.nc") nc nc$dim$lon$vals #output lon values 0.0->357.5 nc$dim$lat$vals #output lat values 90->-90 nc$dim$time$vals #output time values in GMT hours: 1297320, 1298064 nc$dim$time$units #[1] "hours since 1800-01-01 00:00:0.0" #nc$dim$level$vals Lon <- ncvar_get(nc, "lon") Lat <- ncvar_get(nc, "lat") Time<- ncvar_get(nc, "time") #Time is the same as nc$dim$time$vals head(Time) #[1] 1297320 1298064 1298760 1299504 1300224 1300968 library(chron) Tymd<-month.day.year(Time[1]/24,c(month = 1, day = 1, year = 1800)) #c(month = 1, day = 1, year = 1800) is the reference time Tymd #$month #[1] 1 #$day #[1] 1 #$year #[1] 1948 #1948-01-01 precnc<- ncvar_get(nc, "air") dim(precnc) #[1] 144 73 826, i.e., 826 months=1948-01 to 2016-10, 68 years 10 mons #Plot the data along a meridional line #plot a 90S-90N temp along a meridional line at 180E plot(seq(-90,90,length=73),precnc[72,,1], type="l", xlab="Lat", ylab="Temp [oC]", main="90S-90N temperature [mm/day] along a meridional line at 35E: Jan 1948") #Convert NCEp Reanalysis into a space-time dataset #Write the data as space-time matrix with a header precst=matrix(0,nrow=10512,ncol=826) temp=as.vector(precnc[,,1]) head(temp) for (i in 1:826) {precst[,i]=as.vector(precnc[ , , i])} dim(precst) #[1] 10512 826 #Build lat and lon for 10512 spatial positions using rep LAT=rep(Lat, 144) LON=rep(Lon[1],73) for (i in 2:144){LON=c(LON,rep(Lon[i],73))} gpcpst=cbind(LAT, LON, precst) dim(gpcpst) #[1] 10512 828 #The first two columns are lat and lon. 826 mons: 1948.01-2016.10 #Convert the Julian day and hour into calendar mons for time tm=month.day.year(Time/24, c(month = 1, day = 1, year = 1800)) tm1=paste(tm$year,"-",tm$month) #tm1=data.frame(tm1) tm2=c("Lat","Lon",tm1) colnames(gpcpst) <- tm2 setwd("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch12-13-RGraphics") #setwd routes the desired csv file to a given directory write.csv(gpcpst,file="ncepJan1948_Oct2016.csv") #Compute the January climatology and standard deviation is below. monJ=seq(1,816,12) #Select each January gpcpdat=gpcpst[,3:818] gpcpJ=gpcpdat[,monJ] climJ<-rowMeans(gpcpJ) library(matrixStats)# rowSds command is in the matrixStats package sdJ<-rowSds(gpcpJ) #Plot climatology and standard deviation #Plot Jan climatology library(maps) Lat1=seq(90,-90, len=73) Lat=-Lat1 mapmat=matrix(climJ,nrow=144) mapmat= mapmat[,seq(length(mapmat[1,]),1)] plot.new() int=seq(-50,50,length.out=81) rgb.palette=colorRampPalette(c('black','blue', 'darkgreen','green', 'white','yellow','pink','red','maroon'),interpolate='spline') filled.contour(Lon, Lat, mapmat, color.palette=rgb.palette, levels=int, plot.title=title(main="NCEP Jan SAT RA 1948-2015 climatology [deg C]", xlab="Longitude", ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) #--------------------- #Plot Jan Standard Deviation Lat=-Lat1 mapmat=matrix(sdJ,nrow=144) mapmat= mapmat[,seq(length(mapmat[1,]),1)] plot.new() mapmat=pmin(mapmat,5) #Compress the values >5 to 5 int=seq(0,5,length.out=81) rgb.palette=colorRampPalette(c('black','blue', 'green', 'yellow','pink','red','maroon'),interpolate='spline') filled.contour(Lon, Lat, mapmat, color.palette=rgb.palette, levels=int, plot.title=title(main="NCEP Jan SAT RA 1948-2015 Standard Deviation [deg C]", xlab="Longitude", ylab="Latitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="[oC]")) #Compute the Jan EOFs monJ=seq(1,816,12) gpcpdat=gpcpst[,3:818] gpcpJ=gpcpdat[,monJ] climJ<-rowMeans(gpcpJ) library(matrixStats) sdJ<-rowSds(gpcpJ) anomJ=(gpcpdat[,monJ]-climJ)/sdJ #standardized anomalies anomAW=sqrt(cos(gpcpst[,1]*pi/180))*anomJ #Area weighted anormalies svdJ=svd(anomAW) #execute SVD #plot eigenvalues par(mar=c(3,4,2,4)) plot(100*(svdJ$d)^2/sum((svdJ$d)^2), type="o", ylab="Percentage of variance [%]", xlab="Mode number", main="Eigenvalues of covariance matrix") legend(20,5, col=c("black"),lty=1,lwd=2.0, legend=c("Percentange variance"),bty="n", text.font=2,cex=1.0, text.col="black") par(new=TRUE) plot(cumsum(100*(svdJ$d)^2/sum((svdJ$d)^2)),type="o", col="blue",lwd=1.5,axes=FALSE,xlab="",ylab="") legend(20,50, col=c("blue"),lty=1,lwd=2.0, legend=c("Cumulative variance"),bty="n", text.font=2,cex=1.0, text.col="blue") axis(4) mtext("Cumulative variance [%]",side=4,line=2) #plot EOF1: The physical EOF= eigenvector divided by area factor mapmat=matrix(svdJ$u[,1]/sqrt(cos(gpcpst[,1]*pi/180)),nrow=144) rgb.palette=colorRampPalette(c('blue','green','white', 'yellow','red'),interpolate='spline') int=seq(-0.04,0.04,length.out=61) mapmat=mapmat[, seq(length(mapmat[1,]),1)] filled.contour(Lon, Lat, -mapmat, color.palette=rgb.palette, levels=int, plot.title=title(main="January EOF1 from 1948-2016 NCEP Temp Data"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="Scale")) # #plot PC1 pcdat<-svdJ$v[,1] Time<-seq(1948,2015) plot(Time, -pcdat, type="o", main="PC1 of NCEP RA Jan SAT: 1948-2015", xlab="Year",ylab="PC values", lwd=2, ylim=c(-0.3,0.3)) #EOF from de-trended data monJ=seq(1,816,12) gpcpdat=gpcpst[,3:818] gpcpJ=gpcpdat[,monJ] climJ<-rowMeans(gpcpJ) library(matrixStats) sdJ<-rowSds(gpcpJ) anomJ=(gpcpdat[,monJ]-climJ)/sdJ trendM<-matrix(0,nrow=10512, ncol=68)#trend field matrix trendV<-rep(0,len=10512)#trend for each grid box: a vector for (i in 1:10512) { trendM[i,] = (lm(anomJ[i,] ~ Time))$fitted.values trendV[i]<-lm(anomJ[i,] ~ Time)$coefficients[2] } dtanomJ = anomJ - trendM dim(dtanomJ) dtanomAW=sqrt(cos(gpcpst[,1]*pi/180))*dtanomJ svdJ=svd(dtanomAW) #Plot the area-weighted global average Jan temp from 1948-2015 #Begin from the space-time data matrix gpcpst[,1] vArea=cos(gpcpst[,1]*pi/180) anomA=vArea*anomJ dim(anomA) JanSAT<-colSums(anomA)/sum(vArea) plot(Time, JanSAT, type="o", lwd=2, main="Global Average Jan SAT Anomalies from NCEP RA", xlab="Year",ylab="Temperature [oC]") regSAT<-lm(JanSAT ~ Time) #0.61oC/100a trend abline(regSAT, col="red", lwd=4) text(1965,0.4,"Linear trend 0.61oC/100a", col="red", cex=1.3) #plot the trend of Jan SAT non-standardized anomaly data #Begin with the space-time data matrix monJ=seq(1,816,12) gpcpdat=gpcpst[,3:818] gpcpJ=gpcpdat[,monJ] plot(gpcpJ[,23]) climJ<-rowMeans(gpcpJ) anomJ=(gpcpdat[,monJ]-climJ) trendV<-rep(0,len=10512)#trend for each grid box: a vector for (i in 1:10512) { trendV[i]<-lm(anomJ[i,] ~ Time)$coefficients[2] } mapmat1=matrix(10*trendV,nrow=144) mapv1=pmin(mapmat1,1) #Compress the values >5 to 5 mapmat=pmax(mapv1,-1) #compress the values <-5 t -5 rgb.palette=colorRampPalette(c('blue','green','white', 'yellow','red'), interpolate='spline') int=seq(-1,1,length.out=61) mapmat=mapmat[, seq(length(mapmat[1,]),1)] filled.contour(Lon, Lat, mapmat, color.palette=rgb.palette, levels=int, plot.title=title( main="Trend of the NCEP RA1 Jan 1948-2015 Anom Temp", xlab="Latitude",ylab="Longitude"), plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, key.title=title(main="oC/10a")) #EOF analysis of the Tahiti and Darwin standardized SLP data #SOI analysis # Read the txt data #Pta<-read.table("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch5-SOI/PSTANDtahiti", header=F) Pta<-read.table("http://shen.sdsu.edu/data/PSTANDtahiti", header=F) # Remove the first column that is the year ptamon<-Pta[,seq(2,13)] #Convert the matrix into a vector according to mon: Jan 1951, Feb 1951, ..., Dec 2015 ptamonv<-c(t(ptamon)) xtime<-seq(1951, 2016-1/12, 1/12) # Plot the Tahiti standardized SLP anomalies plot(xtime, ptamonv,type="l",xlab="Year",ylab="Presure", main="Standardized Tahiti SLP Anomalies", col="red", xlim=range(xtime), ylim=range(ptamonv)) # Do the same for Darwin SLP #Pda<-read.table("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch5-SOI/PSTANDdarwin.txt", header=F) Pda<-read.table("http://shen.sdsu.edu/data/PSTANDdarwin.txt", header=F) pdamon<-Pda[,seq(2,13)] pdamonv<-c(t(pdamon)) plot(xtime, pdamonv,type="l",xlab="Year",ylab="Presure", main="Standardized Darwin SLP Anomalies", col="blue", xlim=range(xtime), ylim=range(pdamonv)) #Plot the SOI index plot(xtime, ptamonv-pdamonv,type="l",xlab="Year", ylab="SOI index", col="black",xlim=range(xtime), ylim=c(-4,4), lwd=1) #Add ticks on top edge of the plot box axis (3, at=seq(1951,2015,4), labels=seq(1951,2015,4)) #Add ticks on the right edge of the plot box axis (4, at=seq(-4,4,2), labels=seq(-4,4,2)) # If put a line on a plot, use the command below lines(xtime,ptamonv-pdamonv,col="black", lwd=1) #Plot cumulative SOI: CSOI cnegsoi<--cumsum(ptamonv-pdamonv) plot(xtime, cnegsoi,type="l",xlab="Year",ylab="Negative CSOI index", col="black",xlim=range(xtime), ylim=range(cnegsoi), lwd=1) #Space-time SLP/SOI data ptada <-cbind(ptamonv,pdamonv) ptada<-t(ptada) svdptd<-svd(ptada) recontd=svdptd$u%*%diag(svdptd$d[1:2])%*%t(svdptd$v) #Plot WSOI1 D=diag(svdptd$d) U=svdptd$u V=svdptd$v xtime<-seq(1951, 2016-1/12, 1/12) wsoi1=D[1,1]*t(V)[1,] plot(xtime, wsoi1,type="l",xlab="Year",ylab="Weighted SOI 1", col="black",xlim=range(xtime), ylim=range(wsoi1), lwd=1) axis (3, at=seq(1951,2015,4), labels=seq(1951,2015,4)) #Plot WSOI2 wsoi2=D[2,2]*t(V)[2,] plot(xtime, wsoi2,type="l",xlab="Year",ylab="Weighted SOI 2", col="black",xlim=range(xtime), ylim=c(-2,2), lwd=1) axis (3, at=seq(1951,2015,4), labels=seq(1951,2015,4)) #SVD Analysis for the Darwin-Tahiti Standardized SLP #By Sam Shen, August 2015, Revised May 2017 setwd("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/Ch5-SOI") Pta<-read.table("PSTANDtahiti", header=F) # Remove the first column that is the year ptamon<-Pta[,seq(2,13)] #Convert the matrix into a vector according to mon: Jan 1951, Feb 1951, ..., Dec 2015 ptamonv<-c(t(ptamon)) #Generate time ticks from Jan 1951 to Dec 2015 xtime<-seq(1951, 2016-1/12, 1/12) # Plot the Tahiti standardized SLP anomalies plot(xtime, ptamonv,type="l",xlab="Year",ylab="Presure", main="Standardized Tahiti SLP Anomalies", col="red", xlim=range(xtime), ylim=range(ptamonv)) # Do the same for Darwin SLP Pda<-read.table("PSTANDdarwin.txt", header=F) pdamon<-Pda[,seq(2,13)] pdamonv<-c(t(pdamon)) plot(xtime, pdamonv,type="l",xlab="Year",ylab="Presure", main="Standardized Darwin SLP Anomalies", col="blue", xlim=range(xtime), ylim=range(pdamonv)) #Plot the SOI index SOI <- ptamonv-pdamonv plot(xtime, SOI,type="l",xlab="Year",main="Southern Oscillation Index", ylab="SOI index", col="black",xlim=range(xtime), ylim=c(-6,6), lwd=1) #Add ticks on top edge of the plot box axis (3, at=seq(1951,2015,4), labels=seq(1951,2015,4)) #Add ticks on the right edge of the plot box axis (4, at=seq(-4,4,2), labels=seq(-4,4,2)) lines(xtime, rep(0,length(xtime))) abline(lm(SOI ~ xtime), col="red", lwd=2) #Cumulative negative SOI cnegsoi<--cumsum(SOI) plot(xtime, cnegsoi,type="l",xlab="Year",ylab="Negative CSOI index", col="black",xlim=range(xtime), ylim=range(cnegsoi), lwd=2, main="Cumulative southern oscillation index") #Space-time standardized SLP data matrix ptada <-cbind(ptamonv,pdamonv) ptada <- t(ptada) dim(ptada) #[1] 2 780 #SVD analysis of the space-time data svdptd <- svd(ptada) EOF = svdptd$u EOF # [,1] [,2] #[1,] -0.6146784 0.7887779 #[2,] 0.7887779 0.6146784 #Weighted SOI WSOI=t(EOF)%*%ptada dim(WSOI) #The two weighted SOI indices WSOI1=WSOI[1,] WSOI2=WSOI[2,] plot(xtime, WSOI1,type="l",xlab="Year", main="Weighted Southern Oscillation Indices", ylab="WSOI index", col="black",xlim=range(xtime), ylim=c(-4,4), lwd=1) legend(1945,5.5, lty=1, legend=c("WSOI1: 67% variance"), bty="n",col=c("black"),cex=1.3) lines(xtime, WSOI2,type="l", col="blue") legend(1945,4.7, lty=1, legend=c("WSOI2: 33% variance"), bty="n",col=c("blue"), cex=1.3) #Energy for each mode Energy=svdptd$d Energy #[1] 31.34582 22.25421 #Energy ratio Energy^2/sum(Energy^2) #Principal components PC=svdptd$v PC1=PC[,1] PC2=PC[,2] #Plot the PCs plot(xtime, PC1,type="l",xlab="Year", main="Tahiti-Darwin Principal Components", ylab="PC index", col="black",xlim=range(xtime), ylim=c(-0.15,0.15), lwd=1) legend(1945,0.20, lty=1, legend=c("PC1: 67% variance"), bty="n",col=c("black"),cex=1.3) lines(xtime, PC2,type="l", col="blue") legend(1945,0.17, lty=1, legend=c("PC2: 33% variance"), bty="n",col=c("blue"), cex=1.3) #Plot cumulative PCs plot(xtime, cumsum(PC1),type="l",xlab="Year", main="Cumulative Tahiti-Darwin Principal Components", ylab="Cumulative PC index", col="black",xlim=range(xtime), ylim=c(-6,1),lwd=1) legend(1960,2, lty=1, legend=c("Cumulative PC1: 67% variance"), bty="n",col=c("black"),cex=1.3) lines(xtime, cumsum(PC2),type="l", col="blue") legend(1960,1.5, lty=1, legend=c("Cumulative PC2: 33% variance"), bty="n",col=c("blue"), cex=1.3) #Display the two ENSO modes on a world map library(maps) install.packages("mapdata") library(mapdata) plot.new() par(mfrow=c(2,1)) par(mar=c(0,0,0,0)) #Zero space between (a) and (b) map(database="world2Hires",ylim=c(-70,70), mar = c(0,0,0,0)) grid(nx=12,ny=6) points(231, -18,pch=16,cex=2, col="red") text(231, -30, "Tahiti 0.61", col="red") points(131, -12,pch=16,cex=2.6, col="blue") text(131, -24, "Darwin -0.79", col="blue") axis(2, at=seq(-70,70,20), col.axis="black", tck = -0.05, las=2, line=-0.9,lwd=0) axis(1, at=seq(0,360,60), col.axis="black",tck = -0.05, las=1, line=-0.9,lwd=0) text(180,30, "El Nino Southern Oscillation Mode 1",col="purple",cex=1.3) text(10,-60,"(a)", cex=1.4) box() par(mar=c(0,0,0,0)) #Plot mode 2 map(database="world2Hires",ylim=c(-70,70), mar = c(0,0,0,0)) grid(nx=12,ny=6) points(231, -18,pch=16,cex=2.6, col="red") text(231, -30, "Tahiti 0.79", col="red") points(131, -12,pch=16,cex=2, col="red") text(131, -24, "Darwin 0.61", col="red") text(180,30, "El Nino Southern Oscillation Mode 2",col="purple",cex=1.3) axis(2, at=seq(-70,70,20), col.axis="black", tck = -0.05, las=2, line=-0.9,lwd=0) axis(1, at=seq(0,360,60), col.axis="black",tck = -0.05, las=1, line=-0.9,lwd=0) text(10,-60,"(b)", cex=1.4) #Plot principal components of the Darwin and Tahiti SLP tdsvd<-svd(ptada) dat=provideDimnames(tdsvd$u,base=list(c("Tahiti", "Darwin"), c("EOF1","EOF2"))) lon=c(131,231) lat=c(-24,-18) ft=as.data.frame(dat) dp=ggplot(ft,aes(lon,lat)) dp1=dp+geom_point(aes(colour=factor(EOF1)),cex=9) + xlim(0,360) + ylim(-90,90) dp1#Show the plot of EOF1 dev.off() plot(seq(1951, 2015, length=780), tdsvd$v[,1], xlab="Year", ylab="WSOI1", col="red", main="PC1 as the weighted SOI") #Use NOAAGlobalTemp monthly time series as an example #Read the data online d1=read.table("https://www1.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/timeseries/aravg.ann.land_ocean.90S.90N.v4.0.1.201704.asc") dim(d1) #[1] 138 6 d1[1:3,] #V1 V2 V3 V4 V5 V6 #1 1880 -0.366861 0.009744 0.001402 0.000850 0.007492 #2 1881 -0.316453 0.009777 0.001402 0.000877 0.007498 #3 1882 -0.317072 0.009767 0.001402 0.000897 0.007468 tmean15=d1[1:136,2] #Extract the temperature data from 1880-2015 tmean15 #Show the data of the second column below #Linear trend yrtime15=seq(1880,2015) reg8015<-lm(tmean15 ~ yrtime15) # Display regression results reg8015 #Call: lm(formula = tmean15 ~ yrtime15) #Coefficients: #(Intercept) yrtime15 # -13.299762 0.006724 # Plot the temperature time series and its trend line plot(yrtime15,tmean15,xlab="Year",ylab="Temperature deg C", main="Global Annual Mean Land and Ocean Surface Temperature Anomalies 1880-2015", type="l") abline(reg8015, col="red") text(1930, 0.4, "Linear temperature trend 0.67 oC per century", col="red",cex=1.2) #Histogram and its density fit h<-hist(tmean15, main="Histogram of 1880-2015 Temperature Anomalies",xlab="Temperature anomalies") #Plot historgram xfit<-seq(min(tmean15),max(tmean15), length=30) areat=diff(h$mids[1:2])*length(tmean15) #Normalization area yfit<-areat*dnorm(xfit, mean=mean(tmean15), sd=sd(tmean15)) lines(xfit,yfit,col="blue",lwd=2) #Plot the normal fit #Kernal estimate plot(density(tmean15), main="Kernel estimate of density", xlab="Temperature") #Kernel estimate density lines(xfit,dnorm(xfit, mean=mean(tmean15), sd=sd(tmean15)), col="blue") #Moment estimated normal #Scatter plot #Change directory to where the datasets are setwd("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/chap4data-refs") ust=read.csv("USJantemp1951-2016-nohead.csv",header=FALSE) soi=read.csv("soi-data-nohead.csv", header=FALSE) #Read data soid=soi[,2] #Take the second column SOI data soim=matrix(soid,ncol=12,byrow=TRUE) #Make the SOI into a matrix with each month as a column soij=soim[,1] #Take the first column for Jan SOI ustj=ust[,3] #Take the third column: Jan US temp data plot(soij,ustj,main="Scatter plot between Jan SOI and US temp",xlab="SOI[dimensionless]", ylab="US Temp degF", pch=19) # Plot the scatter plot soiust=lm(ustj ~ soij) #Linear regression abline(soiust, col="red") #Linear trend line