---

Python Code Version 3.0 for

Climate Mathematics: Theory and Applications - Chapter 10

---

Samuel S.P. Shen and Richard C.J. Somerville

---

This Version 3.0 is authored by Briana Ramirez, edited by Samuel Shen. Liu Yang, Sandra Villamar, and Joaquin Stawsky contributed codes to this version.

Video tutorial for the python code can be found at the following URL: https://www.youtube.com/channel/UC7D9i0kBMzPTHyEhU0h6W9g

This version is based upon the previous version described in the following box.

######################################################################################################################
#This Python Code for Climate Science is written for the book entitled "Climate Mathematics: Theory and Applications"#
#A Cambridge University Press book authored by SSP Shen and RCJ Somerville in July 2019                              #
#The Python codes were based on the R codes written by Samuel Shen Distinguished Professor,                          #
#San Diego State University, USA and were translated from R by Louis Selstad, Stephen Shen,                          #
#Gregori Clarke, and Dakota Newmann and edited by Samuel Shen.                                                       #
######################################################################################################################

#FIRST TIME Python users*****
#These package need to be installed (on the terminal or anaconda interface) before importing them below. 

#Follow this tutorial for package installation before
# https://towardsdatascience.com/importerror-no-module-named-xyz-45e4a5339e1b

#Change your file path to the folder where your downloaded data is stored
#MAC HELP: https://support.apple.com/guide/mac-help/go-directly-to-a-specific-folder-on-mac-mchlp1236/mac
#PC HELP: https://www.sony.com/electronics/support/articles/00015251
import os
# os.chdir("/Users/sshen/climmath/data")
os.chdir('/Users/HP/Documents/sshen/climmath/data')

%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}

#Style Dictionary to standardize plotting scheme between different python scripts 
import matplotlib.pyplot as plt

styledict = {'xtick.labelsize':20,
             'xtick.major.size':9,
             'xtick.major.width':1,
             'ytick.labelsize':20,
             'ytick.major.size':9,
             'ytick.major.width':1,
             'legend.framealpha':0.0,
             'legend.fontsize':15,
             'axes.labelsize':20,
             'axes.titlesize':25,
             'axes.linewidth':2,
             'figure.figsize':(12,8),
             'savefig.format':'jpg'}
plt.rcParams.update(**styledict)

#Function that creates personalized discrete Colormap
import numpy as np 
from matplotlib import cm as cm1
from matplotlib.colors import ListedColormap, to_rgba

def newColMap(colors):
    """
    This function creates a new color map from a list of colors given
    as a parameter. Recommended length of list of colors is at least 6.
    """
    first = np.repeat([to_rgba(colors[0])], 2, axis = 0)
    last = np.repeat([to_rgba(colors[-1])], 2, axis = 0)
    v = cm1.get_cmap('viridis', 16*(len(colors)-2))
    newcolors = v(np.linspace(0, 1, 16*(len(colors)-2)))
    for (i, col) in enumerate(colors[1:-1]):
        newcolors[16*i : 16*(i+1), :] = to_rgba(col)
    return ListedColormap(np.append(np.append(first,newcolors, axis=0), last, axis=0))

---

Chapter 10: Advanced Python Analysis and Plotting: EOFs, Trends, and Global Data

---

import numpy as np 
import matplotlib
from matplotlib import cm as cm1
import matplotlib.pyplot as plt
import pandas as pd
import warnings
import os

# import cartopy.crs as ccrs
# import cartopy.feature as cfeature

from scipy import optimize as opt
import sympy as sm
import math as m
import cartopy 
import netCDF4 as nc
from matplotlib.colors import LinearSegmentedColormap
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic
from matplotlib.patches import Polygon

from datetime import datetime

from scipy.ndimage import gaussian_filter
from sklearn import datasets, linear_model
from sklearn.neighbors import KernelDensity as kd
from sklearn.linear_model import LinearRegression
from pylab import *
warnings.filterwarnings("ignore")

---

10.2: 2Dim Spatial Domain EOFs and 1Dim Temporal PCs

---

# create two arrays, x and y, ranging from 0 to 2*pi 
# create a 3D array full of zeros for use later
x = y = np.linspace(0.,2*np.pi,100)
mydat = np.zeros((100,100,10))

# define a function that iterates term by term for the two variables into a specific function
def zfunc(x,y):
    # the triple nested for loop accomplishes term by term order we are looking for
    for t in range(0,10):
        for i in range(0,100):
            for j in range(0,100):
                mydat[i,j,t] = (np.sin(t+1)*(1/np.pi)*np.sin(x[i])*np.sin(y[j]) 
                                + np.exp(-0.3*(t+1))*(1/np.pi)*np.sin(8*x[i])*np.sin(8*y[j]))
    # print the shape and maximum of the result
    print(np.shape(mydat))
    print(np.max(mydat))
    return mydat


# plug in our x and y arrays 
mydata = zfunc(x,y)

(100, 100, 10)
0.4909421672431401

#Fig 10.1
# set a range for the contour plot
maxi = np.linspace(-.5,.5,100)

# define the contour plot using previous x, y, and mydata arrays
cmap = cm.get_cmap('hsv', 20)
CS = plt.contourf(x,y,mydata[:,:,0],maxi,cmap=cmap)

# add titles for the graph and axes
plt.title('Original field at t = 1', pad = 10)
plt.xlabel('x')
plt.ylabel('y')

# add a color bar with a label
cbar = plt.colorbar(CS, ticks=[0.4,0.2,0,-0.2,-0.4])
plt.text(6.55, 6.45, "Scale", size=20)

