Precipitation Practice

This is a practice run of the code given to me by Sam Shen.

First, let’s install all necessary packages and set the working dirextory (make sure to change this to your local machine). If already installed, do not remove the #. The program will wite files to the working directory, so make sure you designate an appropritate location on your machine.

Now, we can read in the nc file, using the ncdf4 package, that contains the precipitation data. This can be found at https://www.esrl.noaa.gov/psd/data/gridded/data.gpcp.html. Let’s open this file up and look at some of its variables. It is important to understand the contents of the nc file and the variables that it stores, since we will be interacting with this data.

#change the following line to have the correct file
nc=ncdf4::nc_open("E:/Precip Practice/precip.mon.mean.nc")
nc

File E:/Precip Practice/precip.mon.mean.nc (NC_FORMAT_NETCDF4_CLASSIC):

     4 variables (excluding dimension variables):
        float time_bnds[nv,time]   
            comment: time bounds for each time value
            units: days since 1800-01-01 00:00:00 0:00
        float lat_bnds[nv,lat]   
            units: degrees_north
            comment: latitude values at the north and south bounds of each pixel.
        float lon_bnds[nv,lon]   
            units: degrees_east
            comment: longitude values at the west and east bounds of each pixel.
        float precip[lon,lat,time]   
            long_name: Average Monthly Rate of Precipitation
            valid_range: 0
             valid_range: 100
            units: mm/day
            add_offset: 0
            scale_factor: 1
            missing_value: -9.96920996838687e+36
            precision: 32767
            least_significant_digit: 2
            var_desc: Precipitation
            dataset: GPCP Version 2.3 Combined Precipitation Dataset
            level_desc: Surface
            statistic: Mean
            parent_stat: Mean
            actual_range: 0
             actual_range: 47.3274345397949

     4 dimensions:
        lat  Size:72
            units: degrees_north
            actual_range: 88.75
             actual_range: -88.75
            long_name: Latitude
            standard_name: latitude
            axis: Y
        lon  Size:144
            units: degrees_east
            long_name: Longitude
            actual_range: 1.25
             actual_range: 358.75
            standard_name: longitude
            axis: X
        time  Size:459   *** 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
            actual_range: 65378
             actual_range: 79317
        nv  Size:2

    18 global attributes:
        _NCProperties: version=1|netcdflibversion=4.4.1|hdf5libversion=1.8.16
        Conventions: CF-1.0
        curator: Dr. Jian-Jian Wang
ESSIC, University of Maryland College Park
College Park, MD  20742  USA
Phone: +1 301-405-4887
        citation: Adler, R.F., G.J. Huffman, A. Chang, R. Ferraro, P. Xie, J. Janowiak, B. 
Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. 
Arkin, 2003: The Version 2 Global Precipitation Climatology Project 
(GPCP) Monthly Precipitation Analysis (1979 - Present). J. Hydrometeor., 
4(6), 1147-1167.
        title: GPCP Version 2.3 Combined Precipitation Dataset (Final)
        platform: NOAA POES (Polar Orbiting Environmental Satellites)
        source_obs: CDR RSS SSMI/SSMIS Tbs over ocean 
CDR SSMI/SSMIS rainrates over land (Ferraro) 
Geo-IR (Xie) calibrated by SSMI/SSMIS rainrates for sampling 
TOVS/AIRS empirical precipitation estimates at higher latitudes 
(ocean and land) 
GPCC gauge analysis to bias correct satellite estimates over land and 
merge with satellite based on sampling 
OLR Precipitation Index (OPI) (Xie) used for period before 1988
        documentation: http://www.esrl.noaa.gov/psd/data/gridded/data.gpcp.html
        version: V2.3
        Acknowledgement: 

        contributor_name: Robert Adler    University of Maryland 
