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
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("/Users/sshen/climmath")
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 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)## [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:
#https://www1.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/timeseries
#Data file: aravg.mon.land_ocean.90S.90N.v4.0.1.201703.txt
#setwd("/Users/sshen/Desktop/MyDocs/teach/SIOC290-ClimateMath2016/Rcodes/NOAAGlobalTemp")
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()## null device
## 1
Plot Fig. 11.9: Global average annual mean
par(mar = c(4,4.5,2,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.125,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 orthongonal 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