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
https://shen.sdsu.edu
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 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] 1645
length(tm2)
## [1] 1645
mm1 <- da1[tm1] #Extract months
yy1 <- da1[tm2] #Extract years
head(mm1)
## [1] 1 2 3 4 5 6
head(yy1)
## [1] 1880 1880 1880 1880 1880 1880
length(mm1)
## [1] 1645
length(yy1)
## [1] 1645
rw1 <- paste(yy1, sep = "-", mm1) #Combine YYYY with MM
head(tm1)
## [1]     1  2595  5189  7783 10377 12971
head(tm2)
## [1]     2  2596  5190  7784 10378 12972
tm3 <- 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.3
Lat <- 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 1647
length(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.007624
n2 <- 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.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.009423
lm(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