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: 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 10: Advanced R Analysis and Plotting: EOFs, Trends, and Global Data

Fig. 10.1

x <- seq(0, 2*pi, len = 100)
y <- seq(0, 2*pi, len = 100)
mydat <- array(0,dim = c(100,100,10))
for(t in 1:10){
  z <- function(x,y){z <- sin(t)*(1/pi)*sin(x)*sin(y)+exp(-0.3*t)*(1/pi)*sin(8*x)*sin(8*y)}
  mydat[,,t] <- outer(x,y,z)
}

#Plot the original z(x,y,t) waves for a given t
filled.contour(x,y,mydat[,,10], color.palette = rainbow, 
               plot.title = title(main = "Original field at t = 10", 
                                xlab = "x", ylab = "y", cex.lab = 1.2),
               key.title = title(main = "Scale"),
               plot.axes =  {axis(1,seq(0,3*pi, by = 1), cex = 1.2) 
                 axis(2,seq(0, 2*pi, by = 1), cex = 1.2)}
)

#Space-time data decomposition by SVD
da1 <- matrix(0,nrow = length(x)*length(y),ncol = 10)
for (i in 1:10) {da1[,i] = c(t(mydat[,,i]))}
da2 <- svd(da1)
uda2 <- da2$u
vda2 <- da2$v
dda2 <- da2$d
dda2

##  [1] 3.589047e+01 1.596154e+01 1.685814e-13 7.407470e-14 1.720475e-15
##  [6] 1.216635e-15 1.036319e-15 9.101978e-16 6.840079e-16 6.105709e-16

 

Plot Fig. 10.2

#Plot EOF modes
#Plot EOF1
par(mar = c(4,4,2.5,4))
par(mgp = c(2,1,0))
int <- seq(-0.3,0.3,length.out = 21)
rgb.palette <- colorRampPalette(c('black','purple','blue','white', 
                               'green', 'yellow','pink','red','maroon'),
                             interpolate = 'spline')
filled.contour(x,y,matrix(-uda2[,1],nrow = 100), color.palette = rgb.palette,
               plot.title = title(main = "SVD Mode 1: EOF1", 
                                xlab = "x", ylab = "y", cex.lab = 1.2),
               key.title = title(main = "Scale"),
               plot.axes = {axis(1,seq(0,2*pi, by = 1), cex = 1.2) 
                 axis(2,seq(0, 2*pi, by = 1), cex = 1.2)})

#Plot EOF2
par(mar = c(4,4,2.5,4))
par(mgp = c(2,1,0))
int <- seq(-0.3,0.3,length.out = 21)
rgb.palette <- colorRampPalette(c('black','purple','blue','white', 
                               'green', 'yellow','pink','red','maroon'),
                             interpolate = 'spline')
filled.contour(x,y,matrix(-uda2[,2],nrow = 100), color.palette = rgb.palette,
               plot.title = title(main = "SVD Mode 2: EOF2", 
                                xlab = "x", ylab = "y", cex.lab = 1.2),
               key.title = title(main = "Scale"),
               plot.axes =  {axis(1,seq(0,2*pi, by = 1), cex = 1.2) 
                 axis(2,seq(0, 2*pi, by = 1), cex = 1.2)})

#Plot the theoretical orthogonal modes
z1 <- function(x,y){(1/pi)*sin(x)*sin(y)}
z2 <- function(x,y){(1/pi)*sin(8*x)*sin(8*y)}
fcn1 <- outer(x,y,z1)
fcn2 <- outer(x,y,z2)
#Plot accurate Mode 1
par(mar = c(4,4,2.5,4))
par(mgp = c(2,1,0))
int <- seq(-0.3,0.3,length.out = 21)
rgb.palette <- colorRampPalette(c('black','purple','blue','white', 
                               'green', 'yellow','pink','red','maroon'),
                             interpolate = 'spline')
#color.palette  = rainbow
filled.contour(x,y,fcn1, color.palette = rgb.palette,
               plot.title = title(main = "Accurate Mode 1", 
                                xlab = "x", ylab = "y", cex.lab = 1.2),
               key.title = title(main = "Scale"),
               plot.axes =  {axis(1,seq(0,3*pi, by = 1), cex = 1.2) 
                 axis(2,seq(0, 2*pi, by = 1), cex = 1.2)}
)

 

Plot accurate Mode 2

par(mar = c(4,4,2.5,4))
par(mgp = c(2,1,0))
int = seq(-0.3,0.3,length.out = 21)
rgb.palette = colorRampPalette(c('black','purple','blue','white', 
                               'green', 'yellow','pink','red','maroon'),
                             interpolate = 'spline')
filled.contour(x,y,fcn2, color.palette = rgb.palette, 
               plot.title = title(main = "Accurate Mode 2", 
                                xlab = "x", ylab = "y", cex.lab = 1.2),
               key.title = title(main = "Scale"),
               plot.axes =  {axis(1,seq(0,3*pi, by = 1), cex = 1.2) 
                 axis(2,seq(0, 2*pi, by = 1), cex = 1.2)}
)

 

Plot Fig. 10.3

#Plot PCs and coefficients of the functional patterns
t <- 1:10
plot(t, vda2[,1],type = "o", ylim = c(-1,1), lwd = 2, 
     ylab = "PC or Coefficient", xlab = "Time",
     main = "SVD PCs vs. Accurate Temporal Coefficients")
legend(1,.7, lty = 1, legend = c("PC1: 83% variance"), 
       bty = "n",col = c("black"))
lines(t, vda2[,2],type = "o", col = "red", lwd = 2)
legend(1,.55, lty = 1, legend = c("PC2: 17% variance"), 
       col = "red", bty = "n",text.col = c("red"))
lines(t, -sin(t), col = "blue", type = "o")
legend(1,1, lty = 1, legend = c("Mode 1 Coefficient: 91% variance"), 
       col = "blue", bty = "n",text.col = "blue")
lines(t, -exp(-0.3*t), type = "o",col = "purple")
legend(1,.85, lty = 1, legend = c("Mode 2 Coefficient: 9% variance"), 
       col = "purple", bty = "n",text.col = "purple")

 

Verify orthogonality of PCs

t(vda2[,1])%*%vda2[,2]

##               [,1]
## [1,] -1.665335e-16

#[1,] -1.665335e-16
t <- 1:10
t(-sin(t))%*%(-exp(-0.3*t))

##           [,1]
## [1,] 0.8625048

#[1,] 0.8625048