George Huffman  NASA Goddard Space Flight Center 
David Bolvin    NASA Goddard Space Flight Center/SSAI 
Eric Nelkin     NASA Goddard Space Flight Center/SSAI 
Udo Schneider   GPCC, Deutscher Wetterdienst 
Andreas Becker  GPCC, Deutscher Wetterdienst 
Long Chiu       George Mason University 
Mathew Sapiano  University of Maryland 
Pingping Xie    Climate Prediction Center, NWS, NOAA 
Ralph Ferraro   NESDIS, NOAA 
Jian-Jian Wang  University of Maryland 
Guojun Gu       University of Maryland
        history: Generated at NOAA/ESRL PSD Sep 9 2016 based on data from source 
and stored in single netCDF4 file
        References: http://www.esrl.noaa.gov/psd/data/gridded/data.gpcp.html
        dataset_title: Global Precipitation Climatology Project (GPCP) Monthly Analysis Product
        first_interim_month: September 2016
        description: https://data.nodc.noaa.gov/cgi-bin/iso?id=gov.noaa.ncdc:C00970
        source: https://www.ncei.noaa.gov/data/global-precipitation-climatology-project-gpcp-monthly/access/
        source_documentation: https://www.ncdc.noaa.gov/cdr/atmospheric/precipitation-gpcp-monthly

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  26.25  28.75  31.25  33.75  36.25  38.75  41.25  43.75
 [19]  46.25  48.75  51.25  53.75  56.25  58.75  61.25  63.75  66.25  68.75  71.25  73.75  76.25  78.75  81.25  83.75  86.25  88.75
 [37]  91.25  93.75  96.25  98.75 101.25 103.75 106.25 108.75 111.25 113.75 116.25 118.75 121.25 123.75 126.25 128.75 131.25 133.75
 [55] 136.25 138.75 141.25 143.75 146.25 148.75 151.25 153.75 156.25 158.75 161.25 163.75 166.25 168.75 171.25 173.75 176.25 178.75
 [73] 181.25 183.75 186.25 188.75 191.25 193.75 196.25 198.75 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 251.25 253.75 256.25 258.75 261.25 263.75 266.25 268.75
[109] 271.25 273.75 276.25 278.75 281.25 283.75 286.25 288.75 291.25 293.75 296.25 298.75 301.25 303.75 306.25 308.75 311.25 313.75
[127] 316.25 318.75 321.25 323.75 326.25 328.75 331.25 333.75 336.25 338.75 341.25 343.75 346.25 348.75 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 -63.75 -61.25 -58.75 -56.25 -53.75 -51.25 -48.75 -46.25
[19] -43.75 -41.25 -38.75 -36.25 -33.75 -31.25 -28.75 -26.25 -23.75 -21.25 -18.75 -16.25 -13.75 -11.25  -8.75  -6.25  -3.75  -1.25
[37]   1.25   3.75   6.25   8.75  11.25  13.75  16.25  18.75  21.25  23.75  26.25  28.75  31.25  33.75  36.25  38.75  41.25  43.75
[55]  46.25  48.75  51.25  53.75  56.25  58.75  61.25  63.75  66.25  68.75  71.25  73.75  76.25  78.75  81.25  83.75  86.25  88.75

nc$dim$time$vals

  [1] 65378 65409 65437 65468 65498 65529 65559 65590 65621 65651 65682 65712 65743 65774 65803 65834 65864 65895 65925 65956 65987
 [22] 66017 66048 66078 66109 66140 66168 66199 66229 66260 66290 66321 66352 66382 66413 66443 66474 66505 66533 66564 66594 66625
 [43] 66655 66686 66717 66747 66778 66808 66839 66870 66898 66929 66959 66990 67020 67051 67082 67112 67143 67173 67204 67235 67264
 [64] 67295 67325 67356 67386 67417 67448 67478 67509 67539 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 68300 68331 68359 68390 68420 68451 68481 68512 68543
