Introduction to Modern Mathematical Modeling with R:

A User’s Manual to Train Mathematical Consultants

 

A Cambridge University Press book by

SSP Shen

 

 

R Code written by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
https://shen.sdsu.edu
Email:

 

Compiled and Edited by Joaquin Stawsky
San Diego State University, August 2022

 

 

 

Chapter 9: Visualize Mathematical Models By R

Plot Two Different Time Series on the Same Plot

#Plot US temp and prec. times series on the same figure
#plot.new()
par(mar = c(4.2,4.2,2.5,4.1))
Time <- 2001:2010
Tmean <- c(12.06, 11.78,11.81,11.72,12.02,
           12.36,12.03,11.27,11.33,11.66)
Prec <- c(737.11,737.87,774.95,844.55,
          764.03,757.43,741.17,793.50,820.42,796.80)

plot(Time,Tmean,
     type="o", col="red",
     xlab="Year", ylab=expression(paste(T[mean], " [", degree,"C]")),lwd=1.5,
     main="Contiguous U.S. Annual Mean Temperature \nand Total Precipitation")
grid(nx = NULL, ny = NULL)
legend(2000.5,12.42, 
       col=c("red"),lty=1,lwd=2.0, 
       legend=c(expression(T[mean])),bty="n",
       text.font=2,cex=1.0)

#Allows a figure to be overlaid on the first plot
par(new=TRUE)
plot(Time, Prec, 
     type="o",col="blue",
     lwd=1.5,axes=FALSE,xlab="",ylab="") 
legend(2000.5,839, col=c("blue"),lty=1,lwd=2.0,
       legend=c("Precipitation"),bty="n",text.font=2,cex=1.0)

#Suppress the axes and assign the y-axis to side 4
axis(4)
mtext("Precipitation [mm]",side=4,line=2) 

#legend("topleft",col=c("red","blue"),lty=1,legend=c("Tmean","Prec"),cex=0.6)
#Plot two legends at the same time make it difficult to adjust the font size
#because of different scale

 

Figure setups: Margins, Fonts, Mathematical Symbols, and More

#Margins, math symbol, and figure setups
#plot.new()
par(mar=c(6,4,3,3))
x <- 0.25*(-30:30)
y <- sin(x)
x1 <- x[which(sin(x) >=0)]
y1 <- sin(x1)
x2 <- x[which(sin(x) < 0)]
y2 <- sin(x2)
plot(x1,y1,xaxt="n", xlab="",ylab="",lty=1,type="h",
     lwd=3, tck=-0.02, ylim=c(-1,1), col="red",
     col.lab="purple",cex.axis=1.4) 

lines(x2,y2,xaxt="n", xlab="",ylab="",lty=3,type="h",
      col="blue",lwd=8, tck=-0.02)

axis(1, at=seq(-6,6,2),line=3, cex.axis=1.8)
axis(4, at=seq(-1,1,0.5), lab=c("A", "B", "C", "D","E"),
     cex.axis=2,las=2)

text(0,0.7,font=3,cex=6, "Sine waves", col="darkgreen")
#Itatlic font size 2

mtext(side=2,line=2, 
      expression(y==sin(theta-hat(phi))),cex=1.5, col="blue")
mtext(font=2,
      "Text outside of the figure on side 3",side=3,line=1, cex=1.5)
#Bold font

mtext(font=1, side=1,line=1, 
      expression(paste("Angle in radian: ",
                       theta-phi[0])),cex=1.5, col="red")

par(mar=c(8,6,3,2))
par(mgp=c(2.5,1,0))

plot(1:200/20, rnorm(200),sub="Sub-title: 200 Random Values",
     xlab= "Time", ylab="Random Values", main="Normal Random Values", 
     cex.lab=2.0, cex.axis=1.75, cex.main=2.5, cex.sub=1.5)
grid(nx = NULL, ny = NULL)

par(mgp=c(2,1,0))
plot(sin,xlim=c(10,20))
grid(nx = NULL, ny = NULL)

#A fancy plot of the NOAAGlobalTemp time series
#setwd("~/sshen/mathmodel")
NOAATemp <- read.table('~/Desktop/RMathModel/data/gl_land_oceanHansen.txt', 
                      header = TRUE)
plot.new()
par(mar=c(4,4,3,1))
x <- NOAATemp[,1]
y <- NOAATemp[,2]
z <- rep(-99,length(x))