#Data reconstruction by two modes
B <- uda2[,1:2]%*%diag(dda2)[1:2,1:2]%*%t(vda2[,1:2])

#Data reconstruction by all the ten modes
B1 <- uda2%*%diag(dda2)%*%t(vda2)

max(B-B1)

## [1] 1.10875e-13

#[1] 1.10875e-13 implies B = B1

 

Plot Fig. 10.4: The two-mode reconstruction by SVD

#plot.new()
filled.contour(x,y,matrix(B[,5],nrow = 100), color.palette = rainbow, 
               plot.title = title(main = "2-mode SVD Reconstructed Field t = 5",
                                xlab = "x", ylab = "y", cex.lab = 1.2),
               key.title = title(main = "Scale"),
               plot.axes = {axis(1,seq(0,3*pi, by = 1), cex = 1.2) 
                 axis(2,seq(0, 2*pi, by = 1), cex = 1.2)})

# Read netCDF data files
#install.packages("ncdf")
library(ncdf4)
#Download netCDF dataset air.mon.mean.nc
#https://climatemathematics.sdsu.edu/data/air.mon.mean.nc
#Move air.mon.mean.nc to directory /Users/sshen/data
# 4 dimensions: lon,lat,level,time
setwd("~/sshen/climmath")
nc = nc_open("data/air.mon.mean.nc")
nc

## File data/air.mon.mean.nc (NC_FORMAT_NETCDF4_CLASSIC):
## 
##      1 variables (excluding dimension variables):
##         float air[lon,lat,time]   (Chunking: [144,73,1])  (Compression: shuffle,level 2)
##             long_name: Monthly Mean Air Temperature at sigma level 0.995
##             valid_range: -2000
##              valid_range: 2000
##             units: degC
##             add_offset: 0
##             scale_factor: 1
##             missing_value: -9.96920996838687e+36
##             precision: 1
##             least_significant_digit: 0
##             var_desc: Air Temperature
##             level_desc: Surface
##             statistic: Mean
##             parent_stat: Individual Obs
##             dataset: NCEP Reanalysis Derived Products
##             actual_range: -73.7800064086914
##              actual_range: 41.7490196228027
## 
##      3 dimensions:
##         lat  Size:73 
##             units: degrees_north
##             actual_range: 90
##              actual_range: -90
##             long_name: Latitude
##             standard_name: latitude
##             axis: Y
##         lon  Size:144 
##             units: degrees_east
##             long_name: Longitude
##             actual_range: 0
##              actual_range: 357.5
##             standard_name: longitude
##             axis: X
##         time  Size:826   *** is unlimited *** 
##             long_name: Time
##             delta_t: 0000-01-00 00:00:00
##             avg_period: 0000-01-00 00:00:00
##             prev_avg_period: 0000-00-01 00:00:00
##             standard_name: time
##             axis: T
##             units: hours since 1800-01-01 00:00:0.0
##             actual_range: 1297320
##              actual_range: 1899984
## 
##     8 global attributes:
##         description: Data from NCEP initialized reanalysis (4x/day).  These are the 0.9950 sigma level values
##         platform: Model
##         Conventions: COARDS
##         NCO: 20121012
##         history: Thu May  4 20:11:16 2000: ncrcat -d time,0,623 /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc air.mon.mean.nc
## Thu May  4 18:11:50 2000: ncrcat -d time,0,622 /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc ./surface/air.mon.mean.nc
## Mon Jul  5 23:47:18 1999: ncrcat ./air.mon.mean.nc /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc /dm/dmwork/nmc.rean.ingest/combinedMMs/surface/air.mon.mean.nc
## /home/hoop/crdc/cpreanjuke2farm/cpreanjuke2farm Mon Oct 23 21:04:20 1995 from air.sfc.gauss.85.nc
## created 95/03/13 by Hoop (netCDF2.3)
## Converted to chunked, deflated non-packed NetCDF4 2014/09
##         title: monthly mean air.sig995 from the NCEP Reanalysis
##         References: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.html
##         dataset_title: NCEP-NCAR Reanalysis 1

nc$dim$lon$vals #output lon values 0.0->357.5

##   [1]   0.0   2.5   5.0   7.5  10.0  12.5  15.0  17.5  20.0  22.5  25.0  27.5
##  [13]  30.0  32.5  35.0  37.5  40.0  42.5  45.0  47.5  50.0  52.5  55.0  57.5
##  [25]  60.0  62.5  65.0  67.5  70.0  72.5  75.0  77.5  80.0  82.5  85.0  87.5
##  [37]  90.0  92.5  95.0  97.5 100.0 102.5 105.0 107.5 110.0 112.5 115.0 117.5
##  [49] 120.0 122.5 125.0 127.5 130.0 132.5 135.0 137.5 140.0 142.5 145.0 147.5
##  [61] 150.0 152.5 155.0 157.5 160.0 162.5 165.0 167.5 170.0 172.5 175.0 177.5
##  [73] 180.0 182.5 185.0 187.5 190.0 192.5 195.0 197.5 200.0 202.5 205.0 207.5
##  [85] 210.0 212.5 215.0 217.5 220.0 222.5 225.0 227.5 230.0 232.5 235.0 237.5
##  [97] 240.0 242.5 245.0 247.5 250.0 252.5 255.0 257.5 260.0 262.5 265.0 267.5
## [109] 270.0 272.5 275.0 277.5 280.0 282.5 285.0 287.5 290.0 292.5 295.0 297.5
## [121] 300.0 302.5 305.0 307.5 310.0 312.5 315.0 317.5 320.0 322.5 325.0 327.5
## [133] 330.0 332.5 335.0 337.5 340.0 342.5 345.0 347.5 350.0 352.5 355.0 357.5

nc$dim$lat$vals #output lat values 90->-90

##  [1]  90.0  87.5  85.0  82.5  80.0  77.5  75.0  72.5  70.0  67.5  65.0  62.5
## [13]  60.0  57.5  55.0  52.5  50.0  47.5  45.0  42.5  40.0  37.5  35.0  32.5
## [25]  30.0  27.5  25.0  22.5  20.0  17.5  15.0  12.5  10.0   7.5   5.0   2.5
## [37]   0.0  -2.5  -5.0  -7.5 -10.0 -12.5 -15.0 -17.5 -20.0 -22.5 -25.0 -27.5
## [49] -30.0 -32.5 -35.0 -37.5 -40.0 -42.5 -45.0 -47.5 -50.0 -52.5 -55.0 -57.5
## [61] -60.0 -62.5 -65.0 -67.5 -70.0 -72.5 -75.0 -77.5 -80.0 -82.5 -85.0 -87.5
## [73] -90.0