[106] 68573 68604 68634 68665 68696 68725 68756 68786 68817 68847 68878 68909 68939 68970 69000 69031 69062 69090 69121 69151 69182
[127] 69212 69243 69274 69304 69335 69365 69396 69427 69455 69486 69516 69547 69577 69608 69639 69669 69700 69730 69761 69792 69820
[148] 69851 69881 69912 69942 69973 70004 70034 70065 70095 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 70857 70888 70916 70947 70977 71008 71038 71069 71100
[190] 71130 71161 71191 71222 71253 71281 71312 71342 71373 71403 71434 71465 71495 71526 71556 71587 71618 71647 71678 71708 71739
[211] 71769 71800 71831 71861 71892 71922 71953 71984 72012 72043 72073 72104 72134 72165 72196 72226 72257 72287 72318 72349 72377
[232] 72408 72438 72469 72499 72530 72561 72591 72622 72652 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 73414 73445 73473 73504 73534 73565 73595 73626 73657
[274] 73687 73718 73748 73779 73810 73838 73869 73899 73930 73960 73991 74022 74052 74083 74113 74144 74175 74203 74234 74264 74295
[295] 74325 74356 74387 74417 74448 74478 74509 74540 74569 74600 74630 74661 74691 74722 74753 74783 74814 74844 74875 74906 74934
[316] 74965 74995 75026 75056 75087 75118 75148 75179 75209 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 75970 76001 76030 76061 76091 76122 76152 76183 76214
[358] 76244 76275 76305 76336 76367 76395 76426 76456 76487 76517 76548 76579 76609 76640 76670 76701 76732 76760 76791 76821 76852
[379] 76882 76913 76944 76974 77005 77035 77066 77097 77125 77156 77186 77217 77247 77278 77309 77339 77370 77400 77431 77462 77491
[400] 77522 77552 77583 77613 77644 77675 77705 77736 77766 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 78527 78558 78586 78617 78647 78678 78708 78739 78770
[442] 78800 78831 78861 78892 78923 78952 78983 79013 79044 79074 79105 79136 79166 79197 79227 79258 79289 79317

nc$dim$time$units

[1] "days since 1800-1-1 00:00:0.0"

Now, let us look at the output of the three variables: lon, lat and time.

• lon: or longitude, holds the degrees of longitude 1.25 to 358.75, by an increment of 2.5 degrees.

• lat: or latitude. holds the degrees of latitude -88.75 to 88.75, by an increment of 2.5 degrees.

• time: contains days since 1-1-1800.

Since we want to interpret and analyse this data, storing the data in local variables will be essential. In this next code chunk we will store the lon, lat, precip and time values into variables. As well, the chron package, which stands for chronological, will be used for the month.day.year() function in order to estalish the first and last days of our data. If you haven’t done so already make sure to have the chron and ncdf4 packages installed. The ncdf4 library will be necessary to store the data, using the ncvar_get() function.

library(ncdf4)
LON=ncvar_get(nc, "lon") 
LAT=ncvar_get(nc, "lat") 
Time1<- ncvar_get(nc, "time")
precipNC<- ncvar_get(nc, "precip")
head(Time1) #find the first 6 IDs

[1] 65378 65409 65437 65468 65498 65529

#[1] 65378 65409 65437 65468 65498 65529
tail(Time1) #find the last 6 IDs

[1] 79166 79197 79227 79258 79289 79317

#[1] 79166 79197 79227 79258 79289 79317
library(chron)
#The firt datum time ID 65378
month.day.year(65378,c(month = 1, day = 1, year = 1800))#1979-01-01

$month
[1] 1

$day
[1] 1

$year
[1] 1979

#The last datum time ID 79317
month.day.year(79317,c(month = 1, day = 1, year = 1800))#2017-03-01

$month
[1] 3

$day
[1] 1

$year
[1] 2017

You may be wondering what the percnc variable holds, since we have yet to discuss it. I thought that it was all to important to include with the other varaibles since it is essential to this project. Let’s recall what the nc file described the ‘precip’ variable as:

        float precip[lon,lat,time]   
        long_name: Average Monthly Rate of Precipitation
        valid_range: 0
         valid_range: 100
        units: mm/day
        add_offset: 0
        scale_factor: 1
        missing_value: -9.96920996838687e+36
        precision: 32767
        least_significant_digit: 2
        var_desc: Precipitation
        dataset: GPCP Version 2.3 Combined Precipitation Dataset
        level_desc: Surface
        statistic: Mean
        parent_stat: Mean
        actual_range: 0
         actual_range: 47.327434539794

