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))
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")
# 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()