####################################################################
#
#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
#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,
#Email: sshen@sdsu.edu
#www.climatemathematics.org
#March 2019, San Diego
#
####################################################################
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import rgb2hex
from matplotlib.patches import Polygon
import matplotlib.mlab as mlab
import pandas as pd
from mpl_toolkits.basemap import Basemap
from mpl_toolkits import mplot3d
from netCDF4 import Dataset as ds
from urllib import request
import scipy as sp
from scipy import stats as stt
from sklearn import datasets, linear_model
from sklearn.neighbors import KernelDensity as kd
from sklearn.linear_model import LinearRegression
from random import randint
from scipy.stats import norm
import statsmodels
import sklearn
import math
import statistics
import sympy as sy
from sympy import symbols, diff
import statsmodels.api as sm
from datetime import date
1+4 #The result 5 is shown as an Out[5] below.
2 + np.pi/4 - 0.8 #\pi=3.14...... in Python is given by np.pi
x=1
y=2
z=4
t=2*x**y-z #x**y means x^y: x to the power of y
t
u=2
v=3
u+v
np.sin(u*v) #np.sin is sine function in Python
tmax=[77, 72,75,73,66,64,59] #Enter the temperature data
tmax #Show the temperarure data
np.arange(9) #A Python sequence starts from 0 while R from 1
#A sequence is called an array in Python
# np.arange(x,y) creates an array of numbers incremented by 1
# between integer x and integer y-1
np.arange(1, 9)
np.arange(-5.5, 2) #np.arange(x,y) can take negative and real numbers with increment 1
# np.arange(x,y,b) creates an array of numbers incremented by b
# between x and y-1
#x,y, and b can be real numbers although they are integers in most applications
np.arange(1, 9, 2)
# np.linspace(x, y, n) creates an array of n numbers evenly distributed between x and y
np.linspace(0, 4.5, 5) #linspace means linearly/uniformly distributed space division
# Define the x^2 function and call it samfctn(x)
def samfctn(x):
return x * x;
samfctn(4)
# Define a multivaraite function "x+y-z/2"
def fctn2(x, y, z):
return x + y - z / 2;
fctn2(1, 2, 3)
import math
import numpy
import scipy
import sympy
#After import these packages, the following three return the same pi=3.141592...
math.pi
np.pi
scipy.pi
t = np.arange(-math.pi, 2 *math.pi, 0.1) #math.pi is the same as np.pi and scipy.pi, just using different packages
plt.plot(t, np.sin(t))
plt.show()
t = np.linspace(-3, 2, 20) #uniform division using 20 points
plt.plot(t, t ** 2)
plt.show()
np.size(t) #length of the data array
import matplotlib.pyplot as plt
t = np.linspace(-3, 2, 300) #uniform division using 200 points
plt.plot(t, (np.sin(t ** 2))**2, 'ro', color='blue',linewidth=5) #linestyle ='ro'
plt.xlabel('Time')
plt.ylabel('Function Value')
plt.title('Plot a function')
plt.show()
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
# defines a formula for Z
def f(x, y):
return 1 - x**2 - y**2
X, Y = np.meshgrid(x, y)
Z = f(X, Y) # renders z function on the x, y grid
# creates the 3D plot
fig = plt.figure()
ax = plt.axes(projection = '3d')
ax.contour3D(X, Y, Z, 30, cmap = 'binary')
plt.show()
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
# defines a formula for Z
def f(x, y):
return 1 - x**2 - y**2
X, Y = np.meshgrid(x, y)
Z = f(X, Y) # renders z function on the x, y grid
# creates the 3D plot
fig = plt.figure()
ax = plt.axes(projection = '3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
cmap='viridis', edgecolor='none') #cmap = 'color'
# labels the plot
ax.set_title('Contour surface')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
ax.view_init(30)
plt.show()
x = np.linspace(-1, 1, 30)
y = np.linspace(-1, 1, 30)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)
# creates and colors the contour plot
plt.contourf(X, Y, Z, 20, cmap = 'bwr') #cmap='color-scheme'
# adds colorbar for color map of contours
plt.colorbar()
plt.show()
## Symbolic calculations by Python
# Makes 'x' and 'y' symbols used for derivatives
x = symbols('x')
y = symbols('y')
# 'diff' command takes the derivative of x^2 w.r.t. x
diff(x ** 2, x)
# creates a function of x^2
fx = x ** 2
# takes the derivative of function w.r.t x
rslt = diff(fx, x)
print('The derivative of the function fx w.r.t x is: ', rslt)
# creates a function of x^2*sin(x)
fx = (x ** 2) * sy.sin(x)
# Take the derivative of function w.r.t x
rslt = diff(fx, x)
print('The derivative of the function w.r.t x is: ', rslt)
# Create a multivariate function: x^2+y^2
fxy = x ** 2 + y ** 2
# Take partial derivative of function w.r.t x
rslt = diff(fxy, x)
print('The partial derivative of the function w.r.t x is: ', rslt)
# Take partial derivative of function w.r.t y
rslt = diff(fxy, y)
print('The partial derivative of the function w.r.t y is: ', rslt)
# Integrates x^2 from 0 to 1
sy.integrate(x ** 2, (x, 0, 1))
# integrates cos(x) from 0 to pi/2
sy.integrate(sy.cos(x), (x, 0, math.pi / 2))
# creates a 5X1 matrix
np.mat([[1], [6], [3], [np.pi], [-3]])
# creates a sequence from 2 to 6
np.arange(2, 7)
# creates a sequence incremented by 2 from 1 to 10
np.arange(1, 10, 2)
# creates a 4X1 matrix
x = np.mat([[1], [-1], [1], [-1]])
print(x)
# adds number 1 to each element of x
print(x + 1)
# multiplies number 2 to each element of x
print(x * 2)
#x/2, divides each element of x by 2
print(x / 2)
# creates a 1X4 vector with a sequence incremented by 1 from 1 to 4
y = np.mat(np.arange(1, 5))
# this multiplication multiples each pair of elements and results in 4X1 vector
np.multiply(x, y.T)
print(x / 2)
print(x)
print(y)
# the dot product of two vectors
np.dot(y, x)
x.shape
# transforms x into a row 1X4 vector
x.T.shape
# dot product and results in 1X1 matrix
np.dot(x.T, y.T)
# this column times row yields a 4X4 matrix
np.matmul(x, y)
# convert a vector into a matrix of the same number of elements
# the matrix elements go by column, first column, second, etc
my = y.reshape(2, 2).transpose()
print(my)
# finds dimensions of a matrix
my.shape
# converts a matrix to a vector, again via columns
my.flatten()
# reshapes x to a 2X2 matrix
mx = x.reshape(2, -1).transpose()
# multiplication between each pair of elements
np.multiply(mx, my)
# division between each pair of elements
np.divide(mx, my)
np.subtract(mx, 2 * my)
# this is matrix multiplication in linear algebra
np.matmul(mx, my)
# determinant of a matrix
x = np.linalg.det(my)
print(x)
# inverse of a matrix
z = np.linalg.inv(my)
print(z)
# verifies that the inverse of the matrix is correct by multiplying
# the original matrix and the inverse matrix
a = np.matmul(z, my)
print(a)
# returns a left to right diagonal vector of a matrix
np.diagonal(my)
# returns the eigenvalues of the "my" matrix
myeigval = np.linalg.eigvals(my)
print(myeigval)
# returns the eigenvalues and eigenvectors of the "my" matrix
myeigvect = np.linalg.eig(my)
print(myeigvect)
# returns the singular value decomposition(svd) of the "my" matrix
mysvd = np.linalg.svd(my)
print(mysvd)
# returns a solved linear matrix equation of the "my" matrix and array [1,3]
ysolve = np.linalg.solve(my, [1,3])
print(ysolve)
# verifies the solution is correct by multiplying the "my" matrix with "ysolve"
n = np.matmul(my, ysolve)
print(n)
# generates 10 normally distributed numbers
random = np.random.normal(0, 1, 10)
print(random)
# returns the mean of that random sample
mean = np.mean(random)
print(mean)
# returns the variance of that random sample
var = np.var(random)
print(var)
# returns the standard deviation of that random sample
std = np.std(random)
print(std)
# returns the median of that random sample
median = np.median(random)
print(median)
# returns quantile values in [0,.25,.5,.75,1]
quantile = pd.DataFrame(random).quantile([0, 0.25, 0.5, 0.75, 1])
print(quantile)
# returns range of values (max - min) in random sample
rang = np.ptp(random)
print(rang)
# returns the min value in random
minimum = np.min(random)
print(minimum)
# returns the max value in random
maximum = np.max(random)
print(maximum)
# yields the box plot of random values
fig = plt.figure(1, figsize=(8, 10))
ax = fig.add_subplot(111)
bp = ax.boxplot(random)
plt.show()
# statistical summary of the data sequence
random = np.random.normal(0, 1, 12)
summ = pd.DataFrame(random).describe()
print(summ)
# yields the histogram of 1000 random numbers with a normal distribution
w = np.random.normal(0, 1, 1000)
hist = plt.hist(w, facecolor='black')
plt.show()
# read the 5 deg gridded NOAAGlobalTemp .asc data
file = open('NOAAGlobalTemp.gridded.v4.0.1.201701.asc', "r")
# read float numbers from the file
da1 = []
for line in file:
x = line.strip().split()
for f in x:
da1.append(float(f))
len(da1)
da1[0 : 3]
tm1 = np.arange(0, 4267129, 2594)
tm2 = np.arange(1, 4267130, 2594)
len(tm1)
len(tm2)
# extract months
mm1 = []
for i in range(len(tm1)):
mm1.append(da1[tm1[i]])
# extract years
yy1 = []
for i in range(len(tm2)):
yy1.append(da1[tm2[i]])
# combine YYYY with MM
rw1 = []
for i in range(len(mm1)):
rw1.append(str(int(yy1[i])) + "-" + str(int(mm1[i])))
mm1[0 : 6]
yy1[0 : 6]
len(mm1)
len(yy1)
tm1[0 : 6]
tm2[0 : 6]
tm3 = pd.concat([pd.DataFrame(tm1), pd.DataFrame(tm2)], axis = 1)
# remove the months and years data from the scanned data
tm4 = []
for r in range(len(tm1)):
tm4.append(tm1[r])
tm4.append(tm2[r])
tm4.reverse()
da2 = da1.copy()
for f in tm4:
del da2[f]
len(da2)/(36*72) # months, 137 yrs 1 mon: Jan 1880-Jan 2017
# generates the space-time data
# 2592 (=36*72) rows and 1645 months (=137 yrs 1 mon)
var = [[0 for x in range(1645)] for y in range(2592)]
i = 0
for c in range(1645):
for r in range(2592):
var[r][c] = da2[i]
i += 1
da3 = pd.DataFrame(var, columns = rw1)
lat1=np.linspace(-87.5, 87.5, 36).tolist()
lon1=np.linspace(2.5, 357.5, 72).tolist()
Lat = sorted(lat1 * 72)
Lon = lon1 * 36
rw1.insert(0, "LAT")
rw1.insert(1, "LON")
gpcpst = pd.concat([pd.DataFrame(Lat), pd.DataFrame(Lon),
pd.DataFrame(da3)], axis = 1)
gpcpst.columns = rw1
gpcpst['2015-12'].values
gpcpst.columns
gpcpst.shape
# the first two columns are Lat and Lon
# -87.5 to 87.5 and then 2.5 to 375.5
# the first row for time is header, not counted as data.
gpcpst.to_csv("NOAAGlobalT1.csv") #Output the data as a csv file
# column '2015-12' corresponding to Dec 2015
# convert the vector into a lon-lat matrix for contour map plotting
# compresses numbers to [-6, 6] range
mapmat = np.reshape(gpcpst['2015-12'].values, (-1, 72)) #36 X 72
# this command compresses numbers to the range of -6 to 6
mapmat = np.maximum(np.minimum(mapmat, 6), -6)
dpi = 100
fig = plt.figure(figsize = (1100/dpi, 1100/dpi), dpi = dpi)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.9])
# create map
dmap = Basemap(projection = 'cyl', llcrnrlat = -87.5, urcrnrlat = 87.5,
resolution = 'c', llcrnrlon = 2.5, urcrnrlon = 357.5)
# draw coastlines, state and country boundaries, edge of map
dmap.drawcoastlines()
dmap.drawstates()
dmap.drawcountries()
# create and draw meridians and parallels grid lines
dmap.drawparallels(np.arange( -50, 100, 50.), labels = [1, 0, 0, 0],
fontsize = 10)
dmap.drawmeridians(np.arange(50, 350, 50.), labels = [0, 0, 0, 1],
fontsize = 10)
# convert latitude/longitude values to plot x/y values
# need to reverse sequence of lat1
#lat1_rev = lat1[::-1]
lat1_rev = lat1
x, y = dmap(*np.meshgrid(lon1, lat1_rev))
# contour levels
clevs = np.arange(-7, 7, 0.5)
print(x.shape,y.shape, mapmat.shape)
# draw filled contours
cnplot = dmap.contourf(x, y, mapmat, clevs, cmap = plt.cm.jet)
# add colorbar
# pad: distance between map and colorbar
cbar = dmap.colorbar(cnplot, location = 'right', pad = "2%")
cbar.set_label('[oC]') # add colorbar title string
# add plot title
plt.title('NOAAGlobalTemp Anomalies Dec 2015 [deg C]')
# label x and y
plt.xlabel('Longitude', labelpad = 20)
plt.ylabel('Latitude', labelpad = 40)
# display on screen
plt.show()
# keep only the data for the Pacific region for years from 1951-2000
# get data for region with -20 < Latitude < 20 and 160 < Longitude < 260
n2 = gpcpst[(gpcpst['LAT'] > -20) & (gpcpst['LAT'] < 20) &
(gpcpst['LON'] > 160) & (gpcpst['LON'] < 260)]
print(gpcpst.shape)
print(n2.size)
# from 1951-2000
# note Jan 1951 data starts from column 854 (note python is zero based)
pacificdat = n2.iloc[:, 854 : 1454]
Lat = np.arange(-17.5, 22.5, 5)
Lon = np.arange(162.5, 262.5, 5)
plt.close()
# get Dec 1997 data from pacificdat (tropical region)
mapmat = np.reshape(pacificdat['1997-12'].values, (8, -1)) # 8 * 20
fig = plt.figure(figsize=(1100/dpi, 1100/dpi), dpi = dpi)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.9])
# create map
smap = Basemap(projection = 'cyl', llcrnrlat = -40, urcrnrlat = 40,
resolution = 'c', llcrnrlon = 120, urcrnrlon = 300)
# draw coastlines, state and country boundaries, edge of map
smap.drawcoastlines()
smap.drawstates()
smap.drawcountries()
# create and draw meridians and parallels grid lines
smap.drawparallels(np.arange( -40, 60, 20.),labels=[1, 0, 0, 0],fontsize = 10)
smap.drawmeridians(np.arange(50, 350, 50.),labels=[0, 0, 0, 1],fontsize = 10)
# convert latitude/longitude values to plot x/y values
x, y = smap(*np.meshgrid(Lon, Lat))
# contour levels
clevs = np.arange(-4, 4, 0.1)
# draw filled contours
cnplot = smap.contourf(x, y, mapmat, clevs, cmap = plt.cm.jet)
# add colorbar
# pad: distance between map and colorbar
cbar = smap.colorbar(cnplot, location = 'right', pad = "2%")
# add colorbar title string
cbar.set_label('[oC]', verticalalignment = 'top')
# add plot title
plt.title('Tropic Pacific SAT Anomalies Dec 1997 [deg C]')
# label x and y
plt.xlabel('Longitude', labelpad = 20)
plt.ylabel('Latitude', labelpad = 40)
# display on screen
plt.show()