From the description we see that precip is a three-dimensional array that has the structure: precip(lon, lat, time). The long name is the ‘Average Monthly Rate of Precipitation’ which is measured in mm/day.

The dimension of precipNC, the variable that now holds the 3D array:

dim(precipNC)

[1] 144  72 459

The output shows [1] 144 72 459. But what exactly does that mean?

• The 144 is the number of gridpoints for degrees longitude, from 1.25 to 358.75.

• The 72 is the number of gridpoints for degrees latitude, from-88.75 to 88.75.

• The 459 represents the the number of months in the 38 year span, from January 1979 to March 2017.

Next, we can start to look closer at this data by plotting a 90S-90N precip along a meridional line at 160E over Pacific. This next chunk will plot each latitude gridpoint against the precipNC values found at LAT[15] and January 1979.

plot(LAT,precipNC[15,,1], 
     type="l", xlab="Latitude", ylab="Precipiation [mm/day]",
     main="90S-90N precipitation [mm/day] 
     along a meridional line at 35E: Jan 1979",
     lwd=3)

Now, let’s take another look at this data by computing the climatology and standard deviation of time spam 1981 to 2010. This can be done with a simple nested for loop upon each longitude and latitude gridpoint, along with the mean() and sd() functions. This task is shown as to how the computation is done and a plot of the climatology is created.

climMAT = matrix(0,nrow=144,ncol=72) #stoarge matrix for the climatology 
sdMAT = matrix(0,nrow=144,ncol=72) #storage matrix for the standard deviations
mon=12*seq(0,30,1) + 25 #months to look over
for (i in 1:144){
  for (j in 1:72) {climMAT[i,j]=mean(precipNC[i,j,mon]); 
  sdMAT[i,j]=sd(precipNC[i,j,mon]) 
  }
}
library(maps)
int=seq(0,10,length.out=11)
rgb.palette=colorRampPalette(c('skyblue',  'green', 'blue', 'yellow', 'orange', 'pink','red', 
                               'maroon', 'purple', 'black'),interpolate='spline')
filled.contour(LON, LAT, climMAT, color.palette=rgb.palette, levels=int,
               plot.title=title(main="GPCP 1981-2010 Precipitation Climatotology",
                                xlab="Longitude", ylab="Latitude"),
               plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()},
               key.title=title(main="mm/day"))

This next chunk has the tasks of computing the EOFs from the space-time data and write our data as a space-time matrix. The three-dimensional precipNC array will be used in a for loop, where for each of the 459 months, all longitude and latitude gridpoints will be created as a vector for that month. This for loop will store each vector in a new matrix, precst. Next, a csv file will be written to the working directory that the user specifices. This csv will contain the data fir 75S to 75N from the precst space-time matrix that we created.

The space-time matrix has the structure:

The GPCP space-time data matrix from January 1979 to March 2017 Column1 is Grid ID, Column2 LAT, Column3, LON Row1 is the time 1979-1, 1979-2, …… LAT LON 1979 - 1 1979 - 2 1979 - 3 1979 - 4 1 -88.75 1.25 0 0.336575955 0.27318272 0.171142474 2 -88.75 3.75 0 0.316594929 0.279624701 0.176127926

precst=matrix(0,nrow=10368,ncol=459)
for (i in 1:459) {
  precst[,i]=as.vector(precipNC[ , , i])
}
LAT=rep(LAT, each=144)
LON=rep(LON,72)
yr1=c(rep(1979:2016,each=12), rep(2017,3))
mon1=c(rep(1:12,38),1:3)
timehead=paste(yr1,"-",mon1)
colnames(precst)=timehead
precst1=cbind(LAT,LON,precst)
precst2=precst1[865:9504,]
write.csv(precst2, file="75gpcpst79Jan17Mar.csv")

Now since we have the space-time matrix stored, we can compute the January to Decemeber climatologies and standard devaitions for 1981 to 2010. Thus, each month has a computaiton made! Since we are dealing with a space-time matrix, we can no longer use the mean() and sd() functions that we used in the previous chunk. We need to perform these operations on entire arrays! So, we will need to use the matrixStats library to call the rowMeans() and rowSds() functions, which find the means and standard deviations of entire rows.