# show figure
plt.show()

# set a range for the contour plot
maxi = np.linspace(-.2,.2,100)

# define the contour plot using previous x, y, and mydata arrays
cmap = cm.get_cmap('hsv', 20)
# CS = plt.contourf(x,y,np.flip(mydata[:,:,10], axis=1),maxi,cmap=cmap)
CS = plt.contourf(x,y,mydata[:,:,9],maxi,cmap=cmap)

# add titles for the graph and axes
plt.title('Original field at t = 10', pad = 10)
plt.xlabel('x')
plt.ylabel('y')

# add a color bar with a label
cbar = plt.colorbar(CS, ticks=[0.2,0.1,0,-0.1,-0.2])
plt.text(6.5, 6.5, "Scale", size=20)

# save figure as a .png file
plt.show()

# define variables
x = y = np.linspace(0.,2.*np.pi,100)

# create a matrix of zeros of size (10000,10)
da = np.matrix(np.zeros((np.size(x)*np.size(y),10)))

# define a function that reshapes mydata into a new array of size
def twod_func(x,y,a):
    for i in range(0,10):
        a[:,i] = np.reshape(mydata[:,:,i],(10000,1))
    return a
    
da1 = twod_func(x,y,da)
# SVD for EOFs
u, s, vh = np.linalg.svd(da1,full_matrices=False) 
v = np.transpose(vh)

# to show the accuracy of the svd
x_a = np.dot(np.dot(u,np.diag(s)),vh)

# create an array of zeros of size (100,100)
u2 = np.zeros((100,100))

# define a function that takes all the entries in the first column
# of the u matrix, transposes it and reshapes it into size (100,100)
def dig():
    y1 = np.transpose(u[:,0])
    u2 = np.reshape(-y1,(100,100))
    return u2

def dig2():
    y1 = np.transpose(u[:,1])
    u2 = np.reshape(-y1,(100,100))
    return u2

# create zero matrix, as well as u-mat, from the dig fuction and u-mat2 from the dig2 fuction
zeros = np.zeros((100,100))
u_mat = dig()
u_mat2 = dig2()

# define two more functions that call upon the above z functions to generate spatial data
def fcn1(x,y):
    for i in range(0,100):
        for j in range(0,100):
            zeros[i,j] = z1(x[i],y[j])
    return zeros

def fcn2(x,y):
    for i in range(0,100):
        for j in range(0,100):
            zeros[i,j] = z2(x[i],y[j])
    return zeros

# define two different functions to be used below
def z1(x,y):
    return (1./np.pi)*np.sin(x)*np.sin(y)

def z2(x,y):
    return (1./np.pi)*np.sin(8.*x)*np.sin(8.*y)

# create desired Colormap
colors = ["crimson","red","red","red","tomato","lightpink",
          "pink","peachpuff","yellow","yellow", "yellowgreen","limegreen","limegreen","white","lavenderblush",
          "darkorchid","indigo","indigo","rebeccapurple","darkmagenta","mediumvioletred","mediumvioletred","black"]
colors.reverse()
cmap = newColMap(colors)

colors1 = ["crimson","red","tomato","lightpink",
          "pink","peachpuff","yellow", "yellowgreen","limegreen","white","lavenderblush",
          "darkorchid","indigo","rebeccapurple","darkmagenta","mediumvioletred","black"]
colors1.reverse()
cmap1 = newColMap(colors1)

#Fig 10.2
# graph the EOF
size = np.linspace(-2.5e-2,2.5e-2,25)
bounds = np.array([-0.02,-0.01,0.0,.01,0.02])
bounds_d = np.array([-0.3,-0.2,-0.1,0,0.1,0.2,0.3])

# set up the figure
plt.subplot(221)

# plot the data
dig = plt.contourf(x,y,u_mat,size,cmap=cmap)

# set labels and title
plt.title('SVD Mode 1: EOF1')
plt.ylabel('y')
plt.xticks([0,1,2,3,4,5,6], size=14)
plt.yticks([0,1,2,3,4,5,6],rotation=90, size=14)

# print colorbar
cbar = plt.colorbar(dig, ticks = bounds)
cbar.ax.tick_params(labelsize='large')

# create second plot 
plt.subplot(223)
lvl = np.linspace(-.35,.35,35)

# plot data and set labels
contour = plt.contourf(x,y,fcn1(x,y),lvl,cmap=cmap1)

plt.title('Accurate Mode 1')
plt.xlabel('x')
plt.ylabel('y')
plt.yticks([0,1,2,3,4,5,6], rotation=90, size=14)
plt.xticks([0,1,2,3,4,5,6], size=14)

cbar = plt.colorbar(contour, ticks = bounds_d)
cbar.ax.tick_params(labelsize='large')

# plot third plot
plt.subplot(222)
dig = plt.contourf(x,y,np.flip(u_mat2,axis=1),size,cmap=cmap)

plt.title('SVD Mode 2: EOF2')
plt.ylabel('y')
plt.yticks([0,1,2,3,4,5,6],rotation=90, size=14)
plt.xticks([0,1,2,3,4,5,6], size=14)

cbar = plt.colorbar(dig, ticks = bounds)
cbar.ax.tick_params(labelsize='large')

# plot fourth plot
plt.subplot(224)
lvl = np.linspace(-.35,.35,35)
contour = plt.contourf(x,y,fcn2(x,y),lvl,cmap=cmap1)

plt.title('Accurate Mode 2')
plt.xlabel('x')
plt.ylabel('y')
plt.yticks([0,1,2,3,4,5,6],rotation=90, size=14)
plt.xticks([0,1,2,3,4,5,6], size=14)