nc$dim$time$vals #output time values in GMT hours: 1297320, 1298064, ...

##   [1] 1297320 1298064 1298760 1299504 1300224 1300968 1301688 1302432 1303176
##  [10] 1303896 1304640 1305360 1306104 1306848 1307520 1308264 1308984 1309728
##  [19] 1310448 1311192 1311936 1312656 1313400 1314120 1314864 1315608 1316280
##  [28] 1317024 1317744 1318488 1319208 1319952 1320696 1321416 1322160 1322880
##  [37] 1323624 1324368 1325040 1325784 1326504 1327248 1327968 1328712 1329456
##  [46] 1330176 1330920 1331640 1332384 1333128 1333824 1334568 1335288 1336032
##  [55] 1336752 1337496 1338240 1338960 1339704 1340424 1341168 1341912 1342584
##  [64] 1343328 1344048 1344792 1345512 1346256 1347000 1347720 1348464 1349184
##  [73] 1349928 1350672 1351344 1352088 1352808 1353552 1354272 1355016 1355760
##  [82] 1356480 1357224 1357944 1358688 1359432 1360104 1360848 1361568 1362312
##  [91] 1363032 1363776 1364520 1365240 1365984 1366704 1367448 1368192 1368888
## [100] 1369632 1370352 1371096 1371816 1372560 1373304 1374024 1374768 1375488
## [109] 1376232 1376976 1377648 1378392 1379112 1379856 1380576 1381320 1382064
## [118] 1382784 1383528 1384248 1384992 1385736 1386408 1387152 1387872 1388616
## [127] 1389336 1390080 1390824 1391544 1392288 1393008 1393752 1394496 1395168
## [136] 1395912 1396632 1397376 1398096 1398840 1399584 1400304 1401048 1401768
## [145] 1402512 1403256 1403952 1404696 1405416 1406160 1406880 1407624 1408368
## [154] 1409088 1409832 1410552 1411296 1412040 1412712 1413456 1414176 1414920
## [163] 1415640 1416384 1417128 1417848 1418592 1419312 1420056 1420800 1421472
## [172] 1422216 1422936 1423680 1424400 1425144 1425888 1426608 1427352 1428072
## [181] 1428816 1429560 1430232 1430976 1431696 1432440 1433160 1433904 1434648
## [190] 1435368 1436112 1436832 1437576 1438320 1439016 1439760 1440480 1441224
## [199] 1441944 1442688 1443432 1444152 1444896 1445616 1446360 1447104 1447776
## [208] 1448520 1449240 1449984 1450704 1451448 1452192 1452912 1453656 1454376
## [217] 1455120 1455864 1456536 1457280 1458000 1458744 1459464 1460208 1460952
## [226] 1461672 1462416 1463136 1463880 1464624 1465296 1466040 1466760 1467504
## [235] 1468224 1468968 1469712 1470432 1471176 1471896 1472640 1473384 1474080
## [244] 1474824 1475544 1476288 1477008 1477752 1478496 1479216 1479960 1480680
## [253] 1481424 1482168 1482840 1483584 1484304 1485048 1485768 1486512 1487256
## [262] 1487976 1488720 1489440 1490184 1490928 1491600 1492344 1493064 1493808
## [271] 1494528 1495272 1496016 1496736 1497480 1498200 1498944 1499688 1500360
## [280] 1501104 1501824 1502568 1503288 1504032 1504776 1505496 1506240 1506960
## [289] 1507704 1508448 1509144 1509888 1510608 1511352 1512072 1512816 1513560
## [298] 1514280 1515024 1515744 1516488 1517232 1517904 1518648 1519368 1520112
## [307] 1520832 1521576 1522320 1523040 1523784 1524504 1525248 1525992 1526664
## [316] 1527408 1528128 1528872 1529592 1530336 1531080 1531800 1532544 1533264
## [325] 1534008 1534752 1535424 1536168 1536888 1537632 1538352 1539096 1539840
## [334] 1540560 1541304 1542024 1542768 1543512 1544208 1544952 1545672 1546416
## [343] 1547136 1547880 1548624 1549344 1550088 1550808 1551552 1552296 1552968
## [352] 1553712 1554432 1555176 1555896 1556640 1557384 1558104 1558848 1559568
## [361] 1560312 1561056 1561728 1562472 1563192 1563936 1564656 1565400 1566144
## [370] 1566864 1567608 1568328 1569072 1569816 1570488 1571232 1571952 1572696
## [379] 1573416 1574160 1574904 1575624 1576368 1577088 1577832 1578576 1579272
## [388] 1580016 1580736 1581480 1582200 1582944 1583688 1584408 1585152 1585872
## [397] 1586616 1587360 1588032 1588776 1589496 1590240 1590960 1591704 1592448
## [406] 1593168 1593912 1594632 1595376 1596120 1596792 1597536 1598256 1599000
## [415] 1599720 1600464 1601208 1601928 1602672 1603392 1604136 1604880 1605552
## [424] 1606296 1607016 1607760 1608480 1609224 1609968 1610688 1611432 1612152
## [433] 1612896 1613640 1614336 1615080 1615800 1616544 1617264 1618008 1618752
## [442] 1619472 1620216 1620936 1621680 1622424 1623096 1623840 1624560 1625304
## [451] 1626024 1626768 1627512 1628232 1628976 1629696 1630440 1631184 1631856
## [460] 1632600 1633320 1634064 1634784 1635528 1636272 1636992 1637736 1638456
## [469] 1639200 1639944 1640616 1641360 1642080 1642824 1643544 1644288 1645032
## [478] 1645752 1646496 1647216 1647960 1648704 1649400 1650144 1650864 1651608
## [487] 1652328 1653072 1653816 1654536 1655280 1656000 1656744 1657488 1658160
## [496] 1658904 1659624 1660368 1661088 1661832 1662576 1663296 1664040 1664760
## [505] 1665504 1666248 1666920 1667664 1668384 1669128 1669848 1670592 1671336
## [514] 1672056 1672800 1673520 1674264 1675008 1675680 1676424 1677144 1677888
## [523] 1678608 1679352 1680096 1680816 1681560 1682280 1683024 1683768 1684464
## [532] 1685208 1685928 1686672 1687392 1688136 1688880 1689600 1690344 1691064
## [541] 1691808 1692552 1693224 1693968 1694688 1695432 1696152 1696896 1697640
## [550] 1698360 1699104 1699824 1700568 1701312 1701984 1702728 1703448 1704192
## [559] 1704912 1705656 1706400 1707120 1707864 1708584 1709328 1710072 1710744
## [568] 1711488 1712208 1712952 1713672 1714416 1715160 1715880 1716624 1717344
## [577] 1718088 1718832 1719528 1720272 1720992 1721736 1722456 1723200 1723944
## [586] 1724664 1725408 1726128 1726872 1727616 1728288 1729032 1729752 1730496
## [595] 1731216 1731960 1732704 1733424 1734168 1734888 1735632 1736376 1737048
## [604] 1737792 1738512 1739256 1739976 1740720 1741464 1742184 1742928 1743648
## [613] 1744392 1745136 1745808 1746552 1747272 1748016 1748736 1749480 1750224
## [622] 1750944 1751688 1752408 1753152 1753896 1754592 1755336 1756056 1756800
## [631] 1757520 1758264 1759008 1759728 1760472 1761192 1761936 1762680 1763352
## [640] 1764096 1764816 1765560 1766280 1767024 1767768 1768488 1769232 1769952
## [649] 1770696 1771440 1772112 1772856 1773576 1774320 1775040 1775784 1776528
## [658] 1777248 1777992 1778712 1779456 1780200 1780872 1781616 1782336 1783080
## [667] 1783800 1784544 1785288 1786008 1786752 1787472 1788216 1788960 1789656
## [676] 1790400 1791120 1791864 1792584 1793328 1794072 1794792 1795536 1796256
## [685] 1797000 1797744 1798416 1799160 1799880 1800624 1801344 1802088 1802832
## [694] 1803552 1804296 1805016 1805760 1806504 1807176 1807920 1808640 1809384
## [703] 1810104 1810848 1811592 1812312 1813056 1813776 1814520 1815264 1815936
## [712] 1816680 1817400 1818144 1818864 1819608 1820352 1821072 1821816 1822536
## [721] 1823280 1824024 1824720 1825464 1826184 1826928 1827648 1828392 1829136
## [730] 1829856 1830600 1831320 1832064 1832808 1833480 1834224 1834944 1835688
## [739] 1836408 1837152 1837896 1838616 1839360 1840080 1840824 1841568 1842240
## [748] 1842984 1843704 1844448 1845168 1845912 1846656 1847376 1848120 1848840
## [757] 1849584 1850328 1851000 1851744 1852464 1853208 1853928 1854672 1855416
## [766] 1856136 1856880 1857600 1858344 1859088 1859784 1860528 1861248 1861992
## [775] 1862712 1863456 1864200 1864920 1865664 1866384 1867128 1867872 1868544
## [784] 1869288 1870008 1870752 1871472 1872216 1872960 1873680 1874424 1875144
## [793] 1875888 1876632 1877304 1878048 1878768 1879512 1880232 1880976 1881720
## [802] 1882440 1883184 1883904 1884648 1885392 1886064 1886808 1887528 1888272
## [811] 1888992 1889736 1890480 1891200 1891944 1892664 1893408 1894152 1894848
## [820] 1895592 1896312 1897056 1897776 1898520 1899264 1899984

