Statistics and Data Visualizations in Climate Science

with R and Python

 

A Cambridge University Press Book by

SSP Shen and GR North

 

 

Version 1.0 released in July 2023 and coded by Dr. Samuel Shen, Distinguished Professor

San Diego State University, California, USA

https://shen.sdsu.edu

Email: sshen@sdsu.edu

 

 

 

 

Chapter 6: Covariance Matrices, EOFs, and PCs

Plot Fig. 6.1a, b, c: and Covariance Map by R code

#Covariance from NOAAGlaobalTemp data: December SAT anomalies  
#setwd('/Users/sshen/climstats') 
#install.packages('ncdf4')
library(ncdf4)
nc=ncdf4::nc_open("data/air.mon.anom.nc") #Read data
nc #Check metadata of the dataset

## File data/air.mon.anom.nc (NC_FORMAT_NETCDF4_CLASSIC):
## 
##      2 variables (excluding dimension variables):
##         double time_bnds[nbnds,time]   (Chunking: [2,1])  (Compression: shuffle,level 2)
##             long_name: Time Boundaries
##         float air[lon,lat,time]   (Chunking: [72,36,1])  (Compression: shuffle,level 2)
##             var_desc: Air Temperature
##             level_desc: Surface
##             statistic: Anomaly
##             parent_stat: Observation
##             valid_range: -40
##              valid_range: 40
##             units: degC
##             missing_value: -9.96920996838687e+36
##             long_name: Surface Air Temperature and SST Monthly Anomaly
##             precision: 2
##             cell_methods: time: anomaly (monthly from values)
##             standard_name: air_temperature_anomaly
##             dataset: NOAA Global Temperature
##             actual_range: -20.1061992645264
##              actual_range: 20.0300006866455
##             date_of_file_acquired: 2019-8-2
## 
##      4 dimensions:
##         time  Size:1674   *** is unlimited *** 
##             units: days since 1800-1-1 00:00:0.0
##             long_name: Time
##             delta_t: 0000-01-00 00:00:00
##             avg_period: 0000-01-00 00:00:00
##             standard_name: time
##             axis: T
##             coordinate_defines: start
##             bounds: time_bnds
##             calendar: gregorian
##             coverage_content_type: coordinate
##             actual_range: 29219
##              actual_range: 80139
##         lat  Size:36 
##             long_name: Latitude
##             units: degrees_north
##             actual_range: -87.5
##              actual_range: 87.5
##             axis: Y
##             standard_name: latitude
##             grids: Uniform grid from -87.5 to 87.5 by 5
##             coordinate_defines: center
##             _CoordinateAxisType: Lat
##         lon  Size:72 
##             long_name: Longitude
##             units: degrees_east
##             actual_range: 2.5
##              actual_range: 357.5
##             axis: X
##             standard_name: longitude
##             grids: Uniform grid from 2.5 to 357.5 by 5
##             coordinate_defines: center
##             _CoordinateAxisType: Lon
##         nbnds  Size:2 (no dimvar)
## 
##     22 global attributes:
##         Conventions: CF-1.0
##         Source: ftp://ftp.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/
##         dataset_title: NOAA Global Surface Temperature (NOAAGlobalTemp)
##         References: https://www.esrl.noaa.gov/psd/data/gridded/data.noaaglobaltemp.html
##         keywords_vocabulary: Climate and Forecast (CF) Standard Name Table (Version 46, 25 July 2017)
##         keywords: Earth Science > Oceans > Ocean Temperature > Sea Surface Temperature, Earth Science > Atmosphere > Atmospheric Temperature > Surface Temperature > Air Temperature
##         cdm_data_type: Grid
##         dataset_citation_url: https://doi.org/10.25921/9qth-2p70
##         references: Vose, R. S., et al., 2012: NOAAs merged land-ocean surface temperature analysis. Bulletin of the American Meteorological Society, 93, 1677-1685. doi: 10.1175/BAMS-D-11-00241.1. Huang, B., Peter W. Thorne, et. al, 2017: Extended Reconstructed Sea Surface Temperature version 5 (ERSSTv5), Upgrades, validations, and intercomparisons. J. Climate, 30, 8179-8205. doi: 10.1175/JCLI-D-16-0836.1
##         climatology: Climatology is based on 1971-2000 monthly climatology
##         license: These data are available for use without restriction.
##         source: https://www.ncdc.noaa.gov/noaa-merged-land-ocean-global-surface-temperature-analysis-noaaglobaltemp-v5
##         platform: Ships, moored buoys, surface drifting buoys, Argo floats, and weather stations
##         instrument: Conventional thermometers
##         creator_email: Boyin.Huang@noaa.gov, Xungang.Yin@noaa.gov
##         comment: Merged land ocean surface temperature anomalies. Version 5.0.0, blending ERSST V5 and GHCN-M V4
##         source_institution: DOC/NOAA/NESDIS/National Centers for Environmental Information(NCEI)
##         history: Created 08/2019 using new V5  data from NCEI
##         version: V5
##         geospatial_bounds: POLYGON ((2.5 -87.5, 2.5 87.5, 357.5 87.5, 357.5 -87.5, 2.5 -87.5))
##         title: NOAA Merged Land Ocean Global Surface Temperature Analysis (NOAAGlobalTemp)
##         data_modified: 2019-08-02

nc$dim$lon$vals # 2.5 - 357.5