plt.tight_layout()

cbar = plt.colorbar(contour, ticks = bounds_d)
cbar.ax.tick_params(labelsize='large')

plt.show()

#Fig 10.3
# create variable
t = np.array(range(1,11))

# plot separate lines based on different formulas
plt.plot(t,v[:,0],color='k',marker='o',label='PC1: 83% Variance')
plt.plot(t,v[:,1],color='r',marker='o',label='PC2: 17% Variance')
plt.plot(t,-np.sin(t),color='b',marker='o',label='Original Mode 1 Coefficient: 91% Variance')
plt.plot(t,-np.exp(-.3*t),color='m',marker='o',label='Original Mode 2 Coefficient: 9% Variance')

# set labels
plt.ylim(-1.1,1.1)

plt.xlabel('Time')
plt.ylabel('PC or Coefficient')
plt.yticks([1,0.5,0,-0.5,-1])
plt.title('SVD PCs vs. Accurate Temporal Coefficients')

# show legend and plot 
plt.legend(loc = 'upper left', fontsize=20)
plt.show()

# create arrays to compute a dot product
t = np.array(range(1,11))

# make sine function
def sinarr(t):
    rslt = np.zeros(t)
    for ii in range(1,t+1):
        rslt[ii-1] = -np.sin(ii)
    return rslt

# make exponential function
def exparr(t):
    rslt = np.zeros(t)
    for jj in range(1,t+1):
        rslt[jj-1] = -np.exp(-.3*jj)
    return rslt

# create variables using above functions
sin = sinarr(10)
exp = exparr(10)
print('Dot Product:',np.dot(np.transpose(sin),exp))

# find the variances and compare them
v1 = np.var(sin)
v2 = np.var(exp)
print("Variance:",v1/(v1+v2))

# find the product of the PC's
v_1 = np.transpose(v[:,0])
v_2 = v[:,1]
print("Dot Product of PCs:" , np.dot(v_1,v_2))

Dot Product: 0.8625048359266784
Variance: 0.9098719287510537
Dot Product of PCs: [[2.08166817e-17]]

# create arrays that will be used below 
# and check the sizes of the arrays for matrix multiplication
sdiag = np.diag(s)
vtran = np.transpose(v[:,0:1])

# create two new arrays, B and B1 multiplying the previous
# matrices together, being careful to match their sizes
B = np.transpose(np.matmul(np.matmul(u[:,0:1],sdiag[0:1,0:1]),vtran))
B1 = np.matmul(np.matmul(u,sdiag),np.transpose(v))

#Fig 10.4
# reshaping B to plot it
BB = np.reshape(B[4,:],(100,100))
scale = np.linspace(-.4,.4,25)
bounds = [0.4,0.2,0,-0.2,-0.4]

cmap = cm.get_cmap('hsv', 15)

# set up figure
fig = plt.figure(figsize=(14,6))

# create subplot
plt.subplot(121)
BBplot = plt.contourf(x,y,BB,scale,cmap = cmap)

# set labels 
plt.title('(a) 2-mode SVD reconstructed field t = 5', size=18)
plt.xlabel('x')
plt.ylabel('y')
plt.xticks([0,1,2,3,4,5,6])

# print colorbar
cbar = plt.colorbar(BBplot, ticks=bounds)
cbar.ax.tick_params(labelsize='large')


# reshaping B to plot it
new_B = B1.T
BB = np.reshape(new_B[4,:],(100,100))

# create subplot
plt.subplot(122)
BBplot = plt.contourf(x,y,BB,scale,cmap = cmap)

# set labels
plt.title('(b) All-mode SVD reconstructed field t = 5', size=18)
plt.xlabel('x')
plt.ylabel('y')
plt.xticks([0,1,2,3,4,5,6])

# create colorbar
cbar = plt.colorbar(BBplot, ticks=bounds)
cbar.ax.tick_params(labelsize='large')

plt.tight_layout()

---

10.3 From Climate Data Download to EOFs and PC Visualization: An NCEP/NCAR Reanalysis Example

---

#Download the following data set from NOAA into your active folder
# https://www.esrl.noaa.gov/psd/repository/entry/show?
# entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%
# 3AL25jZXAucmVhbmFseXNpcy5kZXJpdmVkL3N1cmZhY2UvYWlyLm1vbi5tZWFuLm5j

ncd = nc.Dataset("air.mon.mean.nc","r+")

#Define variables
lon_vals = ncd.variables['lon']
lat_vals = ncd.variables['lat']
time = ncd.variables['time']
air = ncd.variables['air']
time_unit = time.units
precnc = ncd.variables['air']

# understanding the time range of the data
start = time.actual_range[0]
end = time.actual_range[1]
hrs_pr_yr = 8760.812777
start_yrs_after_1800 = start/hrs_pr_yr
end_yrs_after_1800 = end/hrs_pr_yr
print(f"The data starts at year {round(start_yrs_after_1800 + 1800, 2)} and ends at about {round(end_yrs_after_1800 + 1800, 2)}")

The data starts at year 1948.08 and ends at about 2016.87

# Fig 10.5
# plot a 90S-90N temp along a meridional line at 180E

# create variables
nc_x = np.linspace(-90,90,73)
y = list(precnc[0,:,71])
y.reverse()

# plot data
plt.plot(nc_x,y,marker="o",color='k' )

# set labels and show figure
plt.xlabel('Latitude')
plt.ylabel('Temperature [$\degree$C]')
plt.yticks(rotation=90)
plt.title('NCEP/NCAR Reanalysis Surface Air Temperature [$\degree$C] \n along a Meridional Line at 180$\degree$E: Jan. 1948'
         , pad = 10)