library(matrixStats)
climMAT=matrix(0,nrow=8640, ncol=14)
climMAT[,1:2]=precst2[,1:2]
sdMAT=climMAT
for(i in 1:12){
  mon=seq(26+i,by=12, len=30)
  monDat=precst2[,mon]
  climMAT[,2+i]=rowMeans(monDat)
  sdMAT[,2+i]<-rowSds(monDat)
}

Since the climatologies and standard deviations have now been calculated, creating plots for each month will allow us to visualize the results. 12 seperate plots will be created foe each the climatology and the standard devaition, one for each month of the year from 1981 to 2010. The first chunk will plot the climatologies and the second will plot the standard deviations.

library(maps)
LAT=seq(-73.75,by=2.5,len=60)
LON=seq(1.25,by=2.5, len=144)
monID=format(ISOdate(2004,1:12,1),"%B")
for (m in 1:12){
mapMatrix=matrix(climMAT[,2+m],nrow=144)
mapMatrix=pmin(mapMatrix,10)
int=seq(0,10,length.out=11)
rgb.palette=colorRampPalette(c('skyblue',  'green', 'blue', 'yellow', 'orange', 'pink','red', 
                               'maroon', 'purple', 'black'),interpolate='spline')
filled.contour(LON, LAT, mapMatrix, color.palette=rgb.palette, levels=int,
               plot.title=title(main=paste(monID[m],"GPCP 1981-2010 Precipitation Climatotology"),
                                xlab="Longitude", ylab="Latitude"),
               plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()},
               key.title=title(main="mm/day"))
}

#Plot standard deviation 
library(maps)
LAT=seq(-73.75,by=2.5,len=60)
LON=seq(1.25,by=2.5, len=144)
monID=format(ISOdate(2004,1:12,1),"%B")
for (m in 1:12){
mapMatrix=matrix(sdMAT[,2+m],nrow=144)
mapMatrix=pmin(mapMatrix,5) #Compress the values >8 to 8
int=seq(0,5,length.out=11)
rgb.palette=colorRampPalette(c('skyblue',  'green', 'blue', 'yellow', 'orange', 'pink',
                               'red', 'maroon', 'purple', 'black'),interpolate='spline')
filled.contour(LON, LAT, mapMatrix, color.palette=rgb.palette, levels=int,
               plot.title=title(main=paste(monID[m], "GPCP 1981-2010 Precipitation Standard Deviation"),
                                xlab="Longitude", ylab="Latitude"),
               plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()},
               key.title=title(main="mm/day"))
}

Again, let’s do some computations. This time we can cimoute the non-standardized anomalies of our space-time matrix. The following chunk will create a new matrix tha computes the standardized anomaly matrix, anomST, and one that contains that area weighted anomaly matrix, anomW.

anomST=matrix(0, nrow=8640,ncol=461)
anomST[,1:2]=climMAT[,1:2]
colnames(anomST)<-colnames(precst2)
anomAW=anomST
colnames(anomAW)<-colnames(precst2)
for (m in 1:12){
  monC=seq(2+m,461, by=12)
  anomST[,monC]=precst2[,monC]-climMAT[,2+m]
}
anomAW[,3:461]=sqrt(cos(anomST[,1]*pi/180))*anomST[,3:461]

The following chunk will compute the EOF array for each month of the year. EOFs are found using the svd() function, when passing in the according space-time matrix, and using the U matrix, an mxm unitary matrix. Once we have these EOF arrays, we can compute the matricies of eigenvalues of non-standardized anomalies, the percentage of eigenvalues, cumulative eigenvalues of non-standardized anomalies, and the percentage of cumulative eigenvalues. We will write each of these to a csv file on the local machine.