##  [1]   2.5   7.5  12.5  17.5  22.5  27.5  32.5  37.5  42.5  47.5  52.5  57.5
## [13]  62.5  67.5  72.5  77.5  82.5  87.5  92.5  97.5 102.5 107.5 112.5 117.5
## [25] 122.5 127.5 132.5 137.5 142.5 147.5 152.5 157.5 162.5 167.5 172.5 177.5
## [37] 182.5 187.5 192.5 197.5 202.5 207.5 212.5 217.5 222.5 227.5 232.5 237.5
## [49] 242.5 247.5 252.5 257.5 262.5 267.5 272.5 277.5 282.5 287.5 292.5 297.5
## [61] 302.5 307.5 312.5 317.5 322.5 327.5 332.5 337.5 342.5 347.5 352.5 357.5

nc$dim$lat$vals # -87.5 - 87.5

##  [1] -87.5 -82.5 -77.5 -72.5 -67.5 -62.5 -57.5 -52.5 -47.5 -42.5 -37.5 -32.5
## [13] -27.5 -22.5 -17.5 -12.5  -7.5  -2.5   2.5   7.5  12.5  17.5  22.5  27.5
## [25]  32.5  37.5  42.5  47.5  52.5  57.5  62.5  67.5  72.5  77.5  82.5  87.5

nc$dim$time$vals