plt.grid()
plt.show()

list(precnc[0,:,71]).reverse()

precnc.shape

(826, 73, 144)

# define variables 
precst = np.zeros((10512,826))

temp = np.reshape(precnc[0,:,:],(144*73))
print(np.shape(temp))
print(temp[0:6])

# create reshaping function
def spmat(x,y):
    for i in range(0,826):
        y[:,i] = np.reshape(x[i,:,:],(144*73))
    return y

# use function and save result
precst2 = spmat(precnc,precst)

print(np.shape(precst2))

# build lat and lon for 10512 spatial positions using rep
def rep(x,y,n):
    for j in range(0,n):
        x = np.append(x,y)
    return x

arra = np.zeros((0))
LAT = rep(arra,lat_vals,144)
LON = rep(arra,lon_vals[0],73)
print(LON)

# create reconstruction function
def recon():
    rslt = LON
    for jj in range(0,143):
        rslt = np.append(rslt,rep(arra,lon_vals[jj],73))
    return rslt
LON2 = recon()
print("The last 6 entries are:",LON2[-6:])

# creating an array with first two columns of lat and lon
# with 826 months spanning from 1948-2016
gpcpst_ = np.column_stack((LAT,(LON2)))
gpcpst = np.column_stack((gpcpst_,precst2))

(10512,)
[-34.926773 -34.926773 -34.926773 -34.926773 -34.926773 -34.926773]
(10512, 826)
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0.]
The last 6 entries are: [355. 355. 355. 355. 355. 355.]

# create an array of all ones for the days of the data
ones = np.ones(826)

# define a function to create a list of months repeating to match the size of our data
def monthfunc():
    montharr = np.zeros(0)
    for i in range(1,71):
        for j in range(1,13):
             montharr = np.append(montharr,int(j))
    for k in range(1,6):
        montharr = np.append(montharr,int(k))
    return montharr

months = monthfunc()

# define a function to create a list of years to match the size of our data
def yearfunc():
    yeararr = np.zeros(0)
    for i in range(1948,2018):
        for j in range(1,13):
            yeararr = np.append(yeararr,i)
    for k in range(1,6):
        yeararr = np.append(yeararr,2018)
    return yeararr

years = yearfunc()

# define a function that pairs a string version of our year and month arrays with a dash
# in between these will serve as our column names for our data
def tmfunc():
    empty = np.zeros(0)
    for i in range(0,826):
        empty = np.append(empty,str(years[i])  + '-'  + str(months[i]))
    return empty

tm1 = tmfunc()

# add the 'Lat' and 'Lon' columns to the beginning
# of our column names array
tm2 = np.append(['Lat','Lon'],tm1)

# use the package pandas.DataFrame to match the column names with the data
GPCPST = pd.DataFrame(gpcpst,columns=tm2)

monJ = np.array(range(0,826,12))
gpcpdat = GPCPST.iloc[:,2:826]
gpcpJ = gpcpdat.iloc[:,monJ]
climJ = gpcpJ.mean(axis=1)

# create desired functions
def rowSDS2():
    rslt = []
    for jj in range(0,10512):
        rslt = np.append(rslt,np.std(gpcpJ.iloc[jj,:]))
    return rslt

sdJ = gpcpJ.std(axis=1)

def anomJ_func2():
    rslt = pd.DataFrame(np.zeros((10512, 69)),columns=tm1[monJ])
    for i in range(0,69):
        rslt.iloc[:,i] = (gpcpdat.iloc[:,monJ[i]]-climJ)/sdJ
    return rslt

anomJ = anomJ_func2()

def anomJ_func3():
    rslt = pd.DataFrame(np.zeros((10512, 69)),columns=tm1[monJ])
    for i in range(0,69):
        rslt.iloc[:,i] = np.sqrt(np.cos(GPCPST.iloc[:,0]*np.pi/180.))*anomJ.iloc[:,i]
    return rslt
        
anomAW = anomJ_func3()

# computing SVD from the area-weighted anomaly matrix anomAW
U, S, Vh = np.linalg.svd(anomAW,full_matrices=False)

print(U.shape,S.shape,Vh.shape)

(10512, 69) (69,) (69, 69)

#Fig 10.6
# create variables
xvals = np.linspace(0.,70.,69)
yvals = 100.*(S**2.)/np.sum(S**2)

# create two subplots
fig, ax1 = plt.subplots()

# plot of first plot
ax1.scatter(xvals,yvals,marker = "o",facecolors='none',edgecolors='k')
ax1.plot(xvals,yvals,color='k',label = "Percentage Variance")
ax1.set_ylabel('Percentage Variance [%]',color='k')
ax1.set_xlabel('Mode number',color='k')
ax1.tick_params(axis='y',labelcolor='k')
ax1.grid()

# twin the axis to prepare second plot
ax2 = ax1.twinx()

# plot data on second axis 
ax2.scatter(xvals,yvals.cumsum(),marker = "o",facecolors='none',edgecolors='b')
ax2.plot(xvals,yvals.cumsum(),color='b',label = "Cumulative Variance")

# set labels
ax2.set_ylabel('Cumulative Variance [%]',color='b')
ax2.tick_params(axis='y',labelcolor='b')

ax1.legend(loc=(.5,.5))
ax2.legend(loc=(.5,.45))
plt.title("Eigenvalues of Covariance Matrix")

# show figure
plt.show()