EOFar=array(0, dim=c(12,8640,30)) #EOFar[month,gridID,modeNo]
PCar=array(0, dim=c(12,30,30))
EigenMat=matrix(0, nrow=30, ncol=12)
for (m in 1:12){
  monD=seq(26+m, by=12, len=30)
  svd8110=svd(anomAW[,monD]) #execute SVD
  EOFar[m,,]=svd8110$u
  colnames(EOFar[m,,])<-paste(rep("E",30),sep="",1:30)
  PCar[m,,]=svd8110$v
  colnames(PCar[m,,])<-paste(rep("PC",30),sep="",1:30)
  EigenMat[,m]=(svd8110$d)^2/30
}
colnames(EigenMat)<-monID
cEigenMat=pEigenMat=pcEigenMat=EigenMat
for(m in 1:12){
  cEigenMat[,m]=cumsum(EigenMat[,m])
  pEigenMat=100*t(t(EigenMat)/colSums(EigenMat))
  pcEigenMat=100*t(t(cEigenMat)/colSums(EigenMat))
}

Now we would like to look more closely at the EOFs of standardizd anomalies from 1981 to 2010. We computed this earlier, but can we write it to a csv file? No, the reason being is that it is a three-dimensional array, and a csv file cannot handle this. So, we will create a nc file!

filename="EOF_GPCP1981_2010.nc"
library(ncdf4)
xvals=seq(1.25,by=2.5,len=144)
yvals=seq(-73.75, by=2.5,len=60) 
nx=length(xvals)
ny=length(yvals)
nmode=30
nmonth=12
lon=ncdim_def("Longitude", "degrees increase", xvals)
lat=ncdim_def("Latitude", "degrees increase", yvals)
mode2=ncdim_def("Mode","numbers", 1:nmode)
time2=ncdim_def("Time","months", 1:nmonth)
var_eof=ncvar_def("EOFvalues", "dimensionless", 
                     list(lon, lat, mode2, time2), 
                     longname="EOF data from GPCP data 1981-2010: 
                               30modes 12months") 
ncnew=nc_create(filename, list(var_eof))
print(paste("The file has", ncnew$ndim,"dimensions"))

[1] "The file has 4 dimensions"

dat1=array(0,dim=c(144,60,30,12))
for (mo in 1:12){
  for (md in 1:30){
    dat1[,,md,mo]=matrix(EOFar[mo,,md],nrow=144)
  }
}
data=as.vector(dat1)
ncvar_put(ncnew, var_eof, data, start=c(1,1,1,1), count=c(nx,ny,nmode,nmonth))
nc_close(ncnew)

As we did in beginning, let’s read in this new nc file and take a look at the variables so we can better understand the data, and so we can check that we correctly created the nc file. We will also grab the EOF files and store in a local variable.

#Read the .nc file
nc=ncdf4::nc_open("E:/Precip Practice/EOF_GPCP1981_2010.nc")
dat3<- ncvar_get(nc, "EOFvalues")

Now, that we have this nc file read in, let’s plot the EOFs that we computed earlier. The following will plot the first EOFs.

LAT=seq(-73.75,by=2.5,len=60)
LON=seq(1.25,by=2.5, len=144)
monID=format(ISOdate(2004,1:12,1),"%B")
for (m in 1:12){
    n=1
    mapMatrix = - matrix(dat3[,,1,1]/sqrt(cos(LAT*pi/180)),nrow=144)
    mapMatrix = pmax(pmin(mapMatrix,0.03),-0.03)
    rgb.palette=colorRampPalette(c('red','yellow','orange', 'white', 
                                   'green','blue','darkblue'),interpolate='spline')
    int=seq(-0.03,0.03,length.out=61)
    mapMatrix=mapMatrix[, seq(length(mapMatrix[1,]),1)]
    filled.contour(LON, LAT, mapMatrix, color.palette=rgb.palette, levels=int,
                   plot.title=title(main=paste(monID[m],"GPCP Precipitation EOF",n),
                                    xlab="", ylab=""),
                   plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()},
                   key.title=title(main="Scale"))
    mtext("Longitude",side=1,line = 2,cex=1.3)
    mtext("Latitude",side=2,line = 2,cex=1.3)
}

Another set of information we have is the matricies that contain the percentage and cumulative of eigenvalues. Plotting those will also allow us to visualize the computed information.

