R Code for Climate Mathematics:
Theory and Applications
A Cambridge University Press book By
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: sshen@sdsu.edu
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
#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
#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()
Plot Fig. 9.3
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
#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
#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
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
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
#Plot a 5-by-5 grid global map of standard normal random values
library(maps)## Warning: package 'maps' was built under R version 4.1.3plot.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
#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 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## File data/air.mon.mean.nc (NC_FORMAT_NETCDF4_CLASSIC):
##
## 1 variables (excluding dimension variables):
## float air[lon,lat,time] (Chunking: [144,73,1]) (Compression: shuffle,level 2)
## long_name: Monthly Mean Air Temperature at sigma level 0.995
## valid_range: -2000
## valid_range: 2000
## units: degC
## add_offset: 0
## scale_factor: 1
## missing_value: -9.96920996838687e+36
## precision: 1
## least_significant_digit: 0
## var_desc: Air Temperature
## level_desc: Surface
## statistic: Mean
## parent_stat: Individual Obs
## dataset: NCEP Reanalysis Derived Products
## actual_range: -73.7800064086914
## actual_range: 41.7490196228027
##
## 3 dimensions:
## lat Size:73
## units: degrees_north
## actual_range: 90
## actual_range: -90
## long_name: Latitude
## standard_name: latitude
## axis: Y
## lon Size:144
## units: degrees_east
## long_name: Longitude
## actual_range: 0
## actual_range: 357.5
## standard_name: longitude
## axis: X
## time Size:826 *** is unlimited ***
## long_name: Time
## delta_t: 0000-01-00 00:00:00
## avg_period: 0000-01-00 00:00:00
## prev_avg_period: 0000-00-01 00:00:00
## standard_name: time
## axis: T
## units: hours since 1800-01-01 00:00:0.0
## actual_range: 1297320
## actual_range: 1899984
##
## 8 global attributes:
## description: Data from NCEP initialized reanalysis (4x/day). These are the 0.9950 sigma level values
## platform: Model
## Conventions: COARDS
## NCO: 20121012
## history: Thu May 4 20:11:16 2000: ncrcat -d time,0,623 /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc air.mon.mean.nc
## Thu May 4 18:11:50 2000: ncrcat -d time,0,622 /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc ./surface/air.mon.mean.nc
## Mon Jul 5 23:47:18 1999: ncrcat ./air.mon.mean.nc /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc /dm/dmwork/nmc.rean.ingest/combinedMMs/surface/air.mon.mean.nc
## /home/hoop/crdc/cpreanjuke2farm/cpreanjuke2farm Mon Oct 23 21:04:20 1995 from air.sfc.gauss.85.nc
## created 95/03/13 by Hoop (netCDF2.3)
## Converted to chunked, deflated non-packed NetCDF4 2014/09
## title: monthly mean air.sig995 from the NCEP Reanalysis
## References: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.html
## dataset_title: NCEP-NCAR Reanalysis 1nc$dim$lon$vals # output values 0.0->357.5## [1] 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5
## [13] 30.0 32.5 35.0 37.5 40.0 42.5 45.0 47.5 50.0 52.5 55.0 57.5
## [25] 60.0 62.5 65.0 67.5 70.0 72.5 75.0 77.5 80.0 82.5 85.0 87.5
## [37] 90.0 92.5 95.0 97.5 100.0 102.5 105.0 107.5 110.0 112.5 115.0 117.5
## [49] 120.0 122.5 125.0 127.5 130.0 132.5 135.0 137.5 140.0 142.5 145.0 147.5
## [61] 150.0 152.5 155.0 157.5 160.0 162.5 165.0 167.5 170.0 172.5 175.0 177.5
## [73] 180.0 182.5 185.0 187.5 190.0 192.5 195.0 197.5 200.0 202.5 205.0 207.5
## [85] 210.0 212.5 215.0 217.5 220.0 222.5 225.0 227.5 230.0 232.5 235.0 237.5
## [97] 240.0 242.5 245.0 247.5 250.0 252.5 255.0 257.5 260.0 262.5 265.0 267.5
## [109] 270.0 272.5 275.0 277.5 280.0 282.5 285.0 287.5 290.0 292.5 295.0 297.5
## [121] 300.0 302.5 305.0 307.5 310.0 312.5 315.0 317.5 320.0 322.5 325.0 327.5
## [133] 330.0 332.5 335.0 337.5 340.0 342.5 345.0 347.5 350.0 352.5 355.0 357.5nc$dim$lat$vals #output values 90->-90## [1] 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 70.0 67.5 65.0 62.5
## [13] 60.0 57.5 55.0 52.5 50.0 47.5 45.0 42.5 40.0 37.5 35.0 32.5
## [25] 30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5
## [37] 0.0 -2.5 -5.0 -7.5 -10.0 -12.5 -15.0 -17.5 -20.0 -22.5 -25.0 -27.5
## [49] -30.0 -32.5 -35.0 -37.5 -40.0 -42.5 -45.0 -47.5 -50.0 -52.5 -55.0 -57.5
## [61] -60.0 -62.5 -65.0 -67.5 -70.0 -72.5 -75.0 -77.5 -80.0 -82.5 -85.0 -87.5
## [73] -90.0# nc$dim$time$vals
# nc$dim$time$units
nc$dim$level$vals## NULLLon <- ncvar_get(nc, "lon")
Lat1 <- ncvar_get(nc, "lat")
Time <- ncvar_get(nc, "time")
head(Time)## [1] 1297320 1298064 1298760 1299504 1300224 1300968#[1] 65378 65409 65437 65468 65498 65529
library(chron)## Warning: package 'chron' was built under R version 4.1.3month.day.year(1297320/24,c(month = 1, day = 1, year = 1800))## $month
## [1] 1
##
## $day
## [1] 1
##
## $year
## [1] 1948#1948-01-01
precnc <- ncvar_get(nc, "air")
dim(precnc)## [1] 144 73 826#[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)
Plot Fig. 9.9
#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]"))
#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
#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
#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
#Wind directions due to the balance between PGF and Coriolis force
#using an arrow plot for vector fields on a map
library(fields)## Warning: package 'fields' was built under R version 4.1.3## Loading required package: spam## Warning: package 'spam' was built under R version 4.1.3## Spam version 2.8-0 (2022-01-05) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: viridis## Warning: package 'viridis' was built under R version 4.1.3## Loading required package: viridisLite## Warning: package 'viridisLite' was built under R version 4.1.3##
## Attaching package: 'viridis'## The following object is masked from 'package:maps':
##
## unemp##
## Try help(fields) to get started.library(maps)
library(mapproj)## Warning: package 'mapproj' was built under R version 4.1.3lat <- 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
#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)## Warning: package 'latticeExtra' was built under R version 4.1.3library(rasterVis)## Warning: package 'rasterVis' was built under R version 4.1.3#install.packages("raster")
library(raster)## Warning: package 'raster' was built under R version 4.1.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.1.3library(sp)
library(rgdal)## Warning: package 'rgdal' was built under R version 4.1.3## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-29, (SVN revision 1165M)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/HP/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/HP/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/HP/Documents/R/win-library/4.1/rgdal/proj#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] 15#[1] 14
u <- ncvar_get(mincwind, "u")
v <- ncvar_get(mincwind, "v")
dim(u)## [1] 1440 719 12#[1] 1440 719 12 #lon, lat, and month
dim(v)## [1] 1440 719 12u9 <- 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)## [1] 0.00 359.75range(wlat)## [1] -89.75 89.75projection(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)
Plot Fig. 9.15
#ggplot for USA States
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.1.3##
## Attaching package: 'ggplot2'## The following object is masked from 'package:latticeExtra':
##
## layerstates <- 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
#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)## Warning: package 'animation' was built under R version 4.1.3## 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")## HTML file created at: index.html