# define function
def mapmatr():
    q = np.zeros((144,73))
    for ii in range(0,144):
        q[ii,:] = U[:,0]/(np.sqrt(np.cos(gpcpst[:,0]
          *np.pi/180.)))
    return q

# reshape data using result of SVD 
mapmatr = (np.reshape(U[:,0]/(np.sqrt(np.cos(gpcpst[:,0]
          *np.pi/180.))),(73,144)))

mapmatr2 = mapmatr[:,range(0,mapmatr[0,:].size)]

# make Colormap and create desired Colormap

# colors = ["darkred","red","darkorange","orange","yellow","lightyellow","white","lightgreen"
#           ,"green","blue","darkblue", "black"]
# myColMap = newColMap(colors)
myColMap = LinearSegmentedColormap.from_list(name='my_list', colors=['blue',
                                            'green','white','yellow','red'], N=256)

#Fig 10.7
# prepare figure
plt.figure(figsize=(12,12))

# plot EOF1
plt.subplot(211)

# using the basemap package to put the cylindrical projection map into the image
eof = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data because 
lons2, data2 = eof.shiftdata(lon_vals, datain = -mapmatr2, lon_0=180)

# draw coastlines, latitudes, and longitudes on the map
eof.drawcoastlines(color = 'black', linewidth = 1)
eof.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=15)
eof.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=15)
eof.drawcountries()

limits = np.linspace(-.04,.04,51)

# plot data on contour map
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data2,limits,cmap = myColMap)

# set labels and colorbar bounds
bounds = np.array([-0.04,-0.02,0,0.02,0.04])
efo = plt.colorbar(eof_plt, ticks = bounds, shrink = 0.85)
efo.ax.tick_params(labelsize='large')
plt.text(375, 93, 'Scale', size=16)
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =40)
plt.title('January EOF1 from 1948-2017 NCEP Temp. Data', pad = 10)

# retrieve PC1 from SVD 
pcdat = Vh[0,:]

# open bottom plot
plt.subplot(212)

#create time variable
TIME = np.array(range(1948,2017))

# plot data
plt.scatter(TIME,-pcdat,marker = "o",facecolors='none',edgecolors='k')
plt.plot(TIME,-pcdat,color='k')
plt.title("PC1 of NCEP RA Jan. SAT: 1948-2017", pad = 10)
plt.xlabel("Year")
plt.ylabel("PC Values")
plt.grid()
plt.tight_layout()
plt.show()

#Fig 10.8
# prepare figure
plt.figure(figsize=(12,12))

# plot EOF2 by reshaping data from SVD
mapmatr = (np.reshape(U[:,1]/(np.sqrt(np.cos(gpcpst[:,0]
          *np.pi/180.))),(73,144)))

mapmatr2 = mapmatr[:,range(0,mapmatr[0,:].size)]

# open first subplot
plt.subplot(211)

# using the basemap package to put the cylindrical projection map into the image
eof = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data
lons2, data2 = eof.shiftdata(lon_vals, datain = -mapmatr2, lon_0=180)

# draw coastlines, latitudes, and longitudes on the map
eof.drawcoastlines(color='black', linewidth=1)
eof.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=15)
eof.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=15)
eof.drawcountries()

# create contour limits
limits = np.linspace(-.04,.04,51)

# plot data
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data2,limits,cmap = myColMap)

# set labels and colorbar bounds
bounds = np.array([-0.04,-0.02,0,0.02,0.04])
efo = plt.colorbar(eof_plt, ticks=bounds, shrink = 0.85)
efo.ax.tick_params(labelsize='large')
plt.text(375, 93, 'Scale', size=16)
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =40)
plt.title('January EOF2 from 1948-2017 NCEP Temp. Data', pad = 10)

pcdat = Vh[1,:]

# open bottom plot
plt.subplot(212)

# create time variable
TIME = np.array(range(1948,2017))

# plot PC2
plt.scatter(TIME,-pcdat,marker = "o",facecolors='none',edgecolors='k')
plt.plot(TIME,-pcdat,color='k')
plt.title("PC2 of NCEP RA Jan. SAT: 1948-2017", pad = 10)
plt.xlabel("Year")
plt.ylabel("PC Values")
plt.grid()
plt.tight_layout()
plt.show()

#Fig 10.9
# set figure up
plt.figure(figsize=(12,12))

# create data from SVD 
mapmatr = (np.reshape(U[:,2]/(np.sqrt(np.cos(gpcpst[:,0]
          *np.pi/180.))),(73,144)))
mapmatr2 = mapmatr[:,range(0,mapmatr[0,:].size)]

# open first plot
plt.subplot(211)

# using the basemap package to put the cylindrical projection map into the image
eof = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data because 
lons2, data2 = eof.shiftdata(lon_vals, datain = -mapmatr2, lon_0=180)

# draw coastlines, latitudes, and longitudes on the map
eof.drawcoastlines(color='black', linewidth=1)
eof.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=15)
eof.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=15)
eof.drawcountries()

limits = np.linspace(-.04,.04,51)

# plot data 
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data2,limits,cmap = myColMap)

# set labels and colorbar bounds
bounds = np.array([-0.04,-0.02,0,0.02,0.04])
efo = plt.colorbar(eof_plt, ticks=bounds, shrink = 0.85)
efo.ax.tick_params(labelsize='large')
plt.text(375, 93, 'Scale', size=16)
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =40)
plt.title('January EOF3 from 1948-2017 NCEP Temp. Data', pad = 10)

pcdat = Vh[2,:]

# open bottom plot
plt.subplot(212)

# create time variable
TIME = np.array(range(1948,2017))