#plot precentage eigenvalues and cumulatives
modeNo=1:30
for (m in 1:12) { 
   par(mar=c(4,5,2,4))
  plot(modeNo, pEigenMat[,i],type="o", ylim=c(0,20),
       col="red", lwd=2.0, xlab="",ylab="", xaxt="n", yaxt="n")
  title(ylab="Percentage Variance [%]",col.lab="red",
        main = paste("GPCP1981-2010 Covariance Eigenvalues for", monID[m], sep = " "),
        cex.lab=1.4)
  axis(1, cex.axis=1.4)
  axis(2, col.axis="red", cex.axis=1.4)
  legend(14,12, col=c("red"),lty=1,lwd=2.0,
         legend=c("Percentange Variance"),bty="n",
         text.font=2,cex=1.0, text.col="red")
  par(new=TRUE)
  plot(modeNo,pcEigenMat[,i], type="o",ylim=c(0,100),
       col="blue",lwd=2.0,axes=FALSE,xlab="",ylab="",xaxt="n", yaxt="n")
  axis(4, col.axis="blue", cex.axis=1.4)
  legend(14,70, col=c("blue"),lty=1,lwd=2.0,
         legend=c("Cumulative Percentage Variance [%]"),bty="n",
         text.font=2,cex=1.0, text.col="blue")
  mtext("Cumulative variance [%]",side=4,line=2.5, cex=1.4, col="blue")
  mtext("EOF Mode Number", side=1,line=2,cex=1.4)
}

Next, we can plot the 12 panels on the same figure one, which is the 12 plots found above, but

modeNo=1:30
par(mfrow = c(4, 3))  # 4 rows and 3 columns
for (m in 1:12) { 
  par(mar=c(4,4,2,4))
  plot(modeNo, pEigenMat[,m],type="o", ylim=c(0,20),
       xlab="", ylab="",
       main = paste("Eigenvalues for", monID[m], split = ""),
       col="red")
  mtext("EOF Mode Number",side=1, line = 2.0, cex=.8)
  mtext("PercentVar[%]",side=2, line = 2.0, cex=.6)
  legend(10,12, col=c("black"),lty=1,lwd=2.0,
         legend=c("PercentVar"),bty="n",
         text.font=1.0,cex=1.0, text.col="red")
  par(new=TRUE)
  plot(modeNo,pcEigenMat[,m], type="o",
       col="blue",lwd=1.5,axes=FALSE,xlab="",ylab="", ylim=c(0,100))
  legend(10,80, col=c("blue"),lty=1,lwd=2.0,
         legend=c("CumulativeVar"),bty="n",
         text.font=1.0,cex=1.0, text.col="blue")
  axis(4)
  mtext("Cum Var [%]",cex=0.6,side=4,line=2)
}

We can now plot the first six physocal EOFs for each month. The physical EOF is found by dividing each eigenvector by the area factor.

LAT=seq(-73.75,by=2.5,len=60)
LON=seq(1.25,by=2.5, len=144)
monID=format(ISOdate(2004,1:12,1),"%B")
for (m in 1:12) { 
  for (n in 1:6){
mapMatrix = - matrix(EOFar[m,,n]/sqrt(cos(climMAT[,1]*pi/180)),nrow=144)
mapMatrix = pmax(pmin(mapMatrix,0.03),-0.03)
rgb.palette=colorRampPalette(c('red','yellow','orange', 'white', 
                               'green','blue','darkblue'),interpolate='spline')
int=seq(-0.03,0.03,length.out=61)
mapMatrix=mapMatrix[, seq(length(mapMatrix[1,]),1)]
filled.contour(LON, LAT, mapMatrix, color.palette=rgb.palette, levels=int,
               plot.title=title(main=paste(monID[m],"GPCP Precipitation EOF",n),
                                xlab="", ylab=""),
               plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()},
               key.title=title(main="Scale"))
mtext("Longitude",side=1,line = 2,cex=1.3)
mtext("Latitude",side=2,line = 2,cex=1.3)
}
}