#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
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)})
#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.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")
#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)