# plot data 
plt.scatter(TIME,-pcdat,marker = "o",facecolors='none',edgecolors='k')
plt.plot(TIME,-pcdat,color='k')

# set labels 
plt.title("PC3 of NCEP RA Jan. SAT: 1948-2017", pad = 10)
plt.xlabel("Year")
plt.ylabel("PC Values")
plt.grid()
plt.tight_layout()
plt.show()

GPCPST.shape

(10512, 828)

# EOF from de-trended data
monJ = np.array(range(0,815,12))
gpcpdat = GPCPST.iloc[:,2:817]
gpcpJ = gpcpdat.iloc[:,monJ]
climJ = gpcpJ.mean(axis=1)
sdJ = gpcpJ.std(axis=1)

anomJ = ((gpcpdat.iloc[:,monJ].T-climJ)/sdJ).T

# create linear regression object
regr = linear_model.LinearRegression()
print(anomJ.iloc[0,:].shape)

(68,)

# create time variable 
TIME = np.array(range(1948,2016))

# create trend V
def trendV():
    trendV = np.zeros(10512)
    for i in range(0,10512):
        regr.fit(np.reshape(TIME,(68,1)),
                 np.reshape(np.array(anomJ.iloc[i,:]),(68,1)))
        trendV[i] = regr.coef_
    return trendV
 
# create trend M
def trendM():
    rslt = np.zeros((10512,68))
    trendM = pd.DataFrame(np.zeros((10512, 68)),columns=tm1[monJ])
    for i in range(0,10512):
        regr.fit(np.reshape(TIME,(68,1)),
                 np.reshape(np.array(anomJ.iloc[i,:]),(68,1)))
        rslt[i,:] = regr.intercept_
    return rslt

# make trend V instance 
trendV = trendV()
print(trendV)

[ 0.02303803  0.02303803  0.02303803 ... -0.03220623 -0.03220623
 -0.03220623]

# make trend M instance
trendM = trendM()
print(trendM)

[[-45.64986127 -45.64986127 -45.64986127 ... -45.64986127 -45.64986127
  -45.64986127]
 [-45.64986127 -45.64986127 -45.64986127 ... -45.64986127 -45.64986127
  -45.64986127]
 [-45.64986127 -45.64986127 -45.64986127 ... -45.64986127 -45.64986127
  -45.64986127]
 ...
 [ 63.81665437  63.81665437  63.81665437 ...  63.81665437  63.81665437
   63.81665437]
 [ 63.81665437  63.81665437  63.81665437 ...  63.81665437  63.81665437
   63.81665437]
 [ 63.81665437  63.81665437  63.81665437 ...  63.81665437  63.81665437
   63.81665437]]

print(trendV.shape)
print(trendM.shape)

(10512,)
(10512, 68)

dtanomJ = anomJ - trendM
dtanomJ.shape

(10512, 68)

vArea = np.cos(gpcpst[:,0]*np.pi/180.)
dtanomAW = np.sqrt(vArea)*dtanomJ.T

# computing SVD from the area-weighted anomaly matrix anomAW
U, S, Vh = np.linalg.svd(dtanomAW.T,full_matrices=False)

print(U.shape,S.shape,Vh.shape)

(10512, 68) (68,) (68, 68)

# define function
def mapmatr():
    q = np.zeros((144,73))
    for ii in range(0,144):
        q[ii,:] = U[:,0]/(np.sqrt(np.cos(gpcpst[:,0]
          *np.pi/180.)))
    return q

#Fig 10.10
# plot EOF1
# prepare figure
plt.figure(figsize=(12,12))
plt.subplot(211)

# create data using SVD 
mapmatr = (np.reshape(U[:,1]/(np.sqrt(np.cos(gpcpst[:,0]
          *np.pi/180.))),(73,144)))

mapmatr2 = mapmatr[:,range(0,mapmatr[0,:].size)]

# using the basemap package to put the cylindrical projection map into the image
eof = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data because 
lons2, data2 = eof.shiftdata(lon_vals, datain = mapmatr2)

# draw coastlines, latitudes, and longitudes on the map
eof.drawcoastlines(color='black', linewidth=1)
eof.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=15)
eof.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=15)
eof.drawcountries()

limits = np.linspace(-.04,.04,51)

# plot EOF1 on contour map
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data2,limits,cmap = myColMap)

# set labels and colorbar
bounds = np.array([-0.04,-0.02,0,0.02,0.04])
efo = plt.colorbar(eof_plt, ticks=bounds, shrink = 0.85)
efo.ax.tick_params(labelsize='large')
plt.text(375, 93, 'Scale', size=16)
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =40)
plt.title('January EOF1 from 1948-2017 NCEP \n De-Trended Standardized Temp. Anomaly Data', pad = 10)

# retrieve PC1 from SVD 
pcdat = Vh[1,:]

# open bottom plot
plt.subplot(212)

# plot PC1
TIME = np.array(range(1948,2016))
plt.scatter(TIME,pcdat,marker = "o",facecolors='none',edgecolors='k')
plt.plot(TIME,pcdat,color='k')

# set labels 
plt.title("PC1 of NCEP RA Jan. De-Trended Standardized SAT: 1948-2017", pad = 10)
plt.xlabel("Year")
plt.ylabel("PC Values")
plt.grid()
plt.tight_layout()
plt.show()

#Fig 10.11
# prepare figure
plt.figure(figsize=(12,12))

# create data from SVD 
mapmatr = (np.reshape(U[:,2]/(np.sqrt(np.cos(gpcpst[:,0]
          *np.pi/180.))),(73,144)))