nc$dim$time$units

## [1] "hours since 1800-01-01 00:00:0.0"

#[1] "hours since 1800-01-01 00:00:0.0"
# nc$dim$level$vals
Lon <- ncvar_get(nc, "lon")
Lat1 <- ncvar_get(nc, "lat")
Time <- ncvar_get(nc, "time")
#Time is the same as nc$dim$time$vals
head(Time)

## [1] 1297320 1298064 1298760 1299504 1300224 1300968

#[1] 1297320 1298064 1298760 1299504 1300224 1300968
library(chron)

## Warning: package 'chron' was built under R version 4.1.3

Tymd <- month.day.year(Time[1]/24,c(month = 1, day = 1, year = 1800))
#c(month = 1, day = 1, year = 1800) is the reference time
Tymd

## $month
## [1] 1
## 
## $day
## [1] 1
## 
## $year
## [1] 1948

#$month
#[1] 1
#$day
#[1] 1
#$year
#[1] 1948
#1948-01-01
precnc <- ncvar_get(nc, "air")
dim(precnc)

## [1] 144  73 826

#[1] 144  73 826, i.e., 826 months=1948-01 to 2016-10, 68 years 10 mons

 

Plot Fig. 10.5

#plot a 90S-90N temp along a meridional line at 180E
par(mar = c(4,4,2.5,4))
plot(seq(90,-90,length = 73),precnc[72,,1], type = "o", lwd = 2,
     xlab = "Latitude", ylab = "Temperature [°C]",
     main = "NCEP/NCAR Reanalysis Surface Air Temperature [°C]\nAlong a Meridional Line at 180°E: Jan. 1948")

#Write data as a CSV space-time matrix with a header
precst <- matrix(0,nrow = 10512,ncol = 826)
temp <- as.vector(precnc[,,1])
head(temp)

## [1] -34.92677 -34.92677 -34.92677 -34.92677 -34.92677 -34.92677

for (i in 1:826) {precst[,i] = as.vector(precnc[ , , i])}
dim(precst)

## [1] 10512   826

#[1] 10512   826
#Build lat and lon for 10512 spatial positions using rep
LAT <- rep(Lat1, 144)
LON <- rep(Lon[1],each = 73)
gpcpst <- cbind(LAT, LON, precst)
dim(gpcpst)

## [1] 10512   828

#[1] 10512   828
#The first two columns are lat and lon. 826 months: 1948.01-2016.10
#Convert the Julian day and hour into calendar months for time 
tm <- month.day.year(Time/24, c(month = 1, day = 1, year = 1800))
tm1 <- paste(tm$year,"-",tm$month)
tm2 <- c("Lat","Lon",tm1) #This is the header
#Assign the header to the space-time data matrix
colnames(gpcpst) <- tm2
#setwd("~/sshen/climmath/data")
#setwd routes the desired csv file to a given directory
write.csv(gpcpst,file = "ncepJan1948_Oct2016.csv")

 

Compute the January climatology and standard deviation: 1948-2015

monJ <- seq(1,816,12)
gpcpdat <- gpcpst[,3:818]
gpcpJ <- gpcpdat[,monJ]
climJ <- rowMeans(gpcpJ) #Use all 68 January's from 1948 to 2015
library(matrixStats)# rowSds command is in the matrixStats package

## Warning: package 'matrixStats' was built under R version 4.1.3

sdJ <- rowSds(gpcpJ)

 

Plot Fig. 9.9