##    [1] 29219 29250 29279 29310 29340 29371 29401 29432 29463 29493 29524 29554
##   [13] 29585 29616 29644 29675 29705 29736 29766 29797 29828 29858 29889 29919
##   [25] 29950 29981 30009 30040 30070 30101 30131 30162 30193 30223 30254 30284
##   [37] 30315 30346 30374 30405 30435 30466 30496 30527 30558 30588 30619 30649
##   [49] 30680 30711 30740 30771 30801 30832 30862 30893 30924 30954 30985 31015
##   [61] 31046 31077 31105 31136 31166 31197 31227 31258 31289 31319 31350 31380
##   [73] 31411 31442 31470 31501 31531 31562 31592 31623 31654 31684 31715 31745
##   [85] 31776 31807 31835 31866 31896 31927 31957 31988 32019 32049 32080 32110
##   [97] 32141 32172 32201 32232 32262 32293 32323 32354 32385 32415 32446 32476
##  [109] 32507 32538 32566 32597 32627 32658 32688 32719 32750 32780 32811 32841
##  [121] 32872 32903 32931 32962 32992 33023 33053 33084 33115 33145 33176 33206
##  [133] 33237 33268 33296 33327 33357 33388 33418 33449 33480 33510 33541 33571
##  [145] 33602 33633 33662 33693 33723 33754 33784 33815 33846 33876 33907 33937
##  [157] 33968 33999 34027 34058 34088 34119 34149 34180 34211 34241 34272 34302
##  [169] 34333 34364 34392 34423 34453 34484 34514 34545 34576 34606 34637 34667
##  [181] 34698 34729 34757 34788 34818 34849 34879 34910 34941 34971 35002 35032
##  [193] 35063 35094 35123 35154 35184 35215 35245 35276 35307 35337 35368 35398
##  [205] 35429 35460 35488 35519 35549 35580 35610 35641 35672 35702 35733 35763
##  [217] 35794 35825 35853 35884 35914 35945 35975 36006 36037 36067 36098 36128
##  [229] 36159 36190 36218 36249 36279 36310 36340 36371 36402 36432 36463 36493
##  [241] 36524 36555 36583 36614 36644 36675 36705 36736 36767 36797 36828 36858
##  [253] 36889 36920 36948 36979 37009 37040 37070 37101 37132 37162 37193 37223
##  [265] 37254 37285 37313 37344 37374 37405 37435 37466 37497 37527 37558 37588
##  [277] 37619 37650 37678 37709 37739 37770 37800 37831 37862 37892 37923 37953
##  [289] 37984 38015 38044 38075 38105 38136 38166 38197 38228 38258 38289 38319
##  [301] 38350 38381 38409 38440 38470 38501 38531 38562 38593 38623 38654 38684
##  [313] 38715 38746 38774 38805 38835 38866 38896 38927 38958 38988 39019 39049
##  [325] 39080 39111 39139 39170 39200 39231 39261 39292 39323 39353 39384 39414
##  [337] 39445 39476 39505 39536 39566 39597 39627 39658 39689 39719 39750 39780
##  [349] 39811 39842 39870 39901 39931 39962 39992 40023 40054 40084 40115 40145
##  [361] 40176 40207 40235 40266 40296 40327 40357 40388 40419 40449 40480 40510
##  [373] 40541 40572 40600 40631 40661 40692 40722 40753 40784 40814 40845 40875
##  [385] 40906 40937 40966 40997 41027 41058 41088 41119 41150 41180 41211 41241
##  [397] 41272 41303 41331 41362 41392 41423 41453 41484 41515 41545 41576 41606
##  [409] 41637 41668 41696 41727 41757 41788 41818 41849 41880 41910 41941 41971
##  [421] 42002 42033 42061 42092 42122 42153 42183 42214 42245 42275 42306 42336
##  [433] 42367 42398 42427 42458 42488 42519 42549 42580 42611 42641 42672 42702
##  [445] 42733 42764 42792 42823 42853 42884 42914 42945 42976 43006 43037 43067
##  [457] 43098 43129 43157 43188 43218 43249 43279 43310 43341 43371 43402 43432
##  [469] 43463 43494 43522 43553 43583 43614 43644 43675 43706 43736 43767 43797
##  [481] 43828 43859 43888 43919 43949 43980 44010 44041 44072 44102 44133 44163
##  [493] 44194 44225 44253 44284 44314 44345 44375 44406 44437 44467 44498 44528
##  [505] 44559 44590 44618 44649 44679 44710 44740 44771 44802 44832 44863 44893
##  [517] 44924 44955 44983 45014 45044 45075 45105 45136 45167 45197 45228 45258
##  [529] 45289 45320 45349 45380 45410 45441 45471 45502 45533 45563 45594 45624
##  [541] 45655 45686 45714 45745 45775 45806 45836 45867 45898 45928 45959 45989
##  [553] 46020 46051 46079 46110 46140 46171 46201 46232 46263 46293 46324 46354
##  [565] 46385 46416 46444 46475 46505 46536 46566 46597 46628 46658 46689 46719
##  [577] 46750 46781 46810 46841 46871 46902 46932 46963 46994 47024 47055 47085
##  [589] 47116 47147 47175 47206 47236 47267 47297 47328 47359 47389 47420 47450
##  [601] 47481 47512 47540 47571 47601 47632 47662 47693 47724 47754 47785 47815
##  [613] 47846 47877 47905 47936 47966 47997 48027 48058 48089 48119 48150 48180
##  [625] 48211 48242 48271 48302 48332 48363 48393 48424 48455 48485 48516 48546
##  [637] 48577 48608 48636 48667 48697 48728 48758 48789 48820 48850 48881 48911
##  [649] 48942 48973 49001 49032 49062 49093 49123 49154 49185 49215 49246 49276
##  [661] 49307 49338 49366 49397 49427 49458 49488 49519 49550 49580 49611 49641
##  [673] 49672 49703 49732 49763 49793 49824 49854 49885 49916 49946 49977 50007
##  [685] 50038 50069 50097 50128 50158 50189 50219 50250 50281 50311 50342 50372
##  [697] 50403 50434 50462 50493 50523 50554 50584 50615 50646 50676 50707 50737
##  [709] 50768 50799 50827 50858 50888 50919 50949 50980 51011 51041 51072 51102
##  [721] 51133 51164 51193 51224 51254 51285 51315 51346 51377 51407 51438 51468
##  [733] 51499 51530 51558 51589 51619 51650 51680 51711 51742 51772 51803 51833
##  [745] 51864 51895 51923 51954 51984 52015 52045 52076 52107 52137 52168 52198
##  [757] 52229 52260 52288 52319 52349 52380 52410 52441 52472 52502 52533 52563
##  [769] 52594 52625 52654 52685 52715 52746 52776 52807 52838 52868 52899 52929
##  [781] 52960 52991 53019 53050 53080 53111 53141 53172 53203 53233 53264 53294
##  [793] 53325 53356 53384 53415 53445 53476 53506 53537 53568 53598 53629 53659
##  [805] 53690 53721 53749 53780 53810 53841 53871 53902 53933 53963 53994 54024
##  [817] 54055 54086 54115 54146 54176 54207 54237 54268 54299 54329 54360 54390
##  [829] 54421 54452 54480 54511 54541 54572 54602 54633 54664 54694 54725 54755
##  [841] 54786 54817 54845 54876 54906 54937 54967 54998 55029 55059 55090 55120
##  [853] 55151 55182 55210 55241 55271 55302 55332 55363 55394 55424 55455 55485
##  [865] 55516 55547 55576 55607 55637 55668 55698 55729 55760 55790 55821 55851
##  [877] 55882 55913 55941 55972 56002 56033 56063 56094 56125 56155 56186 56216
##  [889] 56247 56278 56306 56337 56367 56398 56428 56459 56490 56520 56551 56581
##  [901] 56612 56643 56671 56702 56732 56763 56793 56824 56855 56885 56916 56946
##  [913] 56977 57008 57037 57068 57098 57129 57159 57190 57221 57251 57282 57312
##  [925] 57343 57374 57402 57433 57463 57494 57524 57555 57586 57616 57647 57677
##  [937] 57708 57739 57767 57798 57828 57859 57889 57920 57951 57981 58012 58042
##  [949] 58073 58104 58132 58163 58193 58224 58254 58285 58316 58346 58377 58407
##  [961] 58438 58469 58498 58529 58559 58590 58620 58651 58682 58712 58743 58773
##  [973] 58804 58835 58863 58894 58924 58955 58985 59016 59047 59077 59108 59138
##  [985] 59169 59200 59228 59259 59289 59320 59350 59381 59412 59442 59473 59503
##  [997] 59534 59565 59593 59624 59654 59685 59715 59746 59777 59807 59838 59868
## [1009] 59899 59930 59959 59990 60020 60051 60081 60112 60143 60173 60204 60234
## [1021] 60265 60296 60324 60355 60385 60416 60446 60477 60508 60538 60569 60599
## [1033] 60630 60661 60689 60720 60750 60781 60811 60842 60873 60903 60934 60964
## [1045] 60995 61026 61054 61085 61115 61146 61176 61207 61238 61268 61299 61329
## [1057] 61360 61391 61420 61451 61481 61512 61542 61573 61604 61634 61665 61695
## [1069] 61726 61757 61785 61816 61846 61877 61907 61938 61969 61999 62030 62060
## [1081] 62091 62122 62150 62181 62211 62242 62272 62303 62334 62364 62395 62425
## [1093] 62456 62487 62515 62546 62576 62607 62637 62668 62699 62729 62760 62790
## [1105] 62821 62852 62881 62912 62942 62973 63003 63034 63065 63095 63126 63156
## [1117] 63187 63218 63246 63277 63307 63338 63368 63399 63430 63460 63491 63521
## [1129] 63552 63583 63611 63642 63672 63703 63733 63764 63795 63825 63856 63886
## [1141] 63917 63948 63976 64007 64037 64068 64098 64129 64160 64190 64221 64251
## [1153] 64282 64313 64342 64373 64403 64434 64464 64495 64526 64556 64587 64617
## [1165] 64648 64679 64707 64738 64768 64799 64829 64860 64891 64921 64952 64982
## [1177] 65013 65044 65072 65103 65133 65164 65194 65225 65256 65286 65317 65347
## [1189] 65378 65409 65437 65468 65498 65529 65559 65590 65621 65651 65682 65712
## [1201] 65743 65774 65803 65834 65864 65895 65925 65956 65987 66017 66048 66078
## [1213] 66109 66140 66168 66199 66229 66260 66290 66321 66352 66382 66413 66443
## [1225] 66474 66505 66533 66564 66594 66625 66655 66686 66717 66747 66778 66808
## [1237] 66839 66870 66898 66929 66959 66990 67020 67051 67082 67112 67143 67173
## [1249] 67204 67235 67264 67295 67325 67356 67386 67417 67448 67478 67509 67539
## [1261] 67570 67601 67629 67660 67690 67721 67751 67782 67813 67843 67874 67904
## [1273] 67935 67966 67994 68025 68055 68086 68116 68147 68178 68208 68239 68269
## [1285] 68300 68331 68359 68390 68420 68451 68481 68512 68543 68573 68604 68634
## [1297] 68665 68696 68725 68756 68786 68817 68847 68878 68909 68939 68970 69000
## [1309] 69031 69062 69090 69121 69151 69182 69212 69243 69274 69304 69335 69365
## [1321] 69396 69427 69455 69486 69516 69547 69577 69608 69639 69669 69700 69730
## [1333] 69761 69792 69820 69851 69881 69912 69942 69973 70004 70034 70065 70095
## [1345] 70126 70157 70186 70217 70247 70278 70308 70339 70370 70400 70431 70461
## [1357] 70492 70523 70551 70582 70612 70643 70673 70704 70735 70765 70796 70826
## [1369] 70857 70888 70916 70947 70977 71008 71038 71069 71100 71130 71161 71191
## [1381] 71222 71253 71281 71312 71342 71373 71403 71434 71465 71495 71526 71556
## [1393] 71587 71618 71647 71678 71708 71739 71769 71800 71831 71861 71892 71922
## [1405] 71953 71984 72012 72043 72073 72104 72134 72165 72196 72226 72257 72287
## [1417] 72318 72349 72377 72408 72438 72469 72499 72530 72561 72591 72622 72652
## [1429] 72683 72714 72742 72773 72803 72834 72864 72895 72926 72956 72987 73017
## [1441] 73048 73079 73108 73139 73169 73200 73230 73261 73292 73322 73353 73383
## [1453] 73414 73445 73473 73504 73534 73565 73595 73626 73657 73687 73718 73748
## [1465] 73779 73810 73838 73869 73899 73930 73960 73991 74022 74052 74083 74113
## [1477] 74144 74175 74203 74234 74264 74295 74325 74356 74387 74417 74448 74478
## [1489] 74509 74540 74569 74600 74630 74661 74691 74722 74753 74783 74814 74844
## [1501] 74875 74906 74934 74965 74995 75026 75056 75087 75118 75148 75179 75209
## [1513] 75240 75271 75299 75330 75360 75391 75421 75452 75483 75513 75544 75574
## [1525] 75605 75636 75664 75695 75725 75756 75786 75817 75848 75878 75909 75939
## [1537] 75970 76001 76030 76061 76091 76122 76152 76183 76214 76244 76275 76305
## [1549] 76336 76367 76395 76426 76456 76487 76517 76548 76579 76609 76640 76670
## [1561] 76701 76732 76760 76791 76821 76852 76882 76913 76944 76974 77005 77035
## [1573] 77066 77097 77125 77156 77186 77217 77247 77278 77309 77339 77370 77400
## [1585] 77431 77462 77491 77522 77552 77583 77613 77644 77675 77705 77736 77766
## [1597] 77797 77828 77856 77887 77917 77948 77978 78009 78040 78070 78101 78131
## [1609] 78162 78193 78221 78252 78282 78313 78343 78374 78405 78435 78466 78496
## [1621] 78527 78558 78586 78617 78647 78678 78708 78739 78770 78800 78831 78861
## [1633] 78892 78923 78952 78983 79013 79044 79074 79105 79136 79166 79197 79227
## [1645] 79258 79289 79317 79348 79378 79409 79439 79470 79501 79531 79562 79592
## [1657] 79623 79654 79682 79713 79743 79774 79804 79835 79866 79896 79927 79957
## [1669] 79988 80019 80047 80078 80108 80139

lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
time<- ncvar_get(nc, "time")
library(chron)
month.day.year(29219,c(month = 1, day = 1, year = 1800))

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

#1880-01-01
sat<- ncvar_get(nc, "air")
dim(sat)

## [1]   72   36 1674

#[1] 72   36 1674 
#1674 months=1880-01 to 2019-06, 139 years 6 mons

Dec = seq(12, 1674, by=12)
Decsat=sat[,, Dec]
N = 72*36
P = length(Dec)
STsat = matrix(0, nrow=N, ncol=P)
for (k in 1:P){STsat[,k]=as.vector(Decsat[,,k])}
colnames(STsat)<-1880:2018
STsat[1:4,1:4]

##      1880 1881 1882 1883
## [1,]   NA   NA   NA   NA
## [2,]   NA   NA   NA   NA
## [3,]   NA   NA   NA   NA
## [4,]   NA   NA   NA   NA

LAT=rep(lat, each=72)
LON=rep(lon,36)
STanom=cbind(LAT, LON, STsat)
dim(STanom)

## [1] 2592  141

#[1] 2592  141

#Plot Fig. 6.1a: Zonal covariance matrix 
#Select only the data for the equatorial band -2.5S
n1<-which(STanom[,1]>-4&STanom[,1]<0
          &STanom[,2]>0&STanom[,2]<360)
