################################################################
#
# 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("/rCode")
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()
## png
## 2
###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()
## png
## 2
#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()
## png
## 2
#
#
#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()
## png
## 2
#
#
#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] 6.518766e-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,] 0.9817996
#1.010333 #Approximately normalized
t(pcm[,1])%*%pcm[,2]
## [,1]
## [1,] 0.009436873
#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
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.101, p-value = 0.9414
## 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()
## png
## 2
#
#
#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()
## png
## 2
#
#
#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] -5.576865e-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()
## png
## 2
#
#
#
#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()
## png
## 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 -9.788388e-19 9.513623e-18 4.298203e-18 -2.084895e-18
## [2,] -9.788388e-19 9.801000e-01 -1.681911e-18 2.067985e-19 -1.396056e-18
## [3,] 9.513623e-18 -1.681911e-18 9.801000e-01 -7.331500e-18 7.596827e-20
## [4,] 4.298203e-18 2.067985e-19 -7.331500e-18 9.801000e-01 -8.041466e-18
## [5,] -2.084895e-18 -1.396056e-18 7.596827e-20 -8.041466e-18 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()
## png
## 2