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)