dim(STanom)

## [1] 2592  141

#[1] 2592  141
length(n1)

## [1] 72

#[1] 72 longitude grid boxes at -2.5S
P1=84
P2=133
P= P2 - P1 + 1
dat1=STanom[n1,P1:P2] #1961-2010
dim(dat1)

## [1] 72 50

#[1] 72 50
#72 grid boxes, 50 years of Dec from 1961-2010. 
dat1[1:3,48:50]

##           2008      2009      2010
## [1,] 0.3680965 0.3946191 0.2593701
## [2,] 0.3054363 0.6294945 0.3766301
## [3,] 0.9084032 0.4923795 0.5548739

Lat1=STanom[n1,1]
Lon1=STanom[n1,2]
AreaFac = sqrt(cos(Lat1*pi/180))
dat2 = dat1 - rowMeans(dat1) #Minus the mean
Banddat = AreaFac*dat2
covBand = Banddat%*%t(Banddat)
max(covBand)

## [1] 90.67208

#[1] 90.67199
min(covBand)

## [1] -6.181685

#[1] -6.191734
int=seq(-7, 92,length.out=81)
rgb.palette=colorRampPalette(c('blue', 'darkgreen',
 'green', 'yellow','pink','red','maroon'),
                             interpolate='spline')
#setEPS() #Plot the figure and save the file
#postscript("fig0601a.eps", width = 7, height = 5.5)
par(mar=c(4.2,5.0,1.8,0.0))
par(cex.axis=1.8,cex.lab=1.8, cex.main=1.7)
ticks = c(60,  180, 300)
filled.contour(Lon1, Lon1, covBand, 
               color.palette=rgb.palette, levels=int,
               plot.title=title(main=
 expression("Zonal SAT Covariance: 2.5"*degree*S),
    xlab="Longitude", ylab="Longitude"),
   plot.axes = { axis(1, at= ticks); axis(2, at= ticks, las = 0)},
            key.title=title(main=expression("["*degree*"C]"^2))
      )

#dev.off()

###variance as the diagonal covariance
varzonal = diag(covBand)/(AreaFac)^2
plot(Lon1, varzonal, type = 'h', lwd =2)

#Plot Fig. 6.1b: Meridional covariance matrix
#Select only the data for the meridional band 2.5E
n1<-which(STanom[,1]>-90&STanom[,1]<90
          &STanom[,2]>0&STanom[,2]<4)
P1=84 #time index for Dec 1961
P2=133 #Dec 2010
P= P2 - P1 + 1 #=50 Decembers
dat1=STanom[n1,P1:P2] #1961-2010
Lat1=STanom[n1,1]
Lon1=STanom[n1,2]
AreaFac = sqrt(cos(Lat1*pi/180)) 
dat2 = dat1 - rowMeans(dat1)
Banddat = AreaFac*dat2
covBand = Banddat%*%t(Banddat)
max(covBand, na.rm=TRUE)

## [1] 116.8833

#[1] 115.6266
min(covBand, na.rm=TRUE)

## [1] -34.57829

#[1] -33.72083
int=seq(-35,120,length.out=81)
rgb.palette=colorRampPalette(c('blue', 'darkgreen','green', 
                               'yellow','pink','red','maroon'),
                             interpolate='spline')
ticks = seq(-90, 90, by=30)
#setEPS() #Plot the figure and save the file
#postscript("fig0601b.eps", width = 7, height = 5.5)
par(mar=c(4.2,5.0,1.8,0.0))
par(cex.axis=1.8, cex.lab=1.8, cex.main=1.7)
filled.contour(Lat1, Lat1, covBand, 
               color.palette=rgb.palette, levels=int,
               plot.title=title(main=
      expression("SAT Meridional Covariance: 2.5"*degree*E),
                    xlab="Latitude", ylab="Latitude"),
plot.axes={axis(1, at=ticks); axis(2, at=ticks, las =0); 
           grid()},
      key.title=title(main=expression("["*degree*"C]"^2)))     

#dev.off()

#plot Fig. 6.1(c)
month.day.year(time[1416],c(month = 1, day = 1, year = 1800))

## $month
## [1] 12
## 
## $day
## [1] 1
## 
## $year
## [1] 1997

#Dec 1997
mapmat= sat[,,1416]
mapmat=pmax(pmin(mapmat,5),-5)
int=seq(-5, 5,length.out=51)
rgb.palette=colorRampPalette(c('black','blue',
  'darkgreen','green', 'white','yellow','pink',
  'red','maroon'), interpolate='spline')