#Verify previous calculations and plots in Fig. 9.9
#Plot Jan. climatology
Lat <- -Lat1
mapmat <- matrix(climJ,nrow = 144)
mapmat<- mapmat[,seq(length(mapmat[1,]),1)]
plot.new()
int <- seq(-50,50,length.out = 81)
rgb.palette <- colorRampPalette(c('black','blue', 'darkgreen','green', 
                               'white','yellow','pink','red','maroon'),interpolate = 'spline')
library(maps)

## Warning: package 'maps' was built under R version 4.1.3

filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "NCEP Jan. SAT RA 1948-2015 Climatology [°C]",
                                xlab = "Longitude",ylab = "Latitude",
                                cex.lab = 1.2,cex.axis = 1.2),
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "[°C]"))

 

Plot Jan Standard Deviation

Lat = -Lat1
mapmat = matrix(sdJ,nrow = 144)
mapmat = pmax(pmin(mapmat,6),0)
mapmat = mapmat[,seq(length(mapmat[1,]),1)]
plot.new()
int = seq(0,6,length.out = 91)
rgb.palette = colorRampPalette(c('black','blue', 'green', 
                               'yellow', 'pink','red','maroon'),interpolate = 'spline')
filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "NCEP Jan. SAT RA 1948-2015 Standard Deviation [°C]", 
               cex.lab = 1.2,cex.axis = 1.2, xlab = "Longitude",ylab = "Latitude"),
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "[°C]"))

 

Plot Fig. 10.6 (Eigenvalues and variances explained)

#Compute the Jan. EOFs
monJ <- seq(1,816,12)
gpcpdat <- gpcpst[,3:818]
gpcpJ <- gpcpdat[,monJ]
climJ <- rowMeans(gpcpJ)
library(matrixStats)
sdJ <- rowSds(gpcpJ)
anomJ <- (gpcpdat[,monJ]-climJ)/sdJ #standardized anomalies
anomAW <- sqrt(cos(gpcpst[,1]*pi/180))*anomJ #Area weighted anomalies
svdJ <- svd(anomAW)  #execute SVD

#plot eigenvalues
par(mar = c(4,4,2,4))
plot(100*(svdJ$d)^2/sum((svdJ$d)^2), type = "o", ylab = "Percentage of Variance [%]",
     xlab = "Mode Number", main = "Eigenvalues of Covariance Matrix")
legend(20,5, col = c("black"),lty = 1,lwd = 2.0,
       legend = c("Percentage Variance"),bty = "n",
       text.font = 2,cex = 1.2, text.col = "black")
par(new = TRUE)
plot(cumsum(100*(svdJ$d)^2/sum((svdJ$d)^2)),type = "o",
     col = "blue",lwd = 1.5,axes = FALSE,xlab = "",ylab = "")
legend(20,50, col = c("blue"),lty = 1,lwd = 2.0,
       legend = c("Cumulative Variance"),bty = "n",
       text.font = 2,cex = 1.2, text.col = "blue")
axis(4)
mtext("Cumulative Variance [%]",side = 4,line = 2)

 

Plot Fig. 10.7 (EOFs and PCs)

#plot EOF1: The physical EOF = eigenvector divided by area factor
mapmat <- matrix(svdJ$u[,1]/sqrt(cos(gpcpst[,1]*pi/180)),nrow = 144)
rgb.palette = colorRampPalette(c('blue','green','white', 
                               'yellow','red'),interpolate = 'spline')
int <- seq(-0.04,0.04,length.out = 61)
mapmat <- mapmat[, seq(length(mapmat[1,]),1)]
filled.contour(Lon, Lat, -mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "January EOF1 from 1948-2015 NCEP Temp Data",
               xlab = "Longitude",ylab = "Latitude",cex.lab = 1.2,cex.axis = 1.2),
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "Scale"))

#
#plot PC1
pcdat <- svdJ$v[,1]
Time <- seq(1948,2015)
plot(Time, -pcdat, type = "o", main = "PC1 of NCEP RA Jan. SAT: 1948-2015",
     xlab = "Year",ylab = "PC Values",cex.lab = 1.2,cex.axis = 1.2,
     lwd = 2, ylim = c(-0.3,0.3))

 

Plot Fig. 10.8

#plot EOF2: The physical EOF =  eigenvector divided by area factor
mapmat <- matrix(svdJ$u[,2]/sqrt(cos(gpcpst[,1]*pi/180)),nrow = 144)
rgb.palette <- colorRampPalette(c('blue','green','white', 
                               'yellow','red'),interpolate = 'spline')
int <- seq(-0.04,0.04,length.out = 61)
mapmat <- mapmat[, seq(length(mapmat[1,]),1)]
filled.contour(Lon, Lat, -mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "January EOF2 from 1948-2015 NCEP Temp. Data",
               xlab = "Longitude",ylab = "Latitude",cex.lab = 1.2,cex.axis = 1.2),
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "Scale"))

#plot PC2
pcdat <- svdJ$v[,2]
Time <- seq(1948,2015)
plot(Time, -pcdat, type = "o", main = "PC2 of NCEP RA Jan. SAT: 1948-2015",
     xlab = "Year",ylab = "PC Values",cex.lab = 1.2,cex.axis = 1.2,
     lwd = 2, ylim = c(-0.3,0.3))

 

Plot Fig. 10.9

#plot EOF3: The physical EOF = eigenvector divided by area factor
mapmat <- matrix(svdJ$u[,3]/sqrt(cos(gpcpst[,1]*pi/180)),nrow = 144)
rgb.palette <- colorRampPalette(c('blue','green','white', 
                               'yellow','red'),interpolate = 'spline')
int <- seq(-0.04,0.04,length.out = 61)
mapmat <- mapmat[, seq(length(mapmat[1,]),1)]
filled.contour(Lon, Lat, -mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "January EOF3 from 1948-2015 NCEP Temp. Data",
               xlab = "Longitude",ylab = "Latitude", cex.lab  =  1.2,cex.axis = 1.2),
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "Scale"))

#plot PC3
pcdat <- svdJ$v[,3]
Time <- seq(1948,2015)
plot(Time, -pcdat, type = "o", main = "PC3 of NCEP RA Jan. SAT: 1948-2015",
     xlab = "Year",ylab = "PC Values",cex.lab = 1.2,cex.axis = 1.2,
     lwd = 2, ylim = c(-0.3,0.3))

 

Plot Fig. 10.10 (EOFs and PCs from de-trended standardized data)