mapmatr2 = mapmatr[:,range(0,mapmatr[0,:].size)]

plt.subplot(211)

# using the basemap package to put the cylindrical projection map into the image
eof = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data 
lons2, data2 = eof.shiftdata(lon_vals, datain = mapmatr2)

# draw coastlines, latitudes, and longitudes on the map
eof.drawcoastlines(color='black', linewidth=1)
eof.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=15)
eof.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=15)
eof.drawcountries()

limits = np.linspace(-.04,.04,51)

# plot EOF2 
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),-data2,limits,cmap = myColMap)

# set labels and colorbar
bounds = np.array([-0.04,-0.02,0,0.02,0.04])
efo = plt.colorbar(eof_plt, ticks=bounds, shrink = 0.85)
efo.ax.tick_params(labelsize='large')
plt.text(375, 93, 'Scale', size=16)
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =40)
plt.title('January EOF2 from 1948-2017 NCEP \n De-Trended Standardized Temp. Anomaly Data', pad = 10)

# retrieve PC2 data
pcdat = Vh[2,:]

# open bottom plot
plt.subplot(212)

# plot PC1
TIME = np.array(range(1948,2016))
plt.scatter(TIME,-pcdat,marker = "o",facecolors='none',edgecolors='k')
plt.plot(TIME,-pcdat,color='k')

# set labels 
plt.title("PC2 of NCEP RA Jan. Standardized SAT: 1948-2017", pad = 10)
plt.xlabel("Year")
plt.ylabel("PC Values")
plt.grid()
plt.tight_layout()
plt.show()

---

10.4 Area-Weighted Average and Spatial Distribution of Trend

---

# EOF from de-trended data
monJ = np.array(range(0,816,12))
gpcpdat = GPCPST.iloc[:,2:817]
gpcpJ = gpcpdat.iloc[:,monJ]
climJ = gpcpJ.mean(axis=1)
sdJ = gpcpJ.std(axis=1)

# create desired variables for de-trended data
anomJ = ((gpcpdat.iloc[:,monJ].T-climJ)/sdJ).T
vArea = np.cos(gpcpst[:,0]*np.pi/180.)
anomA = vArea*anomJ.T
anomA.shape

(68, 10512)

# compute the de-trended data
JanSAT = anomA.T.sum(axis=0)/vArea.sum()
JanSAT.shape

(68,)

#Fig 10.12
# set up figure 
plt.figure(figsize=(12,7))

# create time variable
TIME = np.linspace(1948,2017, 68)

# create linear regression object
lm = linear_model.LinearRegression()
lm.fit(TIME.reshape(-1,1), JanSAT)
pred = lm.predict(TIME.reshape(-1,1))

# plot linear regression
plt.plot(TIME, pred, "-",color="red", linewidth=2)

# plot data 
plt.plot(TIME, JanSAT, "-ok")

# set labels 
plt.title("Global Average Jan. SAT Anomalies from NCEP RA", pad = 10)
plt.text(1950,0.35,"Linear Trend ${}\degree$C per Century".format(np.round(lm.coef_,3)[0]*100), color="r", fontsize=18)
plt.ylabel("Temperature [$\degree$C]")
plt.xlabel("Year")
plt.grid()
plt.show()

# recreate desired variables 
monJ = np.array(range(0,826,12))
gpcpdat = GPCPST.iloc[:,2:826]
gpcpJ = gpcpdat.iloc[:,monJ]
climJ = gpcpJ.mean(axis=1)

# define function and make instance of function 
def anomJ_func3():
    rslt = pd.DataFrame(np.zeros((10512, 69)),columns=tm1[monJ])
    for i in range(0,69):
        rslt.iloc[:,i] = gpcpdat.iloc[:,monJ[i]]-climJ
    return rslt

anomJ = anomJ_func3()

# create time variable
TIME = np.linspace(1948,2017, 69)

# create trend V function and instance of trend V
def trendV():
    trendV = np.zeros(10512)
    for i in range(0,10512):
        regr.fit(np.reshape(TIME,(69,1)),
                 np.reshape(np.array(anomJ.iloc[i,:]),(69,1)))
        trendV[i] = regr.coef_
    return trendV

trendV = trendV()

# reshape data into desired form 
mapmat1 = np.reshape([10*trendV],(73,144))

# create and use function
def mapv1_func():
    zeros = np.zeros((73,144))
    for i in range(0,73):
        for j in range(0,144):
            if mapmat1[i,j] >= 1.:
                zeros[i,j] = 1.
            else:
                zeros[i,j] = mapmat1[i,j]
    return zeros
mapv1 = mapv1_func()

# create and use function
def mapmat_func():
    zeros = np.zeros((73,144))
    for i in range(0,73):
        for j in range(0,144):
            if mapv1[i,j] <= -1.:
                zeros[i,j] = -1.
            else:
                zeros[i,j] = mapv1[i,j]
    return zeros

mapmat3 = mapmat_func()

#Fig 10.13
# prepare figure
plt.figure(figsize=(14., 7.))

# using the basemap package to put the cylindrical projection map into the image
trend = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data because 
lons2, data3 = trend.shiftdata(lon_vals, datain = mapmat3, lon_0=180)

# draw coastlines, latitudes, and longitudes on the map
trend.drawcoastlines(color='black', linewidth=1)
trend.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=17)
trend.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=17)
trend.drawcountries()

limits = np.linspace(-1.1,1.1,50)

# plot data 
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data3,limits,cmap = myColMap)