#setEPS() #Plot the figure and save the file
#postscript("fig0601c.eps", width = 7, height = 3.5)
par(mar=c(4.2,5.0,1.8,0.0))
par(cex.axis=0.9,cex.lab=0.9, cex.main=0.8)
library(maps)
filled.contour(lon, lat, mapmat, 
               color.palette=rgb.palette, levels=int,
          plot.title=title(main="December 1997 SAT Anomalies",
                          xlab="Longitude",ylab="Latitude"),
               plot.axes={axis(1); axis(2); 
                 map('world2', add=TRUE);grid()},
               key.title=title(main=expression(degree*"C")))
segments(x0=-20,y0=-2.5, x1=255, y1=-2.5, lwd = 5, lty = 2)
segments(x0=-15,y0=-90, x1=-15, y1=90, lwd = 5, lty = 2)

#dev.off()

#
#
#Zonal covariance matrix 
#Select only the data for the equatorial band -2.5S
n1<-which(STanom[,1]>-4&STanom[,1]<0
          &STanom[,2]>0&STanom[,2]<360)
dim(STanom)

## [1] 2592  141

#[1] 2592  142
length(n1)

## [1] 72

#[1] 72 longitude grid boxes at -2.5S
P1=84
P2=133
P= P2 - P1 + 1
dat1=STanom[n1,P1:P2] #1961-2010
dim(dat1)

## [1] 72 50

#[1] 72 50
#72 grid boxes, 50 years of Dec from 1961-2010. 
dat1[1:3,47:50]

##           2007      2008      2009      2010
## [1,] 0.3385578 0.3680965 0.3946191 0.2593701
## [2,] 0.1402067 0.3054363 0.6294945 0.3766301
## [3,] 0.2994784 0.9084032 0.4923795 0.5548739

Lat1=STanom[n1,1]
Lon1=STanom[n1,2]
AreaFac = sqrt(cos(Lat1*pi/180))
dat2 = dat1 - rowMeans(dat1) #Minus the mean
Banddat = AreaFac*dat2
covBand = Banddat%*%t(Banddat)

 

R plot 6.2: Scree plot

K = 10
eigCov =eigen(covBand) 
#covBand is for equatorial zonal band in Fig 6.1(a)
lam = eigCov$values
lamK=lam[1:K]
#setEPS() #Plot the figure and save the file
#postscript("fig0602.eps", width = 6, height = 4)
par(mar=c(4,4,2,4), mgp=c(2.2,0.7,0))
plot(1:K, 100*lamK/sum(lam), ylim=c(0,80), type="o", 
     ylab="Percentage of Variance [%]",
     xlab="EOF Mode Number", 
     cex.lab=1.2, cex.axis = 1.1, lwd=2, 
     main="Scree Plot of the First 10 Eigenvalues")
legend(3,30, col=c("black"),lty=1, lwd=2.0,
       legend=c("Percentange Variance"),bty="n",
       text.font=2,cex=1.0, text.col="black")
par(new=TRUE)
plot(1:K,cumsum(100*lamK/sum(lam)),
     ylim = c(60,100), type="o",
     col="blue",lwd=2, axes=FALSE,
     xlab="",ylab="")
legend(3,80, col=c("blue"),lty=1,lwd=2.0,
       legend=c("Cumulative Percentage Variance"),bty="n",
       text.font=2,cex=1.0, text.col="blue")
axis(4, col="blue", col.axis="blue", mgp=c(3,0.7,0))
mtext("Cumulative Variance [%]",col="blue", 
      cex=1.2, side=4,line=2)

#dev.off()

 

Generate a Random Space-Time Field

#Step 1: Generate EOFs
N <- 100 #Number of spatial points
eof  <-  function(n, x) (sin(n*x)/sqrt(pi/2))*sqrt((pi/N))
x  <-  seq(0,pi, len=N)
sum(eof(4,x)^2)

## [1] 0.99

#[1] 0.99  #Verify the normalization condition
sum(eof(1,x)*eof(2,x))

## [1] 2.615638e-18

#[1] 3.035766e-18  #Verify orthogonality 

#Step 2: Generate PCs 
#install.packages('mvtnorm')
library(mvtnorm) #Multivariate normal
Mode <- 5
Ms <- 1000
univar <- diag(rep(1,Mode))
zeromean <- rep(0, Mode)
rs <- rmvnorm(Ms, zeromean, univar)
pcm=matrix(0, nrow=Ms, ncol=Mode)
for(m in 1:Mode){pcm[,m] = rs[,m]*sqrt(1/Ms)}
t(pcm[,1])%*%pcm[,1]

##          [,1]
## [1,] 1.061961

#1.010333 #Approximately normalized
t(pcm[,1])%*%pcm[,2]

##             [,1]
## [1,] 0.003702674

#0.008040772 #Approximately independent/orthorgonal

#Step 3: Generate an independent random field
lam=c(10, 9, 4.1, 4, 2)
sqrlam = diag(sqrt(lam))
eofm = matrix(0, nrow=N, ncol=Mode)
for(m in 1:Mode){eofm[,m]=eof(m,x)}
Yxr  <-  eofm%*%sqrlam%*%t(pcm)
dim(Yxr)

## [1]  100 1000

#[1] 100 1000

 

Durbin-Watson (DW) Test for No Serial Correlation

#install.packages('lmtest')
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Ms = 1000
r = 1:Ms
regYxr = lm(Yxr[10,] ~ r)
dwtest(regYxr)

## 
##  Durbin-Watson test
## 
## data:  regYxr
## DW = 2.08, p-value = 0.8916
## alternative hypothesis: true autocorrelation is greater than 0

#DW = 2.0344, p-value = 0.696
#Implying no significant serial correlation

 