#EOF from de-trended data
monJ <- seq(1,816,12)
gpcpdat <- gpcpst[,3:818]
gpcpJ <- gpcpdat[,monJ]
climJ <- rowMeans(gpcpJ)
library(matrixStats)
sdJ <- rowSds(gpcpJ)
anomJ <- (gpcpdat[,monJ]-climJ)/sdJ
trendM <- matrix(0,nrow = 10512, ncol = 68)#trend field matrix
trendV <- rep(0,len = 10512)#trend for each grid box: a vector
for (i in 1:10512) {
  trendM[i,] <- (lm(anomJ[i,] ~ Time))$fitted.values
  trendV[i] <- lm(anomJ[i,] ~ Time)$coefficients[2]
}
dtanomJ <- anomJ - trendM
dim(dtanomJ)

## [1] 10512    68

dtanomAW <- sqrt(cos(gpcpst[,1]*pi/180))*dtanomJ
svdJ <- svd(dtanomAW)

#Plot EOF1
mapmat <- matrix(svdJ$u[,1]/sqrt(cos(gpcpst[,1]*pi/180)),nrow = 144)
rgb.palette <- colorRampPalette(c('blue','green','white', 
                               'yellow','red'),interpolate = 'spline')
int <- seq(-0.04,0.04,length.out = 61)
mapmat <- mapmat[, seq(length(mapmat[1,]),1)]
#par(mar = c(4,4,3,4))
filled.contour(Lon, Lat, -mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(
               main = "January EOF1 from 1948-2015 NCEP De-Trended\nStandardized Temp. Anomaly Data",
               xlab = "Longitude",ylab = "Latitude",cex.lab = 1.2,cex.axis = 1.2),
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "Scale"))

#plot PC1
pcdat <- svdJ$v[,1]
Time <- seq(1948,2015)
par(mar = c(4,4,3,0.5))
plot(Time, -pcdat, type = "o", 
     main = "PC1 of NCEP RA Jan. SAT: 1948-2015\nDe-Trended Standardized Data",
     xlab = "Year",ylab = "PC Values",cex.lab = 1.2,cex.axis = 1.2,
     lwd = 2, ylim = c(-0.3,0.3))

 

Plot Fig. 10.11 (EOF2 and PC2 of de-trended standardized data)

#Plot EOF2
mapmat <- matrix(svdJ$u[,2]/sqrt(cos(gpcpst[,1]*pi/180)),nrow = 144)
rgb.palette <- colorRampPalette(c('blue','green','white', 
                               'yellow','red'),interpolate = 'spline')
int <- seq(-0.04,0.04,length.out = 61)
mapmat <- mapmat[, seq(length(mapmat[1,]),1)]
#par(mar = c(4,4,3,4))
filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(
               main = "January EOF2 from 1948-2015 NCEP\nDe-Trended Standardized Temp. Anomaly Data",
               xlab = "Longitude",ylab = "Latitude",cex.lab = 1.2,cex.axis = 1.2),
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "Scale"))

#
#plot PC2
pcdat <- svdJ$v[,2]
Time <- seq(1948,2015)
par(mar = c(4,4,3,0.5))
plot(Time, pcdat, type = "o", 
     main = "PC2 of NCEP RA Jan. SAT: 1948-2016\nDe-Trended Standardized Data",
     xlab = "Year",ylab = "PC Values",cex.lab = 1.2,cex.axis = 1.2,
     lwd = 2, ylim = c(-0.3,0.3))

 

Plot Fig. 10.12 (Global average January SAT and its trend)

#Plot the area-weighted global average Jan temp from 1948-2015
#Begin from the space-time data matrix gpcpst[,1]
vArea <- cos(gpcpst[,1]*pi/180)
anomA <- vArea*anomJ
dim(anomA)

## [1] 10512    68

JanSAT <- colSums(anomA)/sum(vArea)
plot(Time, JanSAT, type = "o", lwd = 2,
     main = "Global Average Jan. SAT Anomalies from NCEP RA",
     xlab = "Year",ylab = "Temperature [°C]")
regSAT <- lm(JanSAT ~ Time)
regSAT

## 
## Call:
## lm(formula = JanSAT ~ Time)
## 
## Coefficients:
## (Intercept)         Time  
##   -9.494035     0.004791

#0.48°C/100a trend
abline(regSAT, col = "red", lwd = 4)
text(1962,0.35,"Linear trend 0.48°C/Century", col = "red", cex = 1.2)

 

Plot Fig. 10.13 (Temperature trend map)

#plot the trend of Jan SAT non-standardized anomaly data 
#Begin with the space-time data matrix
monJ <- seq(1,816,12)
gpcpdat <- gpcpst[,3:818]
gpcpJ <- gpcpdat[,monJ]
plot(gpcpJ[,23],cex.lab = 1.,cex.axis = 1.)

climJ <- rowMeans(gpcpJ)
anomJ <- (gpcpdat[,monJ]-climJ)
trendV <- rep(0,len = 10512)#trend for each grid box: a vector
for (i in 1:10512) {
  trendV[i] <- lm(anomJ[i,] ~ Time)$coefficients[2]
}

mapmat1 <- matrix(10*trendV,nrow = 144)
mapv1 <- pmin(mapmat1,1) #Compress the values >5 to 5
mapmat <- pmax(mapv1,-1) #compress the values <-5 t -5
rgb.palette <- colorRampPalette(c('blue','green','white', 'yellow','red'),
                             interpolate = 'spline')
int <- seq(-1,1,length.out = 61)
mapmat <- mapmat[, seq(length(mapmat[1,]),1)]
#par(mar = c(4,4,2,4))
par(mgp = c(2,1,0), cex.lab = 1.2, cex.axis = 1.2)
filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "Trend of the NCEP RA1 Jan 1948-2015 Temp. Anomaly",
               xlab = "Longitude",ylab = "Latitude"),cex.lab = 1.2,cex.axis = 1.2,
               plot.axes = {axis(1); axis(2);map('world2', add = TRUE);grid()},
               key.title = title(main = "°C/10a"))

 

Plot Fig. 10.14 (GPCP precipitation from 90S to 90N)

#Download the GPCP data from the ESRL website
#https://www.esrl.noaa.gov/psd/data/gridded/data.gpcp.html
#Select the monthly precipitation data by clicking 
#the Download File: precip.mon.mean.nc
#Move the data file to your assigned working directory 
#/Users/sshen/climmath/data
setwd("~/sshen/climmath")
#install.packages("ncdf")
library(ncdf4)
# 4 dimensions: lon,lat,level,time
nc <- ncdf4::nc_open("data/precip.mon.mean.nc")
#nc
nc$dim$lon$vals

