--- title: " " output: html_document: default pdf_document: default ---

# Theory and Applications

### SSP Shen and RCJ Somerville

Version 1.0 released in July 2019 and coded by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
Email:

Version 2.0 compiled and updated by Momtaza Sayd
San Diego State University May 2021
Version 3.0 compiled and updated by Joaquin Stawsky
San Diego State University June 2022

## Chapter 9: R Graphics for Climate Science ### Plot Fig. 9.1 ```{r} #setEPS() #postscript("fig0901.eps", height = 4, width = 8) par(mar = c(4.2,4.2,2.5,4.1)) 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", lwd = 1.5, xlab = "Year",cex.axis = .85, ylab = expression(paste(T[mean]," [", degree,"C]")), main = "Contiguous U.S. Annual Mean\nTemperature and Total Precipitation") legend(2000.5,12.42, col = c("red"),lty = 1,lwd = 2.0, legend = c(expression(paste(bold(T[mean])))),bty = "n",text.font = 2,cex = 1) #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) #Plot two legends at the same time make it difficult to adjust the font size #because of different scale #dev.off() ```

### Plot Fig. 9.2 ```{r} #Margins, math symbol, and figure setups #setEPS() #postscript("fig0902.eps", height = 4, width = 8) #Margins, math symbol, and figure setups par(mar = c(5,4.5,2.5,2.5)) 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 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 radians: ", theta-phi)),cex = 1.5, col = "red") #dev.off() ```

### Plot Fig. 9.3 ```{r} par(mar = c(8,6,3,2)) par(mgp = c(2.5,1,0)) plot(1:200/20, rnorm(200),sub = "Subtitle: 200 Random Values", xlab = "Time", ylab = "Random Values", main = "Normal Random Values", cex.lab = 1.75, cex.axis = 1.5, cex.main = 2.0, cex.sub = 1.5) #Adjust positions of axis labels par(mgp = c(2,1,0)) plot(sin,xlim = c(10,20),cex.lab = 1.2,cex.axis = 1.2) ```

### Plot Fig. 9.4 ```{r} #A fancy plot of the NOAAGlobalTemp time series setwd("∼/sshen/climmath") NOAATemp <- read.table("data/aravg.ann.land_ocean.90S.90N.v4.0.1.2016.txt", header = F) 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 [°C]", xlab = "Time",col = "red", main = "NOAA Global Average Annual Mean Temperature Anomalies", cex.lab = 1.2,cex.axis = 1.2) lines(x2,y2,type = "h", lwd = 3, tck = -0.02, col = "blue") lines(x3,y3,lwd = 2) ```

### Plot Fig. 9.5 ```{r} #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 = expression(paste("T"[mean]," [°C]"))) text(2006, 12,font = 2,"U.S. Annual Mean Temperature", cex = 1.5) text(2001,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, "U.S. Annual Total Precipitation", cex = 1.5) text(2001,840,"(b)") ```

### Figure layout for multiple panels ```{r} layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE), widths = c(3,3), heights = c(2,2)) plot(sin,type = "l", xlim = c(0,20),cex.lab = 1.2,cex.axis = 1.2) plot(sin,xlim = c(0,10),cex.lab = 1.2,cex.axis = 1.2) plot(sin,xlim = c(10,20),cex.lab = 1.2,cex.axis = 1.2) ```

### Contours and color-filled contours ```{r} x <- y <- seq(-1, 1, len = 25) z <- matrix(rnorm(25*25),nrow = 25) contour(x,y,z, main = "Contour Plot of Normal Random Values") filled.contour(x,y,z, main = "Filled Contour Plot of Normal Random Values") filled.contour(x,y,z, color.palette = heat.colors) filled.contour(x,y,z, color.palette = colorRampPalette(c("red", "white", "blue"))) ```

### Plot Fig. 9.6 ```{r,fig.width = 8, fig.height = 5.3} #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 be increasing Lon <- seq(0,355,length = 72) #Must be 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:\n5 Lat-Lon Grid", cex.lab = 1.2,cex.axis = 1.2), 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 ```