Plot Fig. 6.3: R code

#The first 100 samples of the generated field
r=1:100
#setEPS()
#postscript("fig0603.eps", width = 6, height = 4)
par(mar=c(3.5,3.5,2,0), mgp=c(2.0,0.7,0))
rgb.palette=colorRampPalette(
  c('black','blue','green', 
    'yellow','pink','red','maroon'),
  interpolate='spline')
filled.contour(r,x, t(Yxr[,1:100]), 
         color=rgb.palette, xlab="r", ylab="x", cex.lab=1.2,
          main="A random field realization from given EOFs",
          plot.axes={axis(1, cex.axis =1.1); 
            axis(2, las = 0, cex.axis= 1.2); grid()},
          key.title={par(cex.main=0.9); title(main="Value")}
)

#dev.off()

 

Scree Plot

Mode = 1:5
lam=c(10, 9, 4.1, 4, 2)
samp = rep("S100", 5)
sd = sqrt(2/100)*lam
sd2 = sqrt(2/1000)*lam
par(mar=c(4.5, 4.5, 1, 0.5))
plot(Mode, lam, type="l", ylim = c(0,12),
     xlab="Mode", 
     ylab=expression("Variance   "*lambda),
     cex.lab=1.4, cex.axis=1.4)
points(Mode,lam + sd, pch="-", col="red", cex=2)
points(Mode,lam - sd, pch="-",col="red", cex=2)
segments(Mode,lam + sd, Mode,lam - sd, col="red")
points(Mode+0.06,lam + sd2, pch="-", col="blue", cex=2)
points(Mode+0.06,lam - sd2, pch="-",col="blue", cex=2)
segments(Mode +0.06,lam + sd2, 
         Mode + 0.06,lam - sd2, col="blue")
text(3.5,12, 
     "Red standard error bar: 100 samples", col="red")
text(3.5,11, 
     "Blue standard error bar: 1000 samples", col="blue")

 

Plot Fig. 6.4: Scree Plot by R code

Mode = 1:5
lam=c(10.0, 9.0, 4.1, 4.0, 2.0)
samp = rep("S100", 5)
sd = sqrt(2/100)*lam
sd2 = sqrt(2/1000)*lam
#par(mar=c(4.5, 4.5, 1, 0.5))
#setEPS()
#postscript("fig0604.eps", width = 10, height = 8)
par(mar=c(3.5, 4.5, 0.2, 0.5), mgp=c(2.5, 1.0,0))
plot(Mode, lam, type="l", ylim = c(0,12),
     xlab="Mode", 
     ylab=expression("Variance   "*lambda),
     cex.lab=1.8, cex.axis=1.8)
points(Mode,lam + sd, pch="-", col="red", cex=2)
points(Mode,lam - sd, pch="-",col="red", cex=2)
segments(Mode,lam + sd, Mode,lam - sd, col="red")
points(Mode+0.06,lam + sd2, pch="-", col="blue", cex=2)
points(Mode+0.06,lam - sd2, pch="-",col="blue", cex=2)
segments(Mode +0.06,lam + sd2, Mode + 0.06,lam - sd2, 
         col="blue")
text(3.45,12, "Red standard error bar: 100 samples", 
         cex = 1.5, col="red")
text(3.5,11, "Blue standard error bar: 1000 samples", 
     cex = 1.5, col="blue")

#dev.off()

 

Plot Fig. 6.5: 1D EOF Errors by R code

#Ms= 100 samples for North's rule of thumb
#install.packages("mvtnorm")
library(mvtnorm)
set.seed(112)
M=100 #M samples or M independent time steps
N=100 #N spatial locations
Mode=5 # 5 EOF modes to be considered

#Generate PCs
lam=c(10.0, 9.0, 4.1, 4.0, 2.0)
round(sqrt(2/M)*lam, digits = 2)

## [1] 1.41 1.27 0.58 0.57 0.28

#[1] 1.41 1.27 0.58 0.57 0.28
univar =diag(rep(1,Mode)) #SD = 1
zeromean = rep(0, Mode) #mean = 0
rs <- rmvnorm(M, zeromean, univar) 
dim(rs) # 100 samples and 5 modes

## [1] 100   5

#[1] 100    5
mean(rs[,1])

## [1] -0.02943547

#[1] -0.02524037
var(rs[,1])

## [1] 0.9108917

#[1] 0.9108917
t <- seq(0, 2*pi, len=M)
a51 <-rs[,1]*sqrt(1/M)
sum(a51^2)

## [1] 0.9026492

#[1] 0.9026492 is the variance approximation

pcm=matrix(0, nrow=M, ncol=Mode)
for(m in 1:Mode){pcm[,m] = rs[,m]*sqrt(1/M)} 
dim(pcm) #random and independent PCs

## [1] 100   5

#[1] 100   5
sum(pcm[,1]^2) #verify the normality of a PC

## [1] 0.9026492

#[1] 1.021924

#Generate EOFs for spatial patterns
eof <- function(n, x) sqrt(2)* sin(n*x)*sqrt((1/N))
x <- seq(0, pi,len=N) #N locations within [0,1]
sum(eof(3,x)^2) #verify the normality of an EOF

## [1] 0.99

#[1] 0.99
sum(eof(1,x)*eof(2,x)) #verify the EOF orthogonality

## [1] -8.395791e-18

#[1] 3.035766e-18

eofm <- matrix(0, nrow=N, ncol=Mode)
for (m in 1:Mode){eofm[,m]=eof(m,x)}
dim(eofm) #eofm are the 5 modes of EOF data

