---
title: "&nbsp;"
output:
  html_document: default
  pdf_document: default
---

<center><h1>R Code for Climate Mathematics:</h1></center>
<center><h1>Theory and Applications</h1></center>
<p>&nbsp;</p>
<center><h3>A Cambridge University Press book By</h3></center>
<center><h3>SSP Shen and RCJ Somerville</h3></center>

<p>&nbsp;</p>
<p>&nbsp;</p>

<div align="center">Version 1.0 released in July 2019 and coded by Dr. Samuel Shen, Distinguished Professor</div>
<div align="center">San Diego State University, California, USA</div> 
<div align="center"><https://shen.sdsu.edu></div> 
<div align="center">Email: <sshen@sdsu.edu></div> 
<p>&nbsp;</p>
<div align="center">Version 2.0 compiled and updated by Momtaza Sayd</div> 
<div align="center">San Diego State University May 2021</div>
<div align="center">Version 3.0 compiled and updated by Joaquin Stawsky</div>
<div align="center">San Diego State University June 2022</div>

<p>&nbsp;</p>
<p>&nbsp;</p>
<p>&nbsp;</p>
## 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()
```

<p>&nbsp;</p>
### 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[0])),cex = 1.5, col = "red")
#dev.off()
```

<p>&nbsp;</p>
### 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)
```

<p>&nbsp;</p>
### 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)
```

<p>&nbsp;</p>
### 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)")
```

<p>&nbsp;</p>
### 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)
```

<p>&nbsp;</p>
### 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")))
```

<p>&nbsp;</p>
### 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 
```

<p>&nbsp;</p>
### 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%"
```

<p>&nbsp;</p>
### 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)
#[1] 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)
#[1] 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)
```

<p>&nbsp;</p>
### 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]"))
```


<p>&nbsp;</p>
### 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]"))
```

<p>&nbsp;</p>
### 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]"))
```

<p>&nbsp;</p>
### 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)
```

<p>&nbsp;</p>
### 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)
#[1] 14
u <- ncvar_get(mincwind, "u")
v <- ncvar_get(mincwind, "v")
dim(u)
#[1] 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[[1]],cex.lab = 1.2,cex.axis = 1.2)
plot(w[[2]],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[[1]]^2 + w[[2]]^2)
aspect <- atan2(w[[1]], w[[2]])
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)
```

<p>&nbsp;</p>
### 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")

```

<p>&nbsp;</p>
### 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")
```