for (i in 3:length(x)-2){ 
  z[i]<-mean(c(y[i-2],y[i-1],y[i],y[i+1],y[i+2]))
}

n1 <- which(y >= 0)
x1 <- x[n1]
y1 <- y[n1]
n2 <- which(y < 0)
x2 <- x[n2]
y2 <- y[n2]
x3 <- x[2:length(x)-2]
y3 <- z[2:length(x)-2] 

plot(x1,y1,
     type="h",xlim=c(1880,2016),
     lwd=3, tck=0.02, ylim=c(-0.7,0.7),
     ylab="Temperature [°C]", xlab="Time",col="red",
     main="NOAA Global Average Annual Mean Temperature Anomalies")
grid(nx = NULL, ny = NULL)
lines(x2,y2,type="h",lwd=3, tck=-0.02, col="blue") 
lines(x3,y3,lwd=2)

 

Plot Two or More Panels on the Same Figure

#Plot US temp. and prec. times series on the same figure
par(mfrow=c(2,1))
par(mar=c(0,5,3,1)) #Zero space between (a) and (b)
Time <- 2001:2010
Tmean <- c(12.06, 11.78,11.81,11.72,12.02,
           12.36,12.03,11.27,11.33,11.66) 
Prec <- c(737.11,737.87,774.95,844.55,764.03,
          757.43,741.17,793.50,820.42,796.80)
plot(Time,Tmean,
     type="o",col="red",xaxt="n", xlab="",
     ylab=expression(paste(T[mean]," [",degree,"C]"))) 
grid(nx = NULL, ny = NULL)
text(2006, 12,font=2,"US Annual Mean Temperature", cex=1.5) 
text(2001.5,12.25,"(a)")
#Plot the panel on row 2
par(mar=c(3,5,0,1))

plot(Time, Prec,
     type="o",col="blue",
     xlab="Time",ylab="Prec. [mm]") 
grid(nx = NULL, ny = NULL)
text(2006, 800, font=2, 
     "US Annual Total Precipitation", cex=1.5) 
text(2001.5,840,"(b)")

rm(list=ls()) 
#plot.new()
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE), 
       widths=c(3,3),heights=c(2,2))
plot(sin,type="l", xlim=c(0,20)) 
plot(sin,xlim=c(0,10)) 
plot(sin,xlim=c(10,20))

 

Basic Principles for an R Contour Plot

x <- y <- seq(-1, 1, len=25)
z <- matrix(rnorm(25*25),nrow=25)

par(mar=c(3,3,3,3))
contour(x,y,z, 
        main="Contour Plot of Normal Random Values") 

filled.contour(x,y,z, 
      main="Filled Contour Plot of Normal Random Values") 

filled.contour(x,y,z, color.palette = heat.colors, 
      main="Filled Contour Plot of Normal Random Values \nwith Set Colors")

filled.contour(x,y,z, 
      color.palette = colorRampPalette(c("red", "white", "blue")), 
      main="Filled Contour Plot of Normal Random Values \nwith Specific Colors")

 

Plot Contour Color Maps for Random Values on a Map

## Figure 9.6
#Plot a 5-by-5 grid global map of standard normal random values 
library(maps)
#plot.new()
#Step 1: Generate a 5-by-5 grid (pole-to-pole, lon 0 to 355) 
Lat <- seq(-90,90,length=37) 
#Must increasing 

Lon <- seq(0,355,length=72) 
#Must increasing

#Generate the random values
mapdat <- matrix(rnorm(72*37),nrow=72)
#The matrix uses lon as row going and lat as column
#Each row includes data from south to north
#Define color
int <- seq(-3,3,length.out=81) 
rgb.palette <- colorRampPalette(c('black','purple','blue','white','green', 
                    'yellow','pink','red','maroon'), interpolate='spline')

#Plot the values on the world map
filled.contour(Lon, Lat, mapdat, color.palette=rgb.palette, 
               levels=int, 
               plot.title=title(xlab="Longitude", ylab="Latitude",
               main="Standard Normal Random Values on a World Map:\n5° Lat-Lon Grid"), 
               plot.axes={ axis(1); axis(2);map('world2', add=TRUE);grid()} )

#filled.contour() is a contour plot on an x-y grid. 
#Background maps are added later in plot.axes={} 
#axis(1) means ticks on the lower side
#axis(2) means ticks on the left side
#Save image with width=800, maintain aspect ratio

 