##   [1]   1.25   3.75   6.25   8.75  11.25  13.75  16.25  18.75  21.25  23.75
##  [11]  26.25  28.75  31.25  33.75  36.25  38.75  41.25  43.75  46.25  48.75
##  [21]  51.25  53.75  56.25  58.75  61.25  63.75  66.25  68.75  71.25  73.75
##  [31]  76.25  78.75  81.25  83.75  86.25  88.75  91.25  93.75  96.25  98.75
##  [41] 101.25 103.75 106.25 108.75 111.25 113.75 116.25 118.75 121.25 123.75
##  [51] 126.25 128.75 131.25 133.75 136.25 138.75 141.25 143.75 146.25 148.75
##  [61] 151.25 153.75 156.25 158.75 161.25 163.75 166.25 168.75 171.25 173.75
##  [71] 176.25 178.75 181.25 183.75 186.25 188.75 191.25 193.75 196.25 198.75
##  [81] 201.25 203.75 206.25 208.75 211.25 213.75 216.25 218.75 221.25 223.75
##  [91] 226.25 228.75 231.25 233.75 236.25 238.75 241.25 243.75 246.25 248.75
## [101] 251.25 253.75 256.25 258.75 261.25 263.75 266.25 268.75 271.25 273.75
## [111] 276.25 278.75 281.25 283.75 286.25 288.75 291.25 293.75 296.25 298.75
## [121] 301.25 303.75 306.25 308.75 311.25 313.75 316.25 318.75 321.25 323.75
## [131] 326.25 328.75 331.25 333.75 336.25 338.75 341.25 343.75 346.25 348.75
## [141] 351.25 353.75 356.25 358.75

nc$dim$lat$vals

##  [1] -88.75 -86.25 -83.75 -81.25 -78.75 -76.25 -73.75 -71.25 -68.75 -66.25
## [11] -63.75 -61.25 -58.75 -56.25 -53.75 -51.25 -48.75 -46.25 -43.75 -41.25
## [21] -38.75 -36.25 -33.75 -31.25 -28.75 -26.25 -23.75 -21.25 -18.75 -16.25
## [31] -13.75 -11.25  -8.75  -6.25  -3.75  -1.25   1.25   3.75   6.25   8.75
## [41]  11.25  13.75  16.25  18.75  21.25  23.75  26.25  28.75  31.25  33.75
## [51]  36.25  38.75  41.25  43.75  46.25  48.75  51.25  53.75  56.25  58.75
## [61]  61.25  63.75  66.25  68.75  71.25  73.75  76.25  78.75  81.25  83.75
## [71]  86.25  88.75

nc$dim$time$vals

##   [1] 65378 65409 65437 65468 65498 65529 65559 65590 65621 65651 65682 65712
##  [13] 65743 65774 65803 65834 65864 65895 65925 65956 65987 66017 66048 66078
##  [25] 66109 66140 66168 66199 66229 66260 66290 66321 66352 66382 66413 66443
##  [37] 66474 66505 66533 66564 66594 66625 66655 66686 66717 66747 66778 66808
##  [49] 66839 66870 66898 66929 66959 66990 67020 67051 67082 67112 67143 67173
##  [61] 67204 67235 67264 67295 67325 67356 67386 67417 67448 67478 67509 67539
##  [73] 67570 67601 67629 67660 67690 67721 67751 67782 67813 67843 67874 67904
##  [85] 67935 67966 67994 68025 68055 68086 68116 68147 68178 68208 68239 68269
##  [97] 68300 68331 68359 68390 68420 68451 68481 68512 68543 68573 68604 68634
## [109] 68665 68696 68725 68756 68786 68817 68847 68878 68909 68939 68970 69000
## [121] 69031 69062 69090 69121 69151 69182 69212 69243 69274 69304 69335 69365
## [133] 69396 69427 69455 69486 69516 69547 69577 69608 69639 69669 69700 69730
## [145] 69761 69792 69820 69851 69881 69912 69942 69973 70004 70034 70065 70095
## [157] 70126 70157 70186 70217 70247 70278 70308 70339 70370 70400 70431 70461
## [169] 70492 70523 70551 70582 70612 70643 70673 70704 70735 70765 70796 70826
## [181] 70857 70888 70916 70947 70977 71008 71038 71069 71100 71130 71161 71191
## [193] 71222 71253 71281 71312 71342 71373 71403 71434 71465 71495 71526 71556
## [205] 71587 71618 71647 71678 71708 71739 71769 71800 71831 71861 71892 71922
## [217] 71953 71984 72012 72043 72073 72104 72134 72165 72196 72226 72257 72287
## [229] 72318 72349 72377 72408 72438 72469 72499 72530 72561 72591 72622 72652
## [241] 72683 72714 72742 72773 72803 72834 72864 72895 72926 72956 72987 73017
## [253] 73048 73079 73108 73139 73169 73200 73230 73261 73292 73322 73353 73383
## [265] 73414 73445 73473 73504 73534 73565 73595 73626 73657 73687 73718 73748
## [277] 73779 73810 73838 73869 73899 73930 73960 73991 74022 74052 74083 74113
## [289] 74144 74175 74203 74234 74264 74295 74325 74356 74387 74417 74448 74478
## [301] 74509 74540 74569 74600 74630 74661 74691 74722 74753 74783 74814 74844
## [313] 74875 74906 74934 74965 74995 75026 75056 75087 75118 75148 75179 75209
## [325] 75240 75271 75299 75330 75360 75391 75421 75452 75483 75513 75544 75574
## [337] 75605 75636 75664 75695 75725 75756 75786 75817 75848 75878 75909 75939
## [349] 75970 76001 76030 76061 76091 76122 76152 76183 76214 76244 76275 76305
## [361] 76336 76367 76395 76426 76456 76487 76517 76548 76579 76609 76640 76670
## [373] 76701 76732 76760 76791 76821 76852 76882 76913 76944 76974 77005 77035
## [385] 77066 77097 77125 77156 77186 77217 77247 77278 77309 77339 77370 77400
## [397] 77431 77462 77491 77522 77552 77583 77613 77644 77675 77705 77736 77766
## [409] 77797 77828 77856 77887 77917 77948 77978 78009 78040 78070 78101 78131
## [421] 78162 78193 78221 78252 78282 78313 78343 78374 78405 78435 78466 78496
## [433] 78527 78558 78586 78617 78647 78678 78708 78739 78770 78800 78831 78861
## [445] 78892 78923 78952 78983 79013 79044 79074