### Plot Fig. 9.7 ```{r,fig.width = 5.5,fig.height = 5} #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 over\nCanada and USA: 5° Lat-Lon Grid", xlab = "Longitude", ylab = "Latitude", cex.lab = 1.2,cex.axis = 1.2), plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()}) out.width = "50%" #specifying dimensions of output out.height = "50%" ```

### Plot Fig. 9.8 ```{r} #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 #.nc is the netCDF data format. It is a Network Common Data Form #started in 1989 at UCAR, and in the 2010s netCDF became the dominant #data format for writing big climate data. #rm(list = ls(all = TRUE)) setwd("~/sshen/climmath/") # Download netCDF file # Library #install.packages("ncdf") library(ncdf4) # 4 dimensions: lon,lat,level,time nc <- ncdf4::nc_open("data/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) # 65378 65409 65437 65468 65498 65529 library(chron) month.day.year(1297320/24,c(month = 1, day = 1, year = 1800)) #1948-01-01 precnc <- ncvar_get(nc, "air") dim(precnc) # 144 73 826... 826 months = 1948-01 to 2016-10, 68 years 10 months #plot a 90S-90N precip along a meridional line at 160E over Pacific par(mar = c(4.5,5,3,1)) plot(seq(90,-90,length = 73), precnc[65,,1], type = "l", xlab = "Latitude", ylab = "Temperature [°C]", main = "90S-90N Temperature [°C]\nAlong a Meridional Line at 160E: January 1948", lwd = 3, cex.lab = 1.2, cex.axis = 1.2) grid(nx = NULL, ny = NULL) ```

### Plot Fig. 9.9 ```{r,fig.width = 8, fig.height = 5.3} #Verify previous calculations and plots in Fig. 9.9 #Compute and plot climatology and standard deviation Jan 1948-Dec 2015 # par(mar = c(4,4,2.2,0)) 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 must arrange Lat and Lon this way. #Correspondingly, 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 the center, China on the right, and USA on the left. However, our map should #have the Pacific at the center, and USA on the 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)) #Creates same matrix as above int <- seq(-50,50,length.out = 81) rgb.palette <- colorRampPalette(c('black','blue','darkgreen','green', 'white','yellow','pink','red','maroon'),interpolate = 'spline') par(cex.axis = 1.3,cex.lab = 1.3) filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int, plot.title = title(main = "NCEP RA 1948-2015 January Climatology [°C]", xlab = "Longitude",ylab = "Latitude", cex.lab = 1.2,cex.axis = 1.2), plot.axes = {axis(1); axis(2); map('world2', add = TRUE);grid()}, key.title = title(main = "[°C]")) ``` ```{r,fig.width = 8, fig.height = 5.3} #plot standard deviation #plot.new() par(mgp = c(2,1,0)) # par(mar = c(3,3,2.2,0)) par(cex.axis = 1.3,cex.lab = 1.3) mapmat <- sdmat[,seq(length(sdmat[1,]),1)] # mapmat <- pmax(pmin(mapmat,6),0) 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 [°C]", xlab = "Longitude", ylab = "Latitude", cex.lab = 1.2,cex.axis = 1.2), plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()}, key.title = title(main = "[°C]")) ```

### Plot 9.10 ```{r,fig.width = 8, fig.height = 5.3} #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 [°C]", xlab = "Longitude",ylab = "Latitude", cex.lab = 1.2,cex.axis = 1.2), plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()}, key.title = title(main = "[°C]")) ```

### Plot Fig. 9.12 ```{r,fig.width = 7.5, fig.height = 5.3} #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') matdiff <- precnc[,,421]-climmat matdiff <- matdiff[,length(matdiff[1,]):1] matdiff <- pmax(pmin(matdiff,6),-6) filled.contour(Lon, Lat, matdiff, xlim = c(100,300), ylim = c(-50,50), color.palette = rgb.palette, levels = int, plot.title = title( main = "January 1983 Surface Air Temperature Anomaly [°C]", xlab = "Longitude",ylab = "Latitude", cex.lab = 1.2,cex.axis = 1.2), plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()}, key.title = title(main = "[°C]")) ```

