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 11: Climate Data Matrices and Linear Algebra
Read NOAAGlobalTemp and form the space-time data matrix
#NOAAGlobalTemp dataset since 1880: A merged land SAT and ocean SST anomalies
#The monthly 5-deg gridded data and their global average can be
#downloaded from the NOAA website.
#The gridded data has 43MB
#The data can also be downloaded from the website of this book
#https://climatemathematics.sdsu.edu
#The file name is NOAAGlobalTemp.gridded.v4.0.1.201701.asc
rm(list = ls(all = TRUE))
# Download .asc file
#.asc is an ASCII data format
#This book uses "scan" to read the asc data
#and then write the data into a space-time matrix
#Climate Mathematics treats space-time matrix as
#its data standard while big climate data use .nc
setwd("~/sshen/climmath/data")
da1 <- scan("NOAAGlobalTemp.gridded.v4.0.1.201701.asc")
#From Jan 1880 to Jan 2017
length(da1)## [1] 4267130#[1] 4267130
da1[1:3]## [1] 1.0 1880.0 -999.9#[1] 1.0 1880.0 -999.9 #means month, 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)## [1] 1645length(tm2)## [1] 1645mm1 <- da1[tm1] #Extract months
yy1 <- da1[tm2] #Extract years
head(mm1)## [1] 1 2 3 4 5 6head(yy1)## [1] 1880 1880 1880 1880 1880 1880length(mm1)## [1] 1645length(yy1)## [1] 1645rw1 <- paste(yy1, sep = "-", mm1) #Combine YYYY with MM
head(tm1)## [1] 1 2595 5189 7783 10377 12971head(tm2)## [1] 2 2596 5190 7784 10378 12972tm3 <- cbind(tm1,tm2)
tm4 <- as.vector(t(tm3))
head(tm4)## [1] 1 2 2595 2596 5189 5190#[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#[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)
dim(da3)## [1] 2592 1645#[1] 2592 1645
#Put space-time coordinates in the space-time data da3
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#[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.
gpcpst[1:3,1:6] #Part of the data## LAT LON 1880-1 1880-2 1880-3 1880-4
## [1,] -87.5 2.5 -999.9 -999.9 -999.9 -999.9
## [2,] -87.5 7.5 -999.9 -999.9 -999.9 -999.9
## [3,] -87.5 12.5 -999.9 -999.9 -999.9 -999.9# LAT LON 1880-1 1880-2 1880-3 1880-4
#[1,] -87.5 2.5 -999.9 -999.9 -999.9 -999.9
#[2,] -87.5 7.5 -999.9 -999.9 -999.9 -999.9
#[3,] -87.5 12.5 -999.9 -999.9 -999.9 -999.9
write.csv(gpcpst,file = "NOAAGlobalT.csv")
#Output the data as a csv file
Plot Fig. 11.4: Dec 2015 global surface temperature anomalies map
library(maps) # install maps package first if not done before## Warning: package 'maps' was built under R version 4.1.3Lat <- 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
#Convert the vector into a lon-lat matrix for R map plotting
mapmat <- pmax(pmin(mapmat,6),-6) #Put values between -6 and 6
#Matrix flipping is not needed since the data go from 2.5 to 375.5
#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')
filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
plot.title = title(main = "NOAAGlobalTemp Anomalies Dec. 2015 [°C]",
xlab = "Longitude",ylab = "Latitude", cex.lab = 1.5),
plot.axes = {axis(1, cex.axis = 1.25);
axis(2, cex.axis = 1.25);map('world2', add = TRUE);grid()},
key.title = title(main = "[°C]"),
key.axes = {axis(4, cex.axis = 1.25)})
Plot Fig. 11.5: Dec 1997 tropical Pacific SST anomalies map
#Select only the data for the tropical Pacific region
n2 <- which(gpcpst[,1]>-20&gpcpst[,1]<20&gpcpst[,2]>160&gpcpst[,2]<260)
dim(gpcpst)## [1] 2592 1647length(n2)## [1] 160#[1] 160 ( = 8 latitude bends and 20 longitude bends)
pacificdat <- gpcpst[n2,855:1454]
#plot.new()
Lat <- seq(-17.5,17.5, by = 5)
Lon <- seq(162.5, 257.5, by = 5)
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 [°C]: Dec. 1997",
xlab = "Longitude",ylab = "Latitude", cex.lab = 1.25),
plot.axes = {axis(1, cex.axis = 1.25); axis(2, cex.axis = 1.25);
map('world2', add = TRUE);grid()},
key.title = title(main = "[°C]"),
key.axes = {axis(4, cex.axis = 1.25)})
#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 = "Temperature [°C]",
main = "San Diego Temperature Anomalies History")
lm(SanDiegoData ~ seq(1880,2017, len = length(SanDiegoData)))##
## Call:
## lm(formula = SanDiegoData ~ seq(1880, 2017, len = length(SanDiegoData)))
##
## Coefficients:
## (Intercept)
## -15.060295
## seq(1880, 2017, len = length(SanDiegoData))
## 0.007624n2 <- which(gpcpst[,1] == 52.5&gpcpst[,2] == 247.5)
EdmontonData <- gpcpst[n2,855:1454]
plot(seq(1880,2017, len = length(EdmontonData)),
EdmontonData, type = "l",
xlab = "Year", ylab = "Temperature [°C]",
main = "Temperature Anomalies History of Edmonton, Canada")
lm(EdmontonData ~ seq(1880,2017, len = length(EdmontonData)))##
## Call:
## lm(formula = EdmontonData ~ seq(1880, 2017, len = length(EdmontonData)))
##
## Coefficients:
## (Intercept)
## -23.26528
## seq(1880, 2017, len = length(EdmontonData))
## 0.01178#Compute the area-weighted average of the gridded data
#36-by-72 boxes and Jan1880-Jan2016 = 1633 months + lat and lon
#Compute the area-weight for each box and each month
#that has data. Thus the area-weight is a matrix.
areaw <- matrix(0,nrow = 2592,ncol = 1647)
dim(areaw)## [1] 2592 1647#[1] 2592 1647
temp <- gpcpst
areaw[,1] <- temp[,1]
areaw[,2] <- temp[,2]
veca <- cos(areaw[,1]*pi/180) #convert deg into radian
#create an area-weight matrix equal to cosine for the box with data
#and zero for the box with missing data
for(j in 3:1647) {for (i in 1:2592) {if(temp[i,j]> -290.0) {areaw[i,j] = veca[i]} }}
#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 Fig. 11.6: Global average temperature and its trend
plot.new()
timemo <- seq(1880,2017,length = 1645)
par(mar = c(3.5,3.5,3,0.5))
par(mgp = c(2.3,1.0,0.0))
plot(timemo,avev,type = "l", cex.lab = 1.25,
xlab = "Year", ylab = "Temperature Anomaly [°C]",
main = "Area-weighted Global Average of Monthly\nSAT Anomalies: Jan. 1880-Jan. 2017")
abline(lm(avev ~ timemo),col = "blue",lwd = 2)
text(1930,0.7, "Linear trend: 0.69 [°C]/Century",
cex = 1.25, col = "blue")
Plot Fig. 11.7: Compare the warming trends of Edmonton and San Diego
#Extract data for a specified box with given lat and lon
n2 <- which(gpcpst[,1] == 52.5&gpcpst[,2] == 247.5)
dedm <- gpcpst[n2,855:1454]
t <- seq(1880,2017, len = length(dedm))
par(mar = c(4,4.5,2.5,1))
plot(t, dedm, type = "l", ylim = c(-15,15),col = "red",
cex.lab = 1.15, cex.axis = 1.15,
xlab = "Year", ylab = "Temperature Anomalies [°C]")
title("Monthly Temperature Anomalies History of\nEdmonton, Canada, and San Diego, USA",
line = 0.5,cex.main = 1.25)
regedm <- lm(dedm ~ t)
abline(regedm,col = "red", lwd = 3)
#San Diego data
n3 <- which(gpcpst[,1] == 32.5&gpcpst[,2] == 242.5)
dsan <- gpcpst[n3,855:1454]
lines(t, dsan, type = "l", col = "blue")
regsan <- lm(dsan ~ t)
abline(regsan, col = "blue", lwd = 3)
legend(1880,16, col = c("red","blue"),lty = 1,lwd = 2.0,
legend = c("Edmonton, Canada: Trend 1.18 °C/Century, SD 3.01 °C",
"San Diego, USA: Trend 0.76 °C/Century, SD 0.87 °C"),
bty = "n",text.font = 2,cex = 1.15)
regedm##
## Call:
## lm(formula = dedm ~ t)
##
## Coefficients:
## (Intercept) t
## -23.26528 0.01178#0.01178 #Edmonton trend
regsan##
## Call:
## lm(formula = dsan ~ t)
##
## Coefficients:
## (Intercept) t
## -15.060295 0.007624#0.007624 #San Diego trend
sd(dedm)## [1] 3.006863#[1] 3.006863
sd(dsan)## [1] 0.8653399#[1] 0.8653399
Plot Fig. 11.8: Display the trends of 12 months using 12 panels on the same figure
#Read directly from the NCEI URL for the NOAAGlobalTemp monthly global average data:
#You need to find the link to the data file.
#Data file: aravg.mon.land_ocean.90S.90N.v4.0.1.201703.txt
setwd("~/sshen/climmath/data")
aveNCEI <- read.table("aravg.mon.land_ocean.90S.90N.v4.0.1.201703.txt",
header = FALSE)
dim(aveNCEI) #Jan 1880-Feb 2017 #an extra month to be deleted## [1] 1647 10#[1] 1647 10
par(mar = c(4,5,2,1))
timemo <- seq(aveNCEI[1,1],aveNCEI[length(aveNCEI[,1]),1],len = length(aveNCEI[,1]) )
#create matrix of 136 yrs of data matrix
# row = year from 1880 to 2016, col = mon 1 to 12
ave <- aveNCEI[,3]
myear <- length(ave)/12
nyear <- floor(myear)
nmon <- nyear*12
avem <- matrix(ave[1:nmon], ncol = 12, byrow = TRUE)
#compute annual average
annv <- seq(0,length = nyear)
for(y in 1:nyear){annv[y] <- mean(avem[y,])}
#Put the monthly averages and annual ave in a matrix
avemy <- cbind(avem,annv)
#Plot 12 panels on the same figure: Trend for each month
dev.off()## null device
## 1#quartz(width = 10,height = 16, pointsize = 16)
#quartz(display = "name", width = 5, height = 5, pointsize = 12,
# family = "Helvetica", antialias = TRUE, autorefresh = TRUE)
plot.new()
#png(file = 'monthtrend.png') #Automatically saving figure
timeyr <- seq(aveNCEI[1,1], aveNCEI[1,1]+nyear-1)
par(mfrow = c(4, 3)) # 4 rows and 3 columns
par(mgp = c(2,1,0))
for (i in 1:12) {
plot(timeyr, avemy[,i],type = "l", ylim = c(-1.0,1.0),
xlab = "Year",ylab = "Temperature [°C]", cex.lab = 1.15,
main = paste("Month =", i, split = ""))
abline(lm(avemy[,i]~timeyr),col = "red")
text(1945,0.7, paste("Trend = ",
round(digits = 3,
(100*coefficients(lm(avemy[,i]~timeyr))[2]))," °C/Century" ),
col = "red", cex = 1.25)
}
#dev.off()
Plot Fig. 11.9: Global average annual mean SAT anomalies
par(mar = c(4,4.5,3,1.0))
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.,cex.axis = 1, lwd = 2,
xlab = "Year", ylab = "Temperature Anomaly [°C]",
)
title("Area-weighted Global Average of Annual\nSAT Anomalies: 1880-2016", line = 0.5,cex.main = 1.15)
abline(lm(annv ~ timeyr),col = "blue",lwd = 2)
text(1923,0.4, "Linear trend: 0.69 [°C]/Century",
cex = 1, col = "blue")
text(1900,0.07, "Base line",cex = 1, col = "red")
lines(timeyr,rep(0,137), type = "l",col = "red")
Plot Fig. 11.10: Fit polynomials to the global average annual mean
#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 orthogonal polynomial of 9th order
#Orthogonal polynomial fitting is usually better
par(mar = c(4,4,3,2))
plot(timeyr,annv,type = "s",
cex.lab = 1.05,cex.axis = 1.05,cex.main = 1.25, lwd = 2,
xlab = "Year", ylab = "Temperature Anomaly [°C]",
main = "Annual SAT Time Series and its Orthogonal\nPolynomial Fits: 1880-2016")
lines(timeyr,predict(polyor9),col = "blue", lwd = 3)
legend(1880, 0.6, col = c("blue"),lty = 1,lwd = 1.0,
legend = c("9th Order Orthogonal Polynomial Fit"),
bty = "n",text.font = 2,cex = 1.05)
lines(timeyr,predict(polyor20),col = "red", lwd = 3)
legend(1880, 0.7, col = c("red"),lty = 1,lwd = 1.0,
legend = c("20th Order Orthogonal Polynomial Fit"),
bty = "n",text.font = 2,cex = 1.05)
#Deal with missing data NA
x <- 1:8
y <- c(2,4,NA,3,6.8,NA,NA,9)
fitted(lm(y ~ x, na.action = na.exclude))## 1 2 3 4 5 6 7 8
## 2.08 3.04 NA 4.96 5.92 NA NA 8.80# 1 2 3 4 5 6 7 8
#2.08 3.04 NA 4.96 5.92 NA NA 8.80
##
fitted(lm(y ~ x, na.action = na.omit))## 1 2 4 5 8
## 2.08 3.04 4.96 5.92 8.80# 1 2 4 5 8
#2.08 3.04 4.96 5.92 8.80
Plot Fig. 11.11: Spatial pattern of the linear temporal trend
#No missing data for the first month (January 1900) and
#the last month (December 1999)
#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,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: [°C/Century]",
xlab = "Longitude",ylab = "Latitude", cex.lab = 1.25),
plot.axes = {axis(1, cex.axis = 1.25); axis(2, cex.axis = 1.25);
map('world2', add = TRUE);grid()},
key.title = title(main = ""),
key.axes = {axis(4, cex.axis = 1.25)})
#Example of computing the trend for a grid bbox
i <- 600
timemo1 <- seq(1900, 1999, length = 1200)
lm(temp1[i,243:1442] ~ timemo1, na.action = na.omit)##
## Call:
## lm(formula = temp1[i, 243:1442] ~ timemo1, na.action = na.omit)
##
## Coefficients:
## (Intercept) timemo1
## -18.795817 0.009423lm(temp1[i,243:1442] ~ timemo1, na.action = na.exclude)##
## Call:
## lm(formula = temp1[i, 243:1442] ~ timemo1, na.action = na.exclude)
##
## Coefficients:
## (Intercept) timemo1
## -18.795817 0.009423#temp1[i,243:1442]
# 1900-1 1900-2 1900-3 1900-4 1900-5
#-0.7457 NA -1.4406 -1.0936 -0.8193
Plot Fig. 11.12: Spatial pattern of the linear temporal trend
#under a relaxed condition: allowing up to 1/3 data missing
#Trend for each box for the 20th century: Version 2: Allow 2/3 of data
#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: [°C/Decade]",
xlab = "Longitude",ylab = "Latitude", cex.lab = 1.25),
plot.axes = {axis(1, cex.axis = 1.25); axis(2, cex.axis = 1.25);
map('world2', add = TRUE);grid()},
key.title = title(main = ""),
key.axes = {axis(4, cex.axis = 1.25)})
Plot Fig. 11.13: Map of the 1976-2016 trend of monhtly data
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 the 1976-2016 trend map
plot.new()
#par(mar = c(4.5,4.5,3,0))
mapmat <- matrix(120*trend7616,nrow = 72)
mapmat <- pmax(pmin(mapmat,4),-4)
int <- seq(-4,4,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. 1976-Dec. 2016 Temperature Trends: [°C/Decade]",
xlab = "Longitude",ylab = "Latitude", cex.lab = 1.25),
plot.axes = {axis(1, cex.axis = 1.25); axis(2, cex.axis = 1.25);
map('world2', add = TRUE);grid()}, key.title = title(main = ""),
key.axes = {axis(4, cex.axis = 1.25)})
Appendix B: Cross Product
#GPS-Planimeter example: Compute the area of an ellipse based on
#the coordinate data on the boundary using Green's Theorem
n <- 1000
a <- 4
b <- 2
t <- seq(0,2*pi,length = n+1)
x <- a*cos(t)
y <- b*sin(t)
s <- rep(0,n)
for (i in 1:n){s[i] <- -y[i]*x[i+1] + x[i]*y[i+1]}
A <- 0.5*sum(s)
A## [1] 25.13258#[1] 25.13258
#The aera of an ellipse according to the formula: pi ab
pi*a*b## [1] 25.13274#[1] 25.13274