#nc$dim$time$units
#nc$dim$level$vals
Lon <- ncvar_get(nc, "lon")
Lat <- ncvar_get(nc, "lat")
Time <- ncvar_get(nc, "time")
head(Time)

## [1] 65378 65409 65437 65468 65498 65529

#[1] 65378 65409 65437 65468 65498 65529
library(chron)
month.day.year(65378,c(month = 1, day = 1, year = 1800))

## $month
## [1] 1
## 
## $day
## [1] 1
## 
## $year
## [1] 1979

#1979-01-01
precnc <- ncvar_get(nc, "precip")
dim(precnc)

## [1] 144  72 451

#[1] 144  72 451 #2.5-by-2.5, 451 months from Jan 1979-July 2016
#plot a 90S-90N precip along a meridional line at 160E
par(mar = c(4.5,4.5,2,0.5))
plot(seq(-90,90,length = 72), precnc[64,,1], type = "l", lwd = 2,
     xlab = "Latitude", ylab = "Precipitation [mm/day]",
     main = "90°S-90°N Precipitation Along the Meridian at 160°E: Jan. 1979",
     cex.lab = 1.2, cex.axis = 1.2)

 

Plot Fig. 10.15 (GPCP climatology and standard deviation maps)

#Write the data as space-time matrix with a header
precst <- matrix(0,nrow = 10368,ncol = 451)
temp <- as.vector(precnc[,,1])
head(temp)

## [1] 0 0 0 0 0 0

for (i in 1:451) {precst[,i] <- as.vector(precnc[ , , i])}
#precst is the space-time GPCP data 
LAT <- rep(Lat, 144)
LON <- rep(Lon[1],72)
for (i in 2:144){LON <- c(LON, rep(Lon[i],72))}
gpcpst <- cbind(LAT, LON, precst)

#Convert the Julian day into calender dates for time 
tm <- month.day.year(Time, c(month = 1, day = 1, year = 1800))
tm1 <- paste(tm$year,"-",tm$month)
#tm1 <- data.frame(tm1)
tm2 <- c("Lat","Lon",tm1)
colnames(gpcpst) <- tm2
#Look at a sample section of the space-time data 
gpcpst[890:892,1:5]

##         Lat   Lon   1979 - 1  1979 - 2  1979 - 3
## [1,] -26.25 31.25 0.05752099 0.1805309 0.2104454
## [2,] -23.75 31.25 0.10855579 0.2562788 0.2883888
## [3,] -21.25 31.25 0.08718992 0.2738046 0.2569866

#        Lat   Lon   1979 - 1  1979 - 2  1979 - 3
#[1,] -26.25 31.25 0.05752099 0.1805309 0.2104454
#[2,] -23.75 31.25 0.10855579 0.2562788 0.2883888
#[3,] -21.25 31.25 0.08718992 0.2738046 0.2569866
write.csv(gpcpst,file = "gpcp9716jul.csv")
#
#Compute and plot the GPCP climatology from Jan 1979-July 2016
library(maps)
climmat <- precnc[,,1]
for(i in 2:451){climmat <- climmat + precnc[,,i]}
climmat <- climmat/451
mapmat <- climmat
mapmat <- pmax(pmin(mapmat,10),0)
int <- seq(min(mapmat),max(mapmat),length.out = 11)
rgb.palette <- colorRampPalette(c('bisque2','cyan', 'green', 'yellow',
                               'pink','indianred2', 'red','maroon','black'),
                             interpolate = 'spline')
#plot.new()
#par(mar = c(4, 4.5, 2, 1))
filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "GPCP 1979-2016 Precipitation Climatology [mm/day]",
                                xlab = "Longitude",ylab = "Latitude", cex.lab = 1.2),
               plot.axes = {axis(1,seq(0,360, by = 30), cex.axis = 1.2); 
                 axis(2, seq(-90,90,by = 30), cex.axis = 1.2);
                 map('world2', add = TRUE);grid()},
               key.title = title(main = "mm/day"),
               key.axes = {axis(4, seq(0,10, len = 11), cex.axis = 1.2)})

#Compute and plot standard deviation from Jan 1979-July 2016
sdmat <- (precnc[,,1]-climmat)^2
for(i in 2:451){sdmat <- sdmat + (precnc[,,i]-climmat)^2}
sdmat <- sqrt(sdmat/451)
mapmat <- sdmat
mapmat <- pmax(pmin(mapmat,5),0)
int <- seq(min(mapmat),max(mapmat),length.out = 21)
rgb.palette <- colorRampPalette(c('bisque2','cyan', 'green', 'yellow','pink',
                               'indianred2', 'red','maroon','black'),
                             interpolate = 'spline')
#plot.new()
#par(mar = c(4.5, 4.5, 3, 1))
filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(
               main = "GPCP 1979-2016 Standard Deviation of\nPrecipitation [mm/day]",
               xlab = "Longitude",ylab = "Latitude", cex.lab = 1.2),
               plot.axes = {axis(1,seq(0,360, by = 30), cex.axis = 1.2); 
               axis(2,seq(-90,90,by = 30), cex.axis = 1.2);
               map('world2', add = TRUE);grid()},
               key.title = title(main = "mm/day"),
               key.axes = {axis(4, seq(0,5,len = 11), cex.axis = 1.2)})

#Compute the January precipitation climatology and plot its map
Jmon <- seq(3,453,by = 12)
Jclim <- rowMeans(gpcpst[,Jmon])
mapmat <- matrix(Jclim,nrow = 144)
mapmat <- pmax(pmin(mapmat,10),0)
int <- seq(min(mapmat),max(mapmat),length.out = 11)
rgb.palette <- colorRampPalette(c('bisque2','cyan', 'green', 'yellow','pink','indianred2','red','maroon','black'), interpolate = 'spline')
#plot.new()
#par(mar = c(4, 4.5, 2, 1))
filled.contour(Lon, Lat, mapmat, color.palette = rgb.palette, levels = int,
               plot.title = title(main = "GPCP 1979-2016 Precipitation Jan. Climatology [mm/day]",
               xlab = "Longitude",ylab = "Latitude", cex.lab = 1.2),
               plot.axes = {axis(1,seq(0,360, by = 30), cex.axis = 1.2); 
               axis(2, seq(-90,90,by = 30), cex.axis = 1.2);
               map('world2', add = TRUE);grid()},
               key.title = title(main = "mm/day"),
               key.axes = {axis(4, seq(0,10, len = 11), cex.axis = 1.2)})