# set labels and colorbar
efo = plt.colorbar(eof_plt, ticks=[1,0.5,0,-0.5,-1], shrink = 0.80)
plt.text(375, 95, 'Scale', size=16)
plt.xlabel('Longitude',labelpad = 40)
plt.ylabel('Latitude',labelpad = 65)
plt.title('Trend of the NCEP RA1 Jan. 1948-2017 \n Anomaly Temp. [$\degree$C/Century]', pad = 8)
plt.show()

---

10.5 GPCP Precipitation Data: Analysis and Visualization by Python

---

#Download the following data set from NOAA into your active folder
#https://www.esrl.noaa.gov/psd/repository/entry/show?
# entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%
# 3AL25jZXAucmVhbmFseXNpcy5kZXJpdmVkL3N1cmZhY2UvYWlyLm1vbi5tZWFuLm5j

ncd = nc.Dataset("precip.mon.mean.nc","r+")

# define variables
lon_vals = ncd.variables['lon'][:]
lat_vals = ncd.variables['lat'][:]
time = ncd.variables['time']
time_unit = time.units
precnc = ncd.variables['precip']

#Fig 10.14
# plot a 90S-90N temp along a meridional line at 180E
nc_x = np.linspace(-90,90,72)
plt.plot(nc_x,precnc[0,:,63],color='k' )

# set variables
plt.xlabel('Latitude')
plt.ylabel('Precipitation [mm/day]')
plt.title('90S-90N Precipitation along a Meridional Line at 160E: Jan. 1979', size=20)
plt.grid()
plt.show()

# write the data as space-time matrix with a header
precst = np.zeros((10368,451))
temp = precnc[0,:,:]

# fill precst with data from precnc 
for i in range(451):
    precst[:,i] = precnc[i,:,:].flatten()

# create and reshape Lat and Lon variable for DataFrame
LAT = list(lat_vals.data) * 144
LON1 = [np.repeat(i,72) for i in lon_vals.data]
LON = [i for j in LON1 for i in j]

gpcpst_ = pd.DataFrame({"Lat":LAT,"Lon":LON})
gpcpst = pd.concat([gpcpst_, pd.DataFrame(precst)], axis=1)

gpcpst.iloc[889:892,:5]

	Lat	Lon	0	1	2
889	-26.25	31.25	0.057521	0.180531	0.210445
890	-23.75	31.25	0.108556	0.256279	0.288389
891	-21.25	31.25	0.087190	0.273805	0.256987

# create desired Colormap
colors = ["black", "black","purple","crimson","red","red","red","tomato","lightpink",
          "pink","peachpuff","yellow","yellow", "yellowgreen","limegreen","limegreen","mediumaquamarine",
          "skyblue","lightskyblue","wheat","wheat"]
colors.reverse()
myColMap = newColMap(colors)

#Fig 10.15(a)
# prepare variables for plotting 
climmat = precnc[0,:,:]
for i in range(451):
    climmat = climmat + precnc[i,:,:]

climmat = climmat/451
mapmat = climmat

# suppress values into the (0,10) range 
for i in range(72):
    for j in range(144):
        if mapmat[i,j] > 10:
            mapmat[i,j] = 10
        if mapmat[i,j] < 0:
            mapmat[i,j] = 0

myint = np.linspace(0, 10,11)

# prepare figure
plt.figure(figsize=(15., 7.))

# using the basemap package to put the cylindrical projection map into the image
trend = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data because 
lons2, data3 = trend.shiftdata(lon_vals, datain = mapmat, lon_0=180)

# draw coastlines, latitudes, and longitudes on the map
trend.drawcoastlines(color='black', linewidth=1)
trend.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=15)
trend.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=15)
trend.drawcountries()

# plot data 
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),mapmat,myint,cmap = myColMap)

# set labels and colorbar
efo = plt.colorbar(eof_plt, ticks=[10,9,8,7,6,5,4,3,2,1], shrink = 0.85)
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =40)
plt.title('GPCP 1979-2016 Precipitation Climatology [mm/day]', size=18, fontweight="bold")
plt.text(367, 94, '[mm/day]', size = 15)
plt.show()

mapmat.shape

(72, 144)

# create and prepare data 
sdmat = (precnc[0,:,:]-climmat)**2
for i in range(1,451):
    sdmat = sdmat + (precnc[i,:,:]-climmat)**2

sdmat = np.sqrt(sdmat/451)
mapmat = sdmat

# suppress values into the (0,5) range 
for i in range(72):
    for j in range(144):
        if mapmat[i,j] > 5:
            mapmat[i,j] = 5
        if mapmat[i,j] < 0:
            mapmat[i,j] = 0

#Fig 10.15(b)
myint = np.linspace(0, 5,25)

# prepare figure
plt.figure(figsize=(15., 7.))

# using the basemap package to put the cylindrical projection map into the image
trend = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
            llcrnrlon=0.,urcrnrlon=360.,resolution='c')

# shifting the data because 
lons2, data3 = trend.shiftdata(lon_vals, datain = mapmat, lon_0=180)

# draw coastlines, latitudes, and longitudes on the map
trend.drawcoastlines(color='black', linewidth=1)
trend.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1], fontsize=15)
trend.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1], fontsize=15)
trend.drawcountries()

# plot data 
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),mapmat,myint,cmap = myColMap)

# set labels and colorbar
efo = plt.colorbar(eof_plt, ticks=[5,4.5,4,3.5,3,2.5,2,1.5,1,0.5,0], shrink = 0.85)
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =40)
plt.title('GPCP 1979-2016 Standard Deviation of Precipitation [mm/day]', size=18, fontweight="bold", pad = 10)
plt.text(367, 94, '[mm/day]', size = 15)
plt.show()