## [1] 100   5

#[1] 100   5 #100 spatial locations and 5 modes

#Generate the random data with given EOFs
Lam = diag(sqrt(lam)) #eigenvalue matrix
da <- eofm%*%Lam%*%t(pcm) #spectral decomposition
dim(da) #random data at 100 spatial locations

## [1] 100 100

#[1]  100 100
svdda <- svd(da)
round((svdda$d[1:5])^2, digits =2)

## [1] 10.14  8.12  3.83  3.26  2.10

#[1] 10.14  8.12  3.83  3.26  2.10

png(file="fig0605a.png", width=600,height=300)
k=1
sum(svdda$u[,k]^2)

## [1] 1

#[1] 1
par(mar=c(4.5,4.7,2,0.2))
plot(x,-svdda$u[,k], type="l", ylim=c(-0.2,0.2),
     xlab="x", ylab=paste("EOF", k), 
     main="EOF Mode (Blue Curve) vs Exact Mode (Red Curve)",
     col="blue",
     cex.lab=1.4, cex.axis=1.4,cex.main=1.4)
lines(x, eof(k,x), col='red')
#dev.off()

 

Plot Fig. 6.6: R Code for the EOF4 Error

k = 4 #Use the SVD result from 100 samples R = 100
plot(x,svdda$u[,k], type="l", ylim=c(-0.33,0.33),
     xlab="x", ylab="EOFs [dimensionless]", 
     main="EOF4 error (black) vs Exact EOF3 (Orange)",
     col="blue",
     cex.lab=1.4, cex.axis=1.4, cex.main=1.4)
legend(0.2,0.37, bty = "n", cex=1.4, text.col = 'blue',
       lwd=1.2,legend="Sample EOF4",col="blue")
lines(x, eof(k,x), col='red')
legend(0.2,0.41, bty = "n", cex=1.4, text.col = 'red',
       lwd=1.2,legend="True EOF4",col="red")
lines(x,svdda$u[,k] - eof(k,x), lwd=2, col='black'  )
legend(0.2,0.33, bty = "n", cex=1.4, text.col = 'black',
       lwd=2,legend="Sample EOF4 - True EOF4",col="black")
lines(x, -eof(3,x), col='orange')
legend(0.2,0.29, bty = "n", cex=1.4, text.col = 'orange',
       lwd=1.2,legend="True EOF3",col="orange")

 

Plot Fig. 6.7: R Code for 2D Sample EOFs

library(mvtnorm)
dmvnorm(x=c(0,0))

## [1] 0.1591549

M=100
N=100
Mode=5

#Generate PCs
lam=c(10.0, 9.0, 4.1, 4.0, 2.0)
univar =diag(rep(1,Mode))
zeromean = rep(0, Mode)
rs <- rmvnorm(M, zeromean, univar)
dim(rs)

## [1] 100   5

#[1] 100    5
t <- seq(0, len=M)
a51 <-rs[,1]*sqrt(1/M)
pcm=matrix(0, nrow=M, ncol=Mode)
for(m in 1:Mode){pcm[,m] = rs[,m]*sqrt(1/M)} 
dim(pcm)

## [1] 100   5

#[1] 100   5

#Generate true EOFs
eof <- function(n, x, y){(outer(sin(n*x),sin(n*y)))*(2/N)}
#eof <- function(n, x, y){outer(sin(n*x),sin(n*y))}
x = y <- seq(0,pi,len=100)
sum((eof(2,x, y))^2) #verify the normality

## [1] 0.9801

#[1] 0.9801

#Plot true EOF1 2D
png(file="fig0607a.png", width=200,height=200)
par(mar=c(2,2,0.5,0.5))
contour(x,y, eof(1,x,y))
text(1.5,1.5, 'EOF1', cex=2)
dev.off()

## quartz_off_screen 
##                 2

eofm <- matrix(0, nrow=N^2, ncol=Mode)
for (m in 1:Mode){eofm[,m]=c(eof(n=m,x,y))}
dim(eofm)

## [1] 10000     5

#[1] 10000   5 #10000 = 100*100
t(eofm)%*%eofm 

##               [,1]          [,2]          [,3]          [,4]          [,5]
## [1,]  9.801000e-01  2.300419e-18  6.694571e-18 -3.226918e-18 -9.271371e-20
## [2,]  2.300419e-18  9.801000e-01 -6.181324e-18 -4.347857e-18 -2.043335e-19
## [3,]  6.694571e-18 -6.181324e-18  9.801000e-01 -1.471922e-17 -6.567735e-18
## [4,] -3.226918e-18 -4.347857e-18 -1.471922e-17  9.801000e-01 -1.107744e-17
## [5,] -9.271371e-20 -2.043335e-19 -6.567735e-18 -1.107744e-17  9.801000e-01

#approximately equal to I5 an identify matrix

#Generate the random data with given EOFs
Lam = diag(sqrt(lam))
da <- eofm%*%Lam%*%t(pcm)
dim(da)

## [1] 10000   100

#[1] 10000   100
svdda <- svd(da)

#Plot ssample EOF1-5 2D: R=100
#png(file="fig0607n.png", width=200,height=200)
par(mar=c(2,2,0.5,0.5))
k=5
contour(x,y, matrix(svdda$u[,k], ncol=100), 
        cex.lab=1.4, cex.main=1.4)
text(1.5,1.5, 'EOF5', cex=2)

#dev.off()