Figure 9.7

#Plot a 5-by-5 grid regional map to cover USA and Canada 
Lat3 <- seq(10,70,length=13)
Lon3 <- seq(230,295,length=14) 
mapdat <- matrix(rnorm(13*14),nrow=14) 
int <- seq(-3,3,length.out=81) 
rgb.palette <- colorRampPalette(c('black','purple','blue','white','green',
                    'yellow','pink','red','maroon'), interpolate='spline')

filled.contour(Lon3, Lat3, mapdat, 
               color.palette=rgb.palette, 
               levels=int, plot.title=title(
               main="Standard Normal Random Values on\n a World Map: 5° Lat-Lon Grid",
               xlab="Longitude", ylab="Latitude"),
               plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()})

 

Read .nc File

#R plot of NCEP/NCAR Reanalysis PSD monthly temp data .nc file 
#http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.surface.html
rm(list=ls(all=TRUE)) 
#setwd("~/sshen/mathmodel")
# Download netCDF file
# Library
#install.packages("ncdf") 
library(ncdf4)
# 4 dimensions: lon,lat,level,time 
nc <- ncdf4::nc_open("~/Desktop/RMathModel/data/air.mon.mean.nc")
nc
## File ~/Desktop/RMathModel/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 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 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
##   [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
#nc$dim$level$vals
Lon <- ncvar_get(nc, "lon")
Lat1 <- ncvar_get(nc, "lat")
Time <- ncvar_get(nc, "time")
head(Time)
## [1] 1297320 1298064 1298760 1299504 1300224 1300968
library(chron)
month.day.year(1297320/24,
            c(month = 1, day = 1, year = 1800)) 
## $month
## [1] 1
## 
## $day
## [1] 1
## 
## $year
## [1] 1948
precnc <- ncvar_get(nc, "air")
dim(precnc) #826 months=1948-01 to 2016-10, 68 years 10 months
## [1] 144  73 826
#plot a 90S-90N precip. along a meridional line at 160E over Pacific 
par(mfrow=c(1,1)) 
plot(seq(90,-90,length=73),precnc[65,,1],
     type="l", xlab="Latitude", 
     ylab="Temperature [°C]",
     main="90S-90N Temperature [°C] Along a Meridional Line at 160E: Jan. 1948",
     lwd=3)
grid(nx = NULL, ny = NULL)

 

Figure 9.9 (a)

#Compute and plot climatology and standard deviation Jan 1948-Dec 2015 
library(maps)
climmat <- matrix(0, nrow=144, ncol=73)
sdmat <- matrix(0, nrow=144, ncol=73)

Jmon <- 12*seq(0,67,1) 
for (i in 1:144){
  for (j in 1:73) {
    climmat[i,j] <- mean(precnc[i,j,Jmon])
    sdmat[i,j] <- sd(precnc[i,j,])
  }
}
mapmat <- climmat
#Note that R requires coordinates increasing from south to north -90->90 
#and from west to east from 0->360. We must arrange Lat and Lon this way. 

#Correspondingly, we have to flip the data matrix left to right according to
#the data matrix precnc[i,j,]: 360 (i.e. 180W) lon and from North Pole 
#and South Pole, then lon 178.75W, 176.75W, ..., 0E. This puts Greenwich 
#at the center, China on the right, and USA on the left. However, our map should
#have the Pacific at the center, and USA on the right. Thus, we make a flip .

Lat <- -Lat1

mapmat <- mapmat[,length(mapmat[1,]):1]
#Matrix flip around a column 
#mapmat<- t(apply(t(mapmat),2,rev))
int <- seq(-50,50,length.out=81) 
rgb.palette <- colorRampPalette(c('black','blue','darkgreen',
          'green','white', 'yellow','pink','red','maroon'),
                             interpolate='spline') 

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

 

Figure 9.9 (b)

#Plot standard deviation
#plot.new()
par(mgp=c(2,1,0))
par(mar=c(3,3,2,2))
mapmat <- sdmat[,seq(length(sdmat[1,]),1)]
int <- seq(0,20,length.out=81)
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 1948-2015 Jan. SAT RA Standard Deviation [°C]",
                    xlab="Longitude", ylab="Latitude"), 
        plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()}, 
               key.title=title(main="[°C]"))