### Plot Fig. 9.13 ```{r,fig.width = 8, fig.height = 5.3} #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), 12) v <- rep(c(1,-1,1,-1,1,-1), 12) par(mfrow = c(1,1)) par(mar = c(3,3,2,0.5)) 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("165°W","105°W","45°W","15°E","75°E","135°E"), 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("75°S","45°S","15°S","15°N","45°N","75°N"), 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(0, -30, "Subtropical High", col = "red") text(0, 30, "Subtropical High", col = "red") text(0, 0, "Intertropical Convergence Zone (ITCZ)", col = "red") mtext(side = 3, "Polar High", col = "red", line = 0.0, font = 1) ```

### Plot Fig. 9.14 ```{r,fig.width = 8, fig.height = 5.3} #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) #install.packages("rasterVis") #install.packages("latticeExtra") library(latticeExtra) library(rasterVis) #install.packages("raster") library(raster) library(sp) library(rgdal) #download.file("ftp://eclipse.ncdc.noaa.gov/pub/seawinds/SI/uv/clm/uvclm95to05.nc", # "uvclm95to05.nc", method = "curl") setwd("~/sshen/climmath") mincwind <- nc_open("data/uvclm95to05.nc") length(mincwind) # 14 u <- ncvar_get(mincwind, "u") v <- ncvar_get(mincwind, "v") dim(u) # 1440 719 12 #lon, lat, and month dim(v) u9 <- raster(t(u[, , 9])[ncol(u):1, ]) v9 <- raster(t(v[, , 9])[ncol(v):1, ]) filled.contour(u[, , 9],cex.lab = 1.2,cex.axis = 1.2) filled.contour(u[, , 9],cex.lab = 1.2,cex.axis = 1.2,color.palette = heat.colors) filled.contour(u[, , 9],cex.lab = 1.2,cex.axis = 1.2,color.palette = colorRampPalette(c("red", "white", "blue"))) contourplot(u[, , 9]) 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)) plot(w[],cex.lab = 1.2,cex.axis = 1.2) plot(w[],cex.lab = 1.2,cex.axis = 1.2) vectorplot(w * 10, isField = "dXY", cex.lab = 1.2,cex.axis = 1.2,region = FALSE, margin = FALSE, narrows = 10000) slope <- sqrt(w[]^2 + w[]^2) aspect <- atan2(w[], w[]) vectorplot(w*6, isField = "dXY", region = slope, cex.lab = 1.2,cex.axis = 1.2, margin = FALSE, par.settings = BuRdTheme, narrows = 10000, at = 0:10) vectorplot(stack(slope * 10, aspect), isField = TRUE, region = FALSE, margin = FALSE) ```

### Plot Fig. 9.15 ```{r,fig.width = 8, fig.height = 5.3} #ggplot for USA States library(ggplot2) states <- map_data("state") #"states" is in a data.frame p <- ggplot(data = states) + geom_polygon(aes(x = long, y = lat, fill = region, group = group), color = "white") + coord_fixed(1.3) #if fill = FALSE, the large color legend on the right is off. p <- p + xlab("Latitude")+ ylab("Longitude") p + ggtitle("Color Map of the 48 Lower States") + theme(plot.title = element_text(hjust = 0.5)) # p + labs(title = "Color Map of the 48 Lower States", # x ="Latitude", y = "Longitude") ```

### Plot Fig. 9.16 ```{r} #Plot Fig. 9.16 and animation #Free fall animation by 21 frames g <- 9.8 n <- 21 t <- seq(0,10,len = n) #install.packages("animation") library(animation) ## 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 for (i in 1:n) { plot(0, 490-(1/2)*g*(t[i])^2, pch = 19, lwd = 12, col = "black", xlab = "Horizontal location", xlim = c(-10,10), ylim = c(0,500), ylab = "Vertical Position [m]", main = paste("Free Fall Time = ", format(t[i],digits = 2, nsmall = 1), "sec") ) ani.record() # is: function (reset = FALSE, replay.cur = FALSE) } ## Now we can replay it, with an appropriate pause between frames: ## Smaller interval means faster animation. Default: interval = 1 oopts <- ani.options(interval = 0.5, ani.width = 200, ani.height = 400, title = "Free Fall" ) #Animate the frames in the plot window of R Studio ani.replay() #is: function (list) ## Show the animation on an HTML page saveHTML(ani.replay(), img.name = "Fall_animation") ```