Python Code for Climate Science¶
In [1]:
####################################################################
#
#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
#
####################################################################
In [2]:
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
Python As a Smart Calculator¶
In [3]:
1+4 #The result 5 is shown as an Out[5] below.
Out[3]:
In [4]:
2 + np.pi/4 - 0.8 #\pi=3.14...... in Python is given by np.pi
Out[4]:
In [5]:
x=1
y=2
z=4
t=2*x**y-z #x**y means x^y: x to the power of y
t
Out[5]:
In [6]:
u=2
v=3
u+v
Out[6]:
In [7]:
np.sin(u*v) #np.sin is sine function in Python
Out[7]:
In [8]:
tmax=[77, 72,75,73,66,64,59] #Enter the temperature data
tmax #Show the temperarure data
Out[8]:
Define a Sequence in Python¶
In [9]:
np.arange(9) #A Python sequence starts from 0 while R from 1
#A sequence is called an array in Python
Out[9]:
In [10]:
# np.arange(x,y) creates an array of numbers incremented by 1
# between integer x and integer y-1
np.arange(1, 9)
Out[10]:
In [11]:
np.arange(-5.5, 2) #np.arange(x,y) can take negative and real numbers with increment 1
Out[11]:
In [12]:
# 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)
Out[12]:
In [13]:
# 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
Out[13]:
Define a Function in R¶
In [14]:
# Define the x^2 function and call it samfctn(x)
def samfctn(x):
return x * x;
samfctn(4)
Out[14]:
In [15]:
# Define a multivaraite function "x+y-z/2"
def fctn2(x, y, z):
return x + y - z / 2;
fctn2(1, 2, 3)
Out[15]:
Plot with Python¶
In [16]:
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
Out[16]:
Plot the curve of y=sin(x) from -pi to 2 pi¶
In [17]:
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()
Plot x-y data based on a formula¶
In [18]:
t = np.linspace(-3, 2, 20) #uniform division using 20 points
plt.plot(t, t ** 2)
plt.show()
In [19]:
np.size(t) #length of the data array
Out[19]:
In [20]:
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()
Plot a 3D surface, yields a 3D surface with a perspective angle¶
In [21]:
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()
In [22]:
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()
Contour plot¶
In [23]:
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()
In [24]:
## Symbolic calculations by Python
In [25]:
# 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)
Out[25]:
In [26]:
# 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)
In [27]:
# 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)
In [28]:
# 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)
In [29]:
# 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)
In [30]:
# Integrates x^2 from 0 to 1
sy.integrate(x ** 2, (x, 0, 1))
Out[30]:
In [31]:
# integrates cos(x) from 0 to pi/2
sy.integrate(sy.cos(x), (x, 0, math.pi / 2))
Out[31]:
Vectors, Matrices and Solve Linear Equations¶
In [32]:
# creates a 5X1 matrix
np.mat([[1], [6], [3], [np.pi], [-3]])
Out[32]:
In [33]:
# creates a sequence from 2 to 6
np.arange(2, 7)
Out[33]:
In [34]:
# creates a sequence incremented by 2 from 1 to 10
np.arange(1, 10, 2)
Out[34]:
In [35]:
# creates a 4X1 matrix
x = np.mat([[1], [-1], [1], [-1]])
print(x)
In [36]:
# adds number 1 to each element of x
print(x + 1)
In [37]:
# multiplies number 2 to each element of x
print(x * 2)
In [38]:
#x/2, divides each element of x by 2
print(x / 2)
In [39]:
# 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)
Out[39]:
In [40]:
print(x / 2)
print(x)
print(y)
In [41]:
# the dot product of two vectors
np.dot(y, x)
Out[41]:
In [42]:
x.shape
Out[42]:
In [43]:
# transforms x into a row 1X4 vector
x.T.shape
Out[43]:
In [44]:
# dot product and results in 1X1 matrix
np.dot(x.T, y.T)
Out[44]:
In [45]:
# this column times row yields a 4X4 matrix
np.matmul(x, y)
Out[45]:
In [46]:
# 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)
In [47]:
# finds dimensions of a matrix
my.shape
Out[47]:
In [48]:
# converts a matrix to a vector, again via columns
my.flatten()
Out[48]:
In [49]:
# reshapes x to a 2X2 matrix
mx = x.reshape(2, -1).transpose()
# multiplication between each pair of elements
np.multiply(mx, my)
Out[49]:
In [50]:
# division between each pair of elements
np.divide(mx, my)
Out[50]:
In [51]:
np.subtract(mx, 2 * my)
Out[51]:
In [52]:
# this is matrix multiplication in linear algebra
np.matmul(mx, my)
Out[52]:
In [53]:
# determinant of a matrix
x = np.linalg.det(my)
print(x)
In [54]:
# inverse of a matrix
z = np.linalg.inv(my)
print(z)
In [55]:
# 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)
In [56]:
# returns a left to right diagonal vector of a matrix
np.diagonal(my)
Out[56]:
In [57]:
# returns the eigenvalues of the "my" matrix
myeigval = np.linalg.eigvals(my)
print(myeigval)
In [58]:
# returns the eigenvalues and eigenvectors of the "my" matrix
myeigvect = np.linalg.eig(my)
print(myeigvect)
In [59]:
# returns the singular value decomposition(svd) of the "my" matrix
mysvd = np.linalg.svd(my)
print(mysvd)
In [60]:
# returns a solved linear matrix equation of the "my" matrix and array [1,3]
ysolve = np.linalg.solve(my, [1,3])
print(ysolve)
In [61]:
# verifies the solution is correct by multiplying the "my" matrix with "ysolve"
n = np.matmul(my, ysolve)
print(n)
Simple Statistics by Python¶
In [62]:
# generates 10 normally distributed numbers
random = np.random.normal(0, 1, 10)
print(random)
In [63]:
# returns the mean of that random sample
mean = np.mean(random)
print(mean)
In [64]:
# returns the variance of that random sample
var = np.var(random)
print(var)
In [65]:
# returns the standard deviation of that random sample
std = np.std(random)
print(std)
In [66]:
# returns the median of that random sample
median = np.median(random)
print(median)
In [67]:
# returns quantile values in [0,.25,.5,.75,1]
quantile = pd.DataFrame(random).quantile([0, 0.25, 0.5, 0.75, 1])
print(quantile)
In [68]:
# returns range of values (max - min) in random sample
rang = np.ptp(random)
print(rang)
In [69]:
# returns the min value in random
minimum = np.min(random)
print(minimum)
In [70]:
# returns the max value in random
maximum = np.max(random)
print(maximum)
In [71]:
# 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()
In [72]:
# statistical summary of the data sequence
random = np.random.normal(0, 1, 12)
summ = pd.DataFrame(random).describe()
print(summ)
In [73]:
# 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 data from files, process and plot¶
In [74]:
# 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)
Out[74]:
In [75]:
da1[0 : 3]
Out[75]:
In [76]:
tm1 = np.arange(0, 4267129, 2594)
tm2 = np.arange(1, 4267130, 2594)
len(tm1)
Out[76]:
In [77]:
len(tm2)
Out[77]:
In [78]:
# 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])))
In [79]:
mm1[0 : 6]
Out[79]:
In [80]:
yy1[0 : 6]
Out[80]:
In [81]:
len(mm1)
Out[81]:
In [82]:
len(yy1)
Out[82]:
In [83]:
tm1[0 : 6]
Out[83]:
In [84]:
tm2[0 : 6]
Out[84]:
In [85]:
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
Out[85]:
In [86]:
# 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
In [87]:
gpcpst['2015-12'].values
Out[87]:
In [88]:
gpcpst.columns
Out[88]:
In [89]:
gpcpst.shape
Out[89]:
In [90]:
# 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
Plot the climate data map for a given time¶
In [91]:
# 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()
Extract the data for a specific region¶
In [92]:
# 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]
Plot a climate map for a specific region¶
In [93]:
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()
Extract and plot data for a specified box with given lat and lon¶
In [94]:
# get region with Latitude = 32.5 and Longitude = 242.5
n2 = gpcpst[(gpcpst['LAT'] == 32.5) & (gpcpst['LON'] == 242.5)]
SanDiegoData = n2.iloc[:, 854 : 1454] #from 1951 - 2000
# note that using plt.subplots below is equivalent to using
fig, ax = plt.subplots()
ax.plot(np.linspace(1951, 2000,
len(SanDiegoData.values.T)), SanDiegoData.values[0, :])
ax.set(xlabel = 'Year', ylabel = 'Temp [oC]',
title = 'San Diego temperature anomalies history 1951-2000')
ax.grid()
plt.show()
Fill the few missing data¶
In [95]:
temp = gpcpst.copy()
areaw = np.zeros((2592, 1647))
areaw.shape
areaw[:, 0] = temp['LAT'].values
areaw[:, 1] = temp['LON'].values
# convert degrees to radians
veca=np.cos(temp['LAT'].values * np.pi / 180).tolist()
tempvalues = temp.values # get the underlying ndarray from the temp DataFrame
# area-weight matrix equal to cosine lat of the boxes with data
# and to zero for the boxes of missing data -999.9
for j in range(2, 1647):
for i in range(0, 2592):
if tempvalues[i, j] > -290.0:
areaw[i][j] = veca[i]
Compute an area-weighted temperature data matrix and its average¶
In [96]:
# area-weight data matrix first two columns as lat-lon
# create monthly global average vector for 1645 months (Jan 1880- Jan 2017)
tempw = np.multiply(areaw, tempvalues)
tempw[:, 0] = temp['LAT'].values
tempw[:, 1] = temp['LON'].values
avev= np.divide(tempw.sum(axis = 0)[2 : 1647], areaw.sum(axis = 0)[2 : 1647])
Plot time series of the global ave temp¶
In [97]:
#from sklearn.linear_model import LinearRegression
plt.close()
timemo = np.linspace(1880, 2017, 1645)
plt.plot(timemo, avev, '-')
plt.title('Area-weighted global average of monthly SAT anomalies: Jan 1880-Jan 2017')
plt.xlabel('Year')
plt.ylabel('Temperature anomaly [oC]')
# plot linear regression
lm = LinearRegression()
lm.fit(timemo.reshape(-1, 1), avev)
predictions = lm.predict(timemo.reshape(-1, 1))
plt.plot(timemo, predictions, '-')
font = {'family': 'serif',
'color': 'blue',
'weight': 'normal',
'size': 12,
}
plt.text(1885, .5, 'Linear trend: 0.69 [oC] per century', fontdict = font)
plt.show()
Plot the data coverage time series¶
In [98]:
# plot this time series
motime = np.linspace(1880, 2017, 1645)
rcover = 100 * areaw.sum(axis = 0)[2 : 1647] / sum(veca)
plt.close()
plt.plot(motime, rcover)
plt.ylim(0, 100)
plt.title('NOAAGlobalTemp Data Coverage: Jan 1880-Jan 2017')
plt.xlabel('Year')
plt.ylabel('Percent area covered [%]')
plt.show()
Work with the NOAAGlobalTemp time series data¶
In [99]:
# download the NCEI spatial average time series of monthly data
# https://www1.ncdc.noaa.gov/pub/data/noaaglobaltemp/operational/
# timeseries/aravg.mon.land_ocean.90S.90N.v4.0.1.201702.asc
aveNCEI = pd.read_table('http://shen.sdsu.edu/data/' +
'aravg.mon.land_ocean.90S.90N.v4.0.1.201703.asc.txt',
header = None, delim_whitespace = True)
aveNCEI.shape
aveNCEI = aveNCEI.values
avediff = np.subtract(avev, aveNCEI[0:1645, 2])
plt.close()
plt.plot(timemo, avediff, color = 'black')
plt.title('Difference of R average minus NCEI average of global temp')
plt.xlabel('Year')
plt.ylabel('Diffences [oC]')
plt.show()
Plot the annual time series of the NOAAGlobalTemp¶
In [100]:
# compute annual average, then plot the annual mean global average temp
avem = np.reshape(avev[0 : 1644], (-1, 12))
avem.shape #size 137 * 12, with row for the year, column for the month
# average over columns (annual average over 12 months),
# there are 137 years, 1880/01 - 2016/12
annv = np.mean(avem, axis = 1)
# plot the annual mean global average temp (1880-1996)
timeyr = np.arange(1880, 2017)
plt.close()
plt.plot(timeyr, annv, '-')
plt.title('Area-weighted global average of annual SAT anomalies: 1880-2016')
plt.xlabel('Year')
plt.ylabel('Temperature anomaly [oC]')
# plot linear regression
lm = LinearRegression()
lm.fit(timeyr.reshape(-1, 1), annv)
predictions = lm.predict(timeyr.reshape(-1 , 1))
plt.plot(timeyr, predictions, '-')
font = {'family': 'serif',
'color': 'blue',
'weight': 'normal',
'size': 12,
}
plt.text(1900, .4, 'Linear trend: 0.69 [oC] per century', fontdict = font)
font = {'family': 'serif',
'color': 'red',
'weight': 'normal',
'size': 12,
}
plt.text(1900, 0.05, 'Base line', fontdict = font)
plt.plot(timeyr, np.zeros(137), 'r-')
plt.show()
In [101]:
# one can compare with the NOAAGlobalTemp annual time series
# there are some small differences
aveannNCEI = pd.read_table('http://shen.sdsu.edu/data/' +
'aravg.ann.land_ocean.90S.90N.v4.0.1.201703.asc.txt',
header = None, delim_whitespace = True)
aveannNCEI.shape
aveannNCEI = aveannNCEI.values
diff2 = np.subtract(annv, aveannNCEI[0 : 137, 1])
# find the value range (minimum and maximum values in diff2 array)
[np.amin(diff2, axis = 0), np.amax(diff2, axis = 0)]
Out[101]:
Nonlinear trend from polynomial fitting¶
In [102]:
# polynomial fitting to the global average annual mean
# the following polynomial fitting from numpy.polyfit is NOT orthogonal,
# we know orthogonal polynomial fitting is usually better as we saw in the R code implemention
polyor9 = np.polyfit(timeyr, annv, 9)
polyor20 = np.polyfit(timeyr, annv, 20)
plt.close()
plt.plot(timeyr, annv, '-')
plt.title('Annual SAT time series and its orthogonal polynomial fits: 1880-2016')
plt.xlabel('Year')
plt.ylabel('Temperature anomaly [oC]')
plt.plot(timeyr, np.polyval(polyor9, timeyr), "-b",
label = "9th order orthogonal polynomial fit")
plt.plot(timeyr, np.polyval(polyor20, timeyr), "-r",
label = "20th order orthogonal polynomial fit")
plt.legend(loc = 'upper left')
plt.show()
Non-parametric smooth by LOWESS fitting¶
In [103]:
#import statsmodels.api as sm
lowess = sm.nonparametric.lowess
lowessfit = lowess(annv, timeyr, frac = 2./18)
plt.scatter(timeyr, annv, facecolors = 'none', edgecolors = 'b')
plt.plot(lowessfit[:, 0], lowessfit[:, 1])
plt.title('Annual SAT time series with non-parametric smooth lowess fitting: 1880-2016')
plt.xlabel('Year')
plt.ylabel('Temperature anomaly [oC]')
plt.show()
Compute and plot the trend of gridded data¶
In [104]:
# Compute the trend for each box for the 20th century
timemo1 = np.linspace(1900, 2000, 1200)
temp1 = temp.copy()
# replace missing values (any value less than -490.00) with NaN in temp1
# dataframe
temp1 = temp1.where(temp1 >= -490.0, np.NaN)
trendgl = np.zeros(2592)
# for grids in (Lon: 2.5-257.5, Lat: -87.5-87.5), perform linear regression
# fitting on each grid temperature over months between Jan 1900 and Dec 1999
# linear slope cofficient will be saved in trendgl list to be used for
# subsequent plotting
for i in range(0, 2592):
if not(np.isnan(temp1.iloc[i, 242]) or np.isnan(temp1.iloc[i, 1441])):
# column 242 corresponds to month 1990-01,
# and column 1441 corresponds to month 1999-12
# temperature data for a grid from 1990-01 to 1999-12
gridtemp = temp1.iloc[i, 242 : 1442].values
# filter out any month where grid temperature is missing
# (i.e., having a NaN value)
idx = np.isfinite(timemo1) & np.isfinite(gridtemp)
# check if there is at least one month that has
# non-NaN number for linear regression
if np.any(idx[:]):
ab = np.polyfit(timemo1[idx], gridtemp[idx], 1)
# ab[0] has the linear regression slope coefficient
trendgl[i] = ab[0]
else:
trendgl[i] = np.NaN
else:
trendgl[i] = np.NaN
mapmat = np.reshape([n * 100 for n in trendgl], (-1, 72)) #36 X 72
# this command compresses numbers to the range of -2 to 2
mapmat = np.maximum(np.minimum(mapmat, 2), -2)
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
amap = 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
amap.drawcoastlines()
amap.drawstates()
amap.drawcountries()
# create and draw meridians and parallels grid lines
amap.drawparallels(np.arange( -50, 100, 50.),
labels = [1, 0, 0, 0], fontsize = 10)
amap.drawmeridians(np.arange(50, 350, 50.),
labels = [0, 0, 0, 1], fontsize = 10)
# convert latitude/longitude values to plot x/y values
x, y = amap(*np.meshgrid(lon1, lat1))
# contour levels
clevs = np.linspace(-2, 2, 21)
# draw filled contours
cnplot = amap.contourf(x, y, mapmat, clevs, cmap = plt.cm.jet)
# add colorbar
# pad: distance between map and colorbar
cbar = amap.colorbar(cnplot, location = 'right', pad = "2%")
cbar.set_label('[oC]') # add colorbar title string
cbar.set_ticks(np.arange(-2, 2.6, 0.5))
cbar.set_ticklabels(np.arange(-2, 2.6, 0.5))
# add plot title
plt.title('Jan 1900-Dec 1999 temperature trends: [oC/century]')
# label x and y
plt.xlabel('Longitude', labelpad = 20)
plt.ylabel('Latitude', labelpad = 40)
plt.show()
Compute trend for each box for the 20th century - Version 2: Allow 2/3 of data, i.e., 1/3 missing¶
In [105]:
# compute the trend and plot the 20C V2 trend map
timemo1 = np.linspace(1900, 2000, 1200)
temp1 = temp.iloc[:, 242 : 1442]
# replace missing values (any value less than -490.00)
# with NaN in temp1 dataframe
temp1 = temp1.where(temp1 >= -490.0, np.NaN)
temptf = np.isnan(temp1)
countOfNotMissing = np.zeros(2592)
for i in range(0, 2592):
countOfNotMissing[i] = temptf.iloc[i, :].values.tolist().count(False)
trend20c = np.zeros(2592)
for i in range(0, 2592):
if countOfNotMissing[i] > 800: # allow 2/3 of data, i.e., 1/3 missing
# column 242 corresponds to month 1990-01,
# and column 1441 corresponds to month 1999-12
gridtemp = temp1.iloc[i, :].values # temperature for a grid
# from 1990-01 to 1999-12
# filter out any month where grid temperature is missing
# (i.e., having a NaN value)
idx = np.isfinite(timemo1) & np.isfinite(gridtemp)
ab = np.polyfit(timemo1[idx], gridtemp[idx], 1)
# ab[0] has the linear regression slope coefficient
trend20c[i] = ab[0]
else:
trend20c[i] = np.NaN
mapmat = np.reshape([n * 120 for n in trend20c], (-1, 72)) #36 X 72
# this command compresses numbers to the range of -3 to 3
mapmat = np.maximum(np.minimum(mapmat, 3), -3)
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
fmap = 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
fmap.drawcoastlines()
fmap.drawstates()
fmap.drawcountries()
# create and draw meridians and parallels grid lines
fmap.drawparallels(np.arange( -50, 100, 50.),
labels = [1, 0, 0, 0], fontsize = 10)
fmap.drawmeridians(np.arange(50, 350, 50.),
labels = [0, 0, 0, 1], fontsize = 10)
# convert latitude/longitude values to plot x/y values
x, y = fmap(*np.meshgrid(lon1, lat1))
# contour levels
clevs = np.linspace(-3, 3, 41)
# draw filled contours
cnplot = fmap.contourf(x, y, mapmat, clevs, cmap = plt.cm.jet)
# add colorbar
# pad: distance between map and colorbar
cbar = fmap.colorbar(cnplot, location = 'right', pad = "2%")
cbar.set_label('[oC]') #Add colorbar title string
cbar.set_ticks(np.arange(-3, 3.1, 1))
cbar.set_ticklabels(np.arange(-3, 3.1, 1))
# add plot title
plt.title('Jan 1900-Dec 1999 temperature trends: [oC/decade]')
# label x and y
plt.xlabel('Longitude', labelpad = 20)
plt.ylabel('Latitude', labelpad = 40)
plt.show()
Trend from 1976-2016¶
In [106]:
timemo2 = np.linspace(1976, 2017, 492)
temp1 = temp.copy()
#-- replace missing values (any value < -490.00) with NaN in temp1 dataframe --
temp1 = temp1.where(temp1 >= -490.0, np.NaN)
trend7616 = np.zeros(2592)
# for grids in (Lon: 2.5-257.5, Lat: -87.5-87.5), perform linear regression
# fitting on each grid temperature over months between Jan 1976 and Dec 2016
# linear slope cofficient will be saved in trendgl list to be used for
# subsequent plotting
for i in range(0, 2592):
if not(np.isnan(temp1.iloc[i, 1154]) or np.isnan(temp1.iloc[i, 1645])):
# column 1154 corresponds to month 1976-01, and column 1645
# corresponds to month 2016-12
# temperature data for a grid from 1976-01 to 2016-12
gridtemp = temp1.iloc[i, 1154 : 1646].values
# filter out any month where grid temperature is missing
# (i.e., having a NaN value)
idx = np.isfinite(timemo2) & np.isfinite(gridtemp)
# check if there is at least one month that has non-NaN number
# for linear regression
if np.any(idx[:]):
ab = np.polyfit(timemo2[idx], gridtemp[idx], 1)
# ab[0] has the linear regression slope coefficient
trend7616[i] = ab[0]
else:
trend7616[i] = np.NaN
else:
trend7616[i] = np.NaN
Plot US temp and prec times series on the same figure¶
In [107]:
Time = np.arange(2001, 2011)
Tmean = [12.06, 11.78, 11.81, 11.72, 12.02, 12.36, 12.03, 11.27, 11.33, 11.66]
Prec = [737.11, 737.87, 774.95, 844.55, 764.03, 757.43, \
741.17, 793.50, 820.42, 796.80]
plt.close()
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_title('Contiguous U.S. Annual Mean Temperature and Total Precipitation')
ax1.set_xlabel('Year')
ax1.set_ylabel('Tmean [dec C]', color = color)
ax1.plot(Time, Tmean, color = color, marker = 'o', markerfacecolor = 'white')
ax1.tick_params(axis = 'y', labelcolor = color)
ax1.legend('Tmean', loc = 'upper left')
ax2 = ax1.twinx() # instantiate a second axis that shares the same x-axis
color = 'tab:blue'
# we already handled the x-label with ax1
ax2.set_ylabel('Precipitation [mm]', color = color)
ax2.plot(Time, Prec, color = color, marker = 'o', markerfacecolor = 'white')
ax2.tick_params(axis = 'y', labelcolor = color)
ax2.legend('Prec')
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
Margins, math symbol, and figure setups¶
In [108]:
x = pd.DataFrame(np.arange(-30, 31)) * 0.25
y = np.sin(x)
x1 = x[np.sin(x) >= 0].dropna().transpose().values[0]
y1 = np.sin(x1)
x2 = x[np.sin(x) < 0].dropna().transpose().values[0]
y2 = np.sin(x2)
fig,ax1 = plt.subplots()
bar1 = ax1.bar(x1, y1, width = 0.04, ec = "k")
for i in range(len(bar1)):
bar1[i].set_color('r')
ax1.set_ylim(-1, 1)
ax1.set_xticks(np.arange(-6, 6.1, step = 2))
ax1.set_yticks(np.arange(-1, 1.1, step = 0.5))
ax2 = ax1.twinx() # instantiate a second axis that shares the same x-axis
bar2 = ax2.bar(x2, y2, width = 0.1, linestyle = '--', ec = "k")
for i in range(len(bar2)):
bar2[i].set_color('b')
ax2.set_ylim(-1, 1)
ax2.set_yticks(np.arange(-1, 1.1, step = 0.5))
ax2.set_yticklabels(['A', 'B', 'C', 'D', 'E'], fontsize = 14, fontweight = 2)
font = {'family' : 'arial',
'color' : 'green',
'weight' : 'bold',
'size' : 36,
'style' : 'italic'
}
plt.text(-5, .5, 'Sine waves', fontdict = font)
plt.text(-10.6, 0, "y=sin($\\theta$-$\hat\phi$)", color = 'b', size = 14,
rotation = 90, verticalalignment = 'center')
plt.text(0.5, -1.3, "Angle in radian: $\\theta$-$\phi_0$", color = 'r',
size = 17, horizontalalignment = 'center')
plt.text(0.5, 1.14, "Text outside of the figure on side 3", weight = 'bold',
size = 14, horizontalalignment = 'center')
plt.show()
Change the font sizes for axis labels, main title, and sub-title¶
In [109]:
np.random.seed(1)
plt.scatter([n * (1/20) for n in np.arange(1, 201)], np.random.randn(200),
facecolors = 'none', edgecolors = 'black')
plt.title('Normal random values', fontsize = 20, weight = 'bold')
plt.xlabel('Time', fontsize = 12)
plt.ylabel('Random values', fontsize = 12)
plt.text(5, -4.2, "Sub-title: 200 random values", size = 16,
horizontalalignment = 'center')
plt.show()
Wind directions due to the balance between PGF and Coriolis force using an arrow plot for vector fields on a map¶
In [110]:
NOAATemp = np.loadtxt("aravg.ann.land_ocean.90S.90N.v4.0.1.201805.asc.txt",skiprows=0)
print(NOAATemp)
x=NOAATemp[:,0]
y=NOAATemp[:,1]
z = [-99] * x.size
print(y.shape)
In [111]:
def forz(z,y):
for i in range(2,137):
rslt = [(y[i-2],y[i-1],y[i],y[i+1],y[i+2])]
z[i] = np.mean(rslt)
return z
In [112]:
zz = forz(z,y)
print(zz)
In [113]:
def n1func():
n1 = []
for i in range(0,139):
if y[i] >= 0:
n1.append(i)
return n1
def n2func():
n2 = []
for i in range(0,139):
if y[i] < 0:
n2.append(i)
return n2
In [114]:
n1 = n1func()
n2 = n2func()
print(n1,n2)
In [115]:
x1= x[n1]
y1= y[n1]
x2 = x[n2]
y2 = y[n2]
x3 = x[2:136]
y3 = z[2:136]
In [116]:
plt.bar(x1,y1,color='r',width=1.)
plt.bar(x2,y2,color='b',width=1.)
plt.plot(x3,y3,color='k')
plt.xlim(1877,2020)
plt.xlabel('Time')
plt.ylabel('Temperature $^{o}$C')
plt.title('NOAA Global Average Annual Mean Temperature Anomalies')
Out[116]:
In [117]:
#Plot US temp and prec times series on the same figure
#lines 598-610
Time = [2001,2002,2003,2004,2005,2006,2007,2008,2009,2010]
Tmean = [12.06, 11.78,11.81,11.72,12.02,12.36,12.03,11.27,11.33,11.66]
Prec = [737.11,737.87,774.95,844.55,764.03,757.43,741.17,793.50,820.42,796.80]
plt.subplot(2,1,1)
plt.plot(Time,Tmean,color="r")
plt.ylabel("TMean ($^{o}$C)")
plt.title("US Annual Mean Temperature")
plt.subplot(2,1,2)
plt.plot(Time,Prec,color="b")
plt.ylabel("Prec (mm)")
plt.title("US Annual Total Precipitation")
plt.subplots_adjust(top=1.1, bottom=0.1, left=0.10, right=0.95, hspace=0.25,
wspace=0.35)
In [118]:
#Stack multiple panels together as a single figure
#lines 613-618
x1 = np.linspace(0.,20.,100)
x2 = np.linspace(0.,10.,100)
x3 = np.linspace(10.,20.,100)
y1 = np.sin(x1)
y2 = np.sin(x2)
y3 = np.sin(x3)
plt.subplot(221)
plt.plot(x1,y1)
plt.xlabel("x")
plt.ylabel("sin")
plt.subplot(222)
plt.plot(x2,y1)
plt.xlabel("x")
plt.ylabel("sin")
plt.subplot(223)
plt.plot(x3,y1)
plt.xlabel("x")
plt.ylabel("sin")
plt.subplots_adjust(top=1.05, bottom=0.07, left=0.10, right=0.95, hspace=0.25,
wspace=0.35)
In [119]:
a = np.linspace(-90.,90.,37)
b = np.linspace(0.,359.9,72)
mapdat = np.reshape(np.random.normal(loc=0.,size=72*37),(37,72))
In [120]:
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAn = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data because
lons, data = ReAn.shiftdata(b, datain = mapdat, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAn.drawcoastlines(color='black', linewidth=1)
ReAn.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
ReAn.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
ReAn.drawcountries()
limit=np.linspace(-3.,3.,81)
ReAn_plt = plt.contourf(np.array(lons),np.array(a),data,limit,cmap = 'jet')
ceb = plt.colorbar(ReAn_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 30)
plt.ylabel('Latitude',labelpad =30)
plt.title('Standard Normal Random Values on a World Map: 5 Lat-Lon Grid')
plt.show()
In [121]:
#Plot a regional map of a random dataset on a map
#Plot a 5-by-5 grid regional map to cover USA and Canada
#lines 646-657
c = np.linspace(10.,70.,13)
d = np.linspace(230.,295.,14)
mapdat2 = np.reshape(np.random.normal(loc=0.,size=14*13),(13,14))
In [122]:
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAn = Basemap(projection='cyl',llcrnrlat=10.,urcrnrlat=70.,\
llcrnrlon=-130.,urcrnrlon=-65.,resolution='c')
#shifting the data because
lons2, data = ReAn.shiftdata(d, datain = mapdat2, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAn.drawcoastlines(color='black', linewidth=1)
ReAn.drawparallels(np.arange(-90.,91.,10.), labels = [1,0,0,1])
ReAn.drawmeridians(np.arange(-180.,181.,10.), labels= [1,0,0,1])
ReAn.drawcountries()
limit=np.linspace(-3.,3.,81)
ReAn_plt = plt.contourf(np.array(lons2),np.array(c),data,limit,cmap = 'jet')
ceb = plt.colorbar(ReAn_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =30)
plt.title('Standard Normal Random Values on a World Map to cover USA and Canada: 5 Lat-Lon Grid')
plt.show()
In [123]:
#R plot of NCEP/NCAR Reanalysis PSD monthly temp data .nc file
#http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.surface.html
#lines658-698
nc = ds("air.mon.mean.nc","r+")
In [124]:
#define the variables from the data set
lon_vals = nc.variables['lon']
lat_vals = nc.variables['lat']
time = nc.variables['time']
air = nc.variables['air']
time_unit = time.units
In [125]:
precnc = nc.variables['air']
nc_x = np.linspace(-90,90,73)
plt.plot(nc_x,precnc[0,:,71],ls="-",color='k' )
plt.xlabel('Latitude')
plt.ylabel('Temp[$^{\circ}$C]')
plt.title('90S-90N Temperature [$^{o}$C] Along a Meridional Line at 35E: Jan 1948')
plt.show()
In [126]:
def spmat(x,y):
for i in range(0,845):
y[:,i] = np.reshape(x[i,:,:],(144*73))
return y
In [127]:
precst = np.zeros((10512,845))
temp = np.reshape(precnc[0,:,:],(144*73))
precst2 = spmat(precnc,precst)
In [128]:
#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
In [129]:
arra = np.zeros((0))
LAT = rep(arra,lat_vals,144)
LON = rep(arra,lon_vals[0],73)
In [130]:
def recon():
rslt = LON
for jj in range(0,143):
rslt = np.append(rslt,rep(arra,lon_vals[jj],73))
return rslt
In [131]:
LON2 = recon()
In [132]:
#With 826 months spanning from 1948-2016
gpcpst_ = np.column_stack((LAT,(LON2)))
gpcpst = np.column_stack((gpcpst_,precst2))
#Creating an array with first two columns of lat and lon
ones = np.ones(845)
In [133]:
#define a function to create a list of months repeating to match
#the size of our datadef monthfunc():
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()
In [134]:
#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()
In [135]:
#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,845):
empty = np.append(empty,str(years[i]) + '-' + str(months[i]))
return empty
tm1 = tmfunc()
In [136]:
#add the 'Lat' and 'Lon' columns to the beginning
#of our column names array
tm2 = np.append(['Lat','Lon'],tm1)
In [137]:
#use the package pandas.DataFrame to match the column names
#with the data
GPCPST = pd.DataFrame(gpcpst,columns=tm2)
In [138]:
#save the data as a .csv file
GPCPST.to_csv('ncepJan1948_May2018.csv')
In [139]:
monJ = np.array(range(0,844,12))
gpcpdat = gpcpst[:,2:844]
gpcpJ = gpcpdat[:,monJ]
climJ = gpcpJ.mean(axis=1)
In [140]:
def rowSDS():
rslt = []
for jj in range(0,10512):
rslt = np.append(rslt,np.std(gpcpJ[jj,:]))
return rslt
sdJ = rowSDS()
In [141]:
def repeat():
q = np.zeros((144,73))
u = climJ[0:73]
for ii in range(0,144):
q[ii,:] = u
return q
In [142]:
Lat1 = np.linspace(90.,-90.,73)
Lat = -Lat1
mapmat = np.reshape(climJ,(73,144))
In [143]:
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAn = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data because
lons, data = ReAn.shiftdata(lon_vals, datain = mapmat, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAn.drawcoastlines(color='black', linewidth=1)
ReAn.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
ReAn.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
ReAn.drawcountries()
limit = np.linspace(-50.,50.,81)
ReAn_plt = plt.contourf(np.array(lons),np.array(lat_vals),data,limit,cmap = 'jet')
ceb = plt.colorbar(ReAn_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('NCEP Jan SAT RA 1948-2018 climatology [$^{o}$C]')
Out[143]:
In [144]:
#Compute standard deviation Jan 1948-Dec 2016
#lines730-744
Lat = np.linspace(90.,-90.,73)
Lat = -Lat
mapmat2 = np.reshape(sdJ,(73,144))
In [145]:
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAnSd = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data because
lons, data = ReAnSd.shiftdata(lon_vals, datain = mapmat2, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAnSd.drawcoastlines(color='black', linewidth=1)
ReAnSd.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
ReAnSd.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
ReAnSd.drawcountries()
limit = np.linspace(0.,10,40)
ReAnSd_plt = plt.contourf(np.array(lons),np.array(lat_vals),data,limit,cmap = 'jet')
ceb = plt.colorbar(ReAnSd_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('NCEP Jan SAT RA 1948-2018 Standard Deviation [$^{o}$C]')
plt.show()
In [146]:
#Plot the January 1983 temperature using the above setup
#lines756-770
Lat = np.linspace(90.,-90.,73)
Lat = -Lat
mapmat83J = np.reshape(precnc[421,:,:],(73,144))
In [147]:
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAnSd2 = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data
lons, data3 = ReAnSd2.shiftdata(lon_vals, datain = mapmat83J, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAnSd2.drawcoastlines(color='black', linewidth=1)
ReAnSd2.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
ReAnSd2.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
ReAnSd2.drawcountries()
limit = np.linspace(-50.,50,80)
ReAnSd2_plt = plt.contourf(np.array(lons),np.array(lat_vals),data3,limit,cmap = 'jet')
ceb = plt.colorbar(ReAnSd2_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('January 1983 Surface Air Temperature [$^{o}$C]')
plt.show()
In [148]:
#Plot the January 1983 temperature anomaly from NCEP data
#lines771-782
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAnSd3 = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data because
lons, data4 = ReAnSd3.shiftdata(lon_vals, datain = mapmat83J, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAnSd3.drawcoastlines(color='black', linewidth=1)
ReAnSd3.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
ReAnSd3.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
ReAnSd3.drawcountries()
limit = np.linspace(-50.,50.,81)
ReAnSd3_plt = plt.contourf(np.array(lons),np.array(lat_vals),data4,limit,cmap = 'jet')
ceb = plt.colorbar(ReAnSd3_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('January 1983 Surface Air Temperature Anomaly [$^{o}$C]')
plt.show()
In [149]:
#line 793 -- Creating arrays of longitudes and latituds
a = np.array([[-75, -45, -15, 15, 45, 75]])
b = np.array([[-165,-135,-105,-75,-45,-15,15,45,75,105,135,165]])
#concatenating, or combining each array so that they are the same size
lat = np.concatenate((a,a,a,a,a,a,a,a,a,a,a,a), axis = 0)
lon = np.concatenate((b,b,b,b,b,b), axis = 0)
#transposing the newly created arrays
lat_tran = np.transpose(lat)
lon_tran = np.transpose(lon)
In [150]:
#vectors u and v line 799 -- Creating vectors for wind directions
c = np.array([[-1,1,-1,-1,1,-1]])
d = np.array([[1,-1,1,-1,1,-1]])
#concatenating in a similar way as before to match sizes
u = np.concatenate((c,c,c,c,c,c,c,c,c,c,c,c), axis = 0)
v = np.concatenate((d,d,d,d,d,d,d,d,d,d,d,d), axis = 0)
In [151]:
#choosing an arbitrary size for the image
plt.figure(figsize=(10., 5.))
#using the basemap package to put the cylindrical projection map into the image
wmap = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#draw coastlines, latitudes, and longitutes on the map
wmap.drawcoastlines(color='black', linewidth=1)
wmap.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
wmap.drawmeridians(np.arange(-180.,181.,30.), labels= [1,0,0,1])
#plot points in the center of each box made by latitude and longitude
x, y = lon_tran, lat
wmap.scatter(x,y,marker='D',color='k')
#plot vector arrows originating from the midpoints
wmap.quiver(x, y, u, v, scale=30, color = 'b' )
#add axis titles
plt.text(-95,-4, 'Intertropical Convergence Zone (ITCZ)', color='r', fontsize = 15)
plt.text(15,-35, 'Subtropical High', color = 'r', fontsize = 15)
plt.text(15,27, 'Subtropical High', color = 'r', fontsize = 15)
plt.text(-20,92, 'Polar High', color = 'r', fontsize = 15)
#plot
plt.show()
In [152]:
#The following is heavily based on the following stackoverflow link:
#https://stackoverflow.com/questions/39742305/how-to-use-basemap
#-python-to-plot-us-with-50-states
#Lambert Conformal map of lower 48 states.
ggplot = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
#draw state boundaries.
#data from U.S Census Bureau
#http://www.census.gov/geo/www/cob/st2000.html
shp_info = ggplot.readshapefile('st99_d00','states',drawbounds=True)
#create a popdensity array to fill with random integers
popdensity = {
'New Jersey':0,
'Rhode Island':0,
'Massachusetts':0,
'Connecticut':0,
'Maryland':0,
'New York':0,
'Delaware':0,
'Florida':0,
'Ohio':0,
'Pennsylvania':0,
'Illinois':0,
'California':0,
'Hawaii':0,
'Virginia':0,
'Michigan':0,
'Indiana':0,
'North Carolina':0,
'Georgia':0,
'Tennessee':0,
'New Hampshire':0,
'South Carolina':0,
'Louisiana':0,
'Kentucky':0,
'Wisconsin':0,
'Washington':0,
'Alabama':0,
'Missouri':0,
'Texas':0,
'West Virginia':0,
'Vermont':0,
'Minnesota':0,
'Mississippi':0,
'Iowa':0,
'Arkansas':0,
'Oklahoma':0,
'Arizona':0,
'Colorado':0,
'Maine':0,
'Oregon':0,
'Kansas':0,
'Utah':0,
'Nebraska':0,
'Nevada':0,
'Idaho':0,
'New Mexico':0,
'South Dakota':0,
'North Dakota':0,
'Montana':0,
'Wyoming':0,
'Alaska':0}
colors={}
statenames=[]
cmap = plt.cm.jet # use 'jet' colormap
vmin = 0; vmax = 450 # set range.
#fill popdensity with random integers
for shapedict in ggplot.states_info:
statename = shapedict['NAME']
# skip DC and Puerto Rico.
if statename not in ['District of Columbia','Puerto Rico']:
popdensity[statename] = randint(0,100)
for shapedict in ggplot.states_info:
statename = shapedict['NAME']
# skip DC and Puerto Rico.
if statename not in ['District of Columbia','Puerto Rico']:
pop = popdensity[statename]
# calling colormap with value between 0 and 1 returns
# rgba value. Invert color range (hot colors are high
# population), take sqrt root to spread out colors more.
colors[statename] = cmap(1.-np.sqrt((pop-vmin)/(vmax-vmin)))[:3]
statenames.append(statename)
#cycle through state names, color each one.
ax = plt.gca() # get current axes instance
for nshape,seg in enumerate(ggplot.states):
# skip DC and Puerto Rico.
if statenames[nshape] not in ['District of Columbia','Puerto Rico']:
color = rgb2hex(colors[statenames[nshape]])
poly = Polygon(seg,facecolor=color,edgecolor=color)
ax.add_patch(poly)
plt.title('Filling State Polygons with Random Integer data')
plt.show()
Plot the wind field over the ocean¶
In [153]:
#downloading data from said url
downloaded_data = request.urlopen('ftp://eclipse.ncdc.noaa.gov/pub/seawinds/SI/uv/clm/uvclm95to05.nc','Users/louis/AnacondaProjects/uvclm95to05.nc')
#using netCDF4.Dataset to open the data
mincwind = ds("uvclm95to05.nc",'r+',"NETCDF4")
In [154]:
#print the information from the data along with its array shape
#shape = null
print(np.shape(mincwind))
print(mincwind)
In [155]:
#defining and extracting 'u' and 'v' variable
u_0 = mincwind.variables['u']
v_0 = mincwind.variables['v']
#reshaping the arrays into 3D
u = np.reshape(u_0,(12,719,1440))
v = np.reshape(v_0,(12,719,1440))
print(np.shape(v))
print(np.shape(u))
In [156]:
#basic filled contour plot
plt.figure(figsize=(12,7))
levels =range(-23,23)
ucont = plt.contourf(u[8,:,:],levels)
cbar = plt.colorbar(ucont)
plt.savefig('Figure2.png', bbox_inches='tight')
plt.show()
In [157]:
#a preset color map titled 'hot'
plt.figure(figsize=(12,7))
ucont2 = plt.contourf(u[8,:,:],levels,cmap = 'hot')
cbar = plt.colorbar(ucont2)
plt.show()
In [158]:
#a preset colormap titled 'Red to Blue'
plt.figure(figsize=(12,7))
ucont3 = plt.contourf(u[8,:,:],levels,cmap = 'RdBu')
cbar = plt.colorbar(ucont3)
plt.show()
In [159]:
#basic contour plot without filling
plt.figure(figsize=(12,7))
ucont4 = plt.contour(u[8,:,:],levels)
cbar = plt.colorbar(ucont4)
plt.show()
Generate 3D data¶
In [160]:
#create two arrays, x and y, ranging from 0 to 2pi
#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))
In [161]:
#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(1,11,1):
for i in range(0,100):
for j in range(0,100):
mydat[i,j,t-1] = (np.sin(t)*(1./np.pi)*np.sin(x[i])*np.sin(y[j])
+ np.exp(-.3*t)*(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)
Plot the original z(x,y,t) waves for a given t.¶
In [162]:
#set a range for the contour plot
maxi = np.linspace(-.5,.5,100)
#define the contour plot using previous x, y, and mydata arrays
CS = plt.contourf(x,y,mydata[:,:,0],maxi,cmap='hsv')
#add titles for the graph and axes
plt.title('Original field at t=10')
plt.xlabel('x')
plt.ylabel('y')
#add a color bar with a label
cbar = plt.colorbar(CS)
cbar.ax.set_ylabel('Scale')
#save figure as a .png file
plt.savefig('Figure1.png',bbox_inches='tight')
3D array to 2D space-time matrix.¶
In [163]:
#create two arrays, x and y, ranging from 0 to 2pi
x = y = np.linspace(0.,2.*np.pi,100)
In [164]:
#create a matrix of zeros of size (10000,10)
da = np.matrix(np.zeros((np.size(x)*np.size(y),10)))
print(np.shape(da))
In [165]:
#define a function that reshapes mydata into a new array of size
#(10000,10)
def twod_func(x,y,a):
for i in range(0,10):
a[:,i] = np.reshape(mydata[:,:,i],(10000,1))
return a
In [166]:
da1 = twod_func(x,y,da)
In [167]:
#SVD for EOFs
u, s, vh = np.linalg.svd(da1,full_matrices=False)
v = np.transpose(vh)
In [168]:
#Dimensions of EOFs, PCs, and Energy for the corresponding EOFs
print('Energy for the corresponding EOFs:',s.shape)
print('PCs:',v.shape)
print('EOFs:',u.shape)
In [169]:
#to show the accuracy of the svd
x_a = np.dot(np.dot(u,np.diag(s)),vh)
print("Standard deviation of our original array:",np.std(da1))
print("Standard deviation of recreated array",np.std(x_a))
print("Standard deviation of original array minus recreated:",np.std(da1-x_a))
In [170]:
print((s**2.)/np.sum(s**2.))
print(np.shape(u))
In [171]:
#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
u_mat = dig()
print(u_mat[0,:])
In [172]:
#Graph the EOF
size = np.linspace(-2.5e-2,2.5e-2,25)
dig = plt.contourf(x,y,u_mat,size,cmap='hsv')
plt.title('SVD Mode 1: EOF1')
plt.xlabel('x')
plt.ylabel('y')
cbar = plt.colorbar(dig)
cbar.ax.set_ylabel('Scale')
plt.figure()
Out[172]:
Accurate spatial patterns from functions that generate data.¶
In [173]:
#Create an array of all zeros to fill with data
#line 961
x = y = np.linspace(0.,2.*np.pi,100)
zeros = np.zeros((100,100))
In [174]:
#define two different funcitons 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(5.*x)*np.sin(5.*y)
In [175]:
#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
In [176]:
#plot the data
lvl = np.linspace(-.35,.35,35)
contour = plt.contourf(x,y,fcn1(x,y),lvl,cmap='hsv')
plt.title('Accurate Mode 1')
plt.xlabel('x')
plt.ylabel('y')
cbar = plt.colorbar(contour)
cbar.ax.set_ylabel('Scale')
plt.figure()
Out[176]:
Plot PCs and coefficitents of the functional patterns¶
In [177]:
t = np.array(range(1,11))
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 $-sin(t)$: 91% variance')
plt.plot(t,-np.exp(-.3*t),color='m',marker='o',label='Original mode 2 coefficient $e^{-0.3t}$: 9% variance')
plt.ylim(-1.1,1.1)
plt.xlabel('Time')
plt.ylabel('PC of Coefficient')
plt.title('SVD PCs vs. Accurate Temporal Coefficients')
plt.legend(loc=(1.1,.5))
plt.show()
Orthonormal Verification¶
In [178]:
#Create arrays to compute a dot product
t = np.array(range(1,11))
def sinarr(t):
rslt = np.zeros(t)
for ii in range(1,t+1):
rslt[ii-1] = -np.sin(ii)
return rslt
def exparr(t):
rslt = np.zeros(t)
for jj in range(1,t+1):
rslt[jj-1] = -np.exp(-.3*jj)
return rslt
In [179]:
sin = sinarr(10)
exp = exparr(10)
print('The dot product is:',np.dot(np.transpose(sin),exp))
In [180]:
#Now to find the variances and compare them
v1 = np.var(sin)
v2 = np.var(exp)
print(v1/(v1+v2))
In [181]:
#Verify orthogonality of PCs
v_1 = np.transpose(v[:,0])
v_2 = v[:,1]
print("The dot product of PCs is:" , np.dot(v_1,v_2))
Data reconstruction from 2 modes¶
In [182]:
#create arrays that will be used below
#and check the sizes of the arrays for matrix multiplication
sdiag = np.diag(s)
print(np.shape(sdiag[0:1,0:1]))
print(np.shape(u[:,0:1]))
print(np.shape(v[:,0:1]))
vtran = np.transpose(v[:,0:1])
In [183]:
#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))
In [184]:
#reshaping B to plot it
BB = np.reshape(B[4,:],(100,100))
scale = np.linspace(-.33,.36,25)
BBplot = plt.contourf(x,y,BB,scale,cmap = 'hsv')
plt.title('2-mode SVD reduced field t = 5')
plt.xlabel('x')
plt.ylabel('y')
cbar = plt.colorbar(BBplot)
cbar.ax.set_ylabel('Scale')
plt.figure()
Out[184]:
4 dimensions: lon, lat, level, time¶
In [185]:
#first, 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
#then use the Dataset function, ds, to use the data
nc = ds("air.mon.mean.nc","r+")
In [186]:
print(nc)
In [187]:
#define the variables from the data set
lon_vals = nc.variables['lon']
lat_vals = nc.variables['lat']
time = nc.variables['time']
air = nc.variables['air']
time_unit = time.units
In [188]:
print("These are the first 6 entries for time:",'\n',time[0:6])
In [189]:
precnc = nc.variables['air']
In [190]:
#Plot the data along a meridional line
#Plot a 90S-90N temp along a meridional line at 180E
nc_x = np.linspace(-90,90,73)
plt.plot(nc_x,precnc[0,:,71],ls="-",color='k' )
plt.xlabel('Latitude')
plt.ylabel('Temp[$^{\circ}$C]')
plt.title('90S-90N Temperature [$^{o}$C] Along a Meridional Line at 35E: Jan 1948')
plt.show()
In [191]:
precst = np.zeros((10512,845))
temp = np.reshape(precnc[0,:,:],(144*73))
print(np.shape(temp))
print(temp[0:6])
In [192]:
def spmat(x,y):
for i in range(0,845):
y[:,i] = np.reshape(x[i,:,:],(144*73))
return y
In [193]:
precst2 = spmat(precnc,precst)
print(np.shape(precst2))
In [194]:
#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
In [195]:
arra = np.zeros((0))
LAT = rep(arra,lat_vals,144)
LON = rep(arra,lon_vals[0],73)
print(LON)
In [196]:
def recon():
rslt = LON
for jj in range(0,143):
rslt = np.append(rslt,rep(arra,lon_vals[jj],73))
return rslt
In [197]:
LON2 = recon()
print("The last 6 entries are:",LON2[-6:])
In [198]:
#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))
In [199]:
print(np.shape(gpcpst))
print(gpcpst)
Convert the Julian day and hour into calendar mons for time¶
In [200]:
#create an array of all ones for the days of the data
ones = np.ones(845)
In [201]:
#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()
print(np.shape(months))
In [202]:
#define a functionto 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()
print(np.size(years))
In [203]:
#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,845):
empty = np.append(empty,str(years[i]) + '-' + str(months[i]))
return empty
tm1 = tmfunc()
print(tm1.shape)
In [204]:
#add the 'Lat' and 'Lon' columns to the beginning
#of our column names array
tm2 = np.append(['Lat','Lon'],tm1)
print(np.shape(tm2))
In [205]:
#use the package pandas.DataFrame to match the column names
#with the data
GPCPST = pd.DataFrame(gpcpst,columns=tm2)
GPCPST
Out[205]:
In [206]:
#save the data as a .csv file
GPCPST.to_csv('ncepJan1948_May2018.csv')
Compute the January climatology and standard deviation is below¶
In [207]:
monJ = np.array(range(0,844,12))
In [208]:
#choose the 3rd column and onward
gpcpdat = GPCPST.iloc[:,2:844]
gpcpJ = gpcpdat.iloc[:,monJ]
In [209]:
climJ = gpcpJ.mean(axis=1)
print(np.shape(climJ))
print(np.shape(gpcpJ))
In [210]:
def rowSDS():
rslt = []
for jj in range(0,10512):
rslt = np.append(rslt,np.std(gpcpJ.iloc[jj,:]))
return rslt
In [211]:
sdJ = rowSDS()
print(np.max(sdJ))
In [212]:
def repeat():
q = np.zeros((144,73))
u = climJ[0:73]
for ii in range(0,144):
q[ii,:] = u
return q
In [213]:
Lat1 = np.linspace(90.,-90.,73)
Lat = -Lat1
mapmat = np.reshape(np.array(climJ),(73,144))
print(np.shape(mapmat))
print(mapmat[:,0])
In [214]:
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAn = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data because
lons, data = ReAn.shiftdata(lon_vals, datain = mapmat, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAn.drawcoastlines(color='black', linewidth=1)
ReAn.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
ReAn.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
ReAn.drawcountries()
limit = np.linspace(-48.,48.,51)
ReAn_plt = plt.contourf(np.array(lons),np.array(lat_vals),data,limit,cmap = 'jet')
ceb = plt.colorbar(ReAn_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('NCEP Jan SAT RA 1948-2018 climatology [$^{o}$C]')
plt.show()
Plot Jan Standard Deviation¶
In [215]:
Lat1 = np.linspace(90.,-90.,73)
Lat = -Lat1
mapmat2 = np.reshape(sdJ,(73,144))
In [216]:
print(mapmat2.shape)
In [217]:
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
ReAnSd = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data
lons2, data2 = ReAnSd.shiftdata(lon_vals, datain = mapmat2, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
ReAnSd.drawcoastlines(color='black', linewidth=1)
ReAnSd.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
ReAnSd.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
ReAnSd.drawcountries()
limit = np.linspace(0.,8.5,100)
ReAnSd_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data2,limit,cmap = 'jet')
ceb = plt.colorbar(ReAnSd_plt)
ceb.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('NCEP Jan SAT RA 1948-2018 Standard Deviation [$^{o}$C]')
plt.savefig('NCEPStandardDeviation.png',bbox_inches='tight')
plt.show()
Compute the Jan EOFs¶
In [218]:
monJ = np.array(range(0,846,12))
gpcpdat = GPCPST.iloc[:,2:844]
gpcpJ = gpcpdat.iloc[:,monJ]
climJ = gpcpJ.mean(axis=1)
In [219]:
print(np.shape(gpcpdat.iloc[:,monJ]))
print(climJ.shape)
print(sdJ.shape)
In [220]:
def rowSDS2():
rslt = []
for jj in range(0,10512):
rslt = np.append(rslt,np.std(gpcpJ.iloc[jj,:]))
return rslt
In [221]:
sdJ = gpcpJ.std(axis=1)
In [222]:
sdJ.shape
sdJ.iloc[0:1000]
Out[222]:
In [223]:
def anomJ_func2():
rslt = pd.DataFrame(np.zeros((10512, 71)),columns=tm1[monJ])
for i in range(0,71):
rslt.iloc[:,i] = (gpcpdat.iloc[:,monJ[i]]-climJ)/sdJ
return rslt
In [224]:
anomJ = anomJ_func2()
In [225]:
anomJ.shape
anomJ
Out[225]:
In [226]:
def anomJ_func3():
rslt = pd.DataFrame(np.zeros((10512, 71)),columns=tm1[monJ])
for i in range(0,71):
rslt.iloc[:,i] = np.sqrt(np.cos(GPCPST.iloc[:,0]*np.pi/180.))*anomJ.iloc[:,i]
return rslt
In [227]:
anomAW = anomJ_func3()
anomAW.shape
anomAW
Out[227]:
In [228]:
#Computing SVD from the area-weighted anomaly matrix anomAW
U, S, Vh = np.linalg.svd(anomAW,full_matrices=False)
Plot the eigenvalues¶
In [229]:
xvals = np.linspace(0.,70.,71)
yvals = 100.*(S**2.)/np.sum(S**2)
#this creates two subplots
fig, ax1 = plt.subplots()
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.tick_params(axis='y',labelcolor='k')
ax2 = ax1.twinx()
ax2.scatter(xvals,yvals.cumsum(),marker = "o",facecolors='none',edgecolors='b')
ax2.plot(xvals,yvals.cumsum(),color='b',label = "Cumulative Variance")
ax2.set_ylabel('Cumulative Variance %',color='b')
ax2.tick_params(axis='y',labelcolor='b')
ax1.legend(loc=(1.3,.5))
ax2.legend(loc=(1.3,.4))
plt.show()
Plot the EOF1: The physical EOF by taking out the area-factor¶
In [230]:
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
In [231]:
mapmatr = (np.reshape(U[:,0]/(np.sqrt(np.cos(gpcpst[:,0]
*np.pi/180.))),(73,144)))
mapmatr2 = mapmatr[:,range(0,mapmatr[0,:].size)]
In [232]:
#Plot EOF1
plt.figure(figsize=(15., 7.))
#Using the basemap package to put the cylindrical projection map into the image
eof = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#shifting the data because
lons2, data2 = eof.shiftdata(lon_vals, datain = -mapmatr2, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
eof.drawcoastlines(color='black', linewidth=1)
eof.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
eof.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
eof.drawcountries()
limits = np.linspace(-.04,.04,51)
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data2,limits,cmap = 'jet')
efo = plt.colorbar(eof_plt)
efo.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('January EOF1 from 1948-2015 NCEP Temp Data')
plt.show()
In [233]:
print(GPCPST.iloc[:,0])
print(U[:,0])
In [234]:
pcdat = Vh[0,:]
print(pcdat)
In [235]:
#Plot PC1
TIME = np.array(range(1948,2019))
plt.figure()
plt.scatter(TIME,-pcdat,marker = "o",facecolors='none',edgecolors='k')
plt.plot(TIME,-pcdat,color='k')
plt.title("PC1 of NCEP TA Jan SAT: 1948-2018")
plt.xlabel("Year")
plt.ylabel("PC Values")
plt.show
Out[235]:
EOF from de-trended data¶
In [236]:
monJ = np.array(range(0,846,12))
gpcpdat = GPCPST.iloc[:,2:844]
gpcpJ = gpcpdat.iloc[:,monJ]
climJ = gpcpJ.mean(axis=1)
In [237]:
sdJ = gpcpJ.std(axis=1)
In [238]:
def anomJ_func2():
rslt = pd.DataFrame(np.zeros((10512, 71)),columns=tm1[monJ])
for i in range(0,71):
rslt.iloc[:,i] = (gpcpdat.iloc[:,monJ[i]]-climJ)/sdJ
return rslt
In [239]:
anomJ = anomJ_func2()
anomJ
Out[239]:
In [240]:
# Create linear regression object
regr = linear_model.LinearRegression()
print(anomJ.iloc[0,:].shape)
In [241]:
def trendV():
trendV = np.zeros(10512)
for i in range(0,10512):
regr.fit(np.reshape(TIME,(71,1)),
np.reshape(np.array(anomJ.iloc[i,:]),(71,1)))
trendV[i] = regr.coef_
return trendV
In [242]:
def trendM():
rslt = np.zeros((10512,71))
trendM = pd.DataFrame(np.zeros((10512, 71)),columns=tm1[monJ])
for i in range(0,10512):
regr.fit(np.reshape(TIME,(71,1)),
np.reshape(np.array(anomJ.iloc[i,:]),(71,1)))
rslt[i,:] = regr.intercept_
return rslt
In [243]:
trendV = trendV()
print(trendV)
In [244]:
print(regr.fit(np.reshape(TIME,(71,1)),
np.reshape(np.array(anomJ.iloc[0,:]),(71,1))))
In [245]:
trendM = trendM()
In [246]:
print(trendM)
In [ ]:
In [247]:
#vArea = np.cos(GPCPST.iloc[:,0]*np.pi/180.)
In [248]:
#anomA = vArea*anomJ
In [249]:
#JanSAT = anomA.sum(axis=0)/vArea.sum
In [250]:
#plt.figure(figsize=(12,7))
#plt.plot(TIME, JanSAT)
Plot the trend of Jan SAT non-standardized anomaly data. Begin with the space time data matrix¶
In [251]:
monJ = np.array(range(0,846,12))
gpcpdat = GPCPST.iloc[:,2:844]
gpcpJ = gpcpdat.iloc[:,monJ]
climJ = gpcpJ.mean(axis=1)
In [252]:
def anomJ_func3():
rslt = pd.DataFrame(np.zeros((10512, 71)),columns=tm1[monJ])
for i in range(0,71):
rslt.iloc[:,i] = gpcpdat.iloc[:,monJ[i]]-climJ
return rslt
In [253]:
anomJ = anomJ_func3()
In [254]:
def trendV():
trendV = np.zeros(10512)
for i in range(0,10512):
regr.fit(np.reshape(TIME,(71,1)),
np.reshape(np.array(anomJ.iloc[i,:]),(71,1)))
trendV[i] = regr.coef_
return trendV
In [255]:
trendV = trendV()
print(trendV)
In [256]:
mapmat1 = np.reshape([10*trendV],(73,144))
In [257]:
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
In [258]:
mapv1 = mapv1_func()
In [259]:
print(mapv1[:,0])
In [260]:
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
In [261]:
mapmat3 = mapmat_func()
In [262]:
print(np.max(mapmat3))
print(np.min(mapmat3))
In [263]:
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=-180.,urcrnrlon=180.,resolution='c')
#shifting the data because
lons2, data3 = trend.shiftdata(lon_vals, datain = mapmat3, lon_0=0)
#draw coastlines, latitudes, and longitutes on the map
trend.drawcoastlines(color='black', linewidth=1)
trend.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
trend.drawmeridians(np.arange(-180.,181.,60.), labels= [1,0,0,1])
trend.drawcountries()
limits = np.linspace(-1.1,1.1,50)
eof_plt = plt.contourf(np.array(lons2),np.array(lat_vals),data3,limits,cmap = 'jet')
efo = plt.colorbar(eof_plt)
efo.ax.set_ylabel('$^{o}$C')
plt.xlabel('Longitude',labelpad = 20)
plt.ylabel('Latitude',labelpad =20)
plt.title('Trend of the NCEP RA1 Jan 1948-2018 Anom Temp [$^{o}C/(10a)$]')
plt.show()
EOF analysis of the Tahiti and Darwin standardized SLP data SOI analysis¶
line 1241
Download pressure data for Tahiti from https://shen.sdsu.edu/data/PSTANDtahiti into your directory.¶
In [264]:
Pta = np.loadtxt("PSTANDtahiti.txt",skiprows=0)
In [265]:
PTA = np.reshape(Pta[:,1:],(780))
In [266]:
#make an array of years for our x-axis
xtime = np.linspace(1951,2016,(65*12)+1)[0:780]
In [267]:
plt.figure(figsize= (15.,7.))
plt.plot(xtime,PTA,color='r')
plt.xlabel('Year')
plt.ylabel('Pressure')
plt.title('Standardized Tahiti SLP Anomalies')
plt.show()
Download pressure data for Tahiti from http://shen.sdsu.edu/data/PSTANDdarwin.txt into your directory.¶
In [268]:
Pda = np.loadtxt('PSTANDdarwin.txt',skiprows=0)
PDA = np.reshape(Pda[:,1:],(780))
In [269]:
plt.figure(figsize=(15.,7.))
plt.plot(xtime,PDA,color='b')
plt.xlabel('Year')
plt.ylabel('Pressure')
plt.title('Standardized Darwin SLP Anomalies')
plt.show()
Plot the SOI index¶
In [270]:
fig, soi = plt.subplots(figsize=(15,7))
soi.plot(xtime,PTA-PDA,color='k')
soi.set_xlabel('Year')
soi.set_ylabel('SOI Index')
ax2 = soi.twiny()
ax2.plot(xtime,PTA-PDA,color='k')
ax2.set_xlabel('Year')
ax3 = soi.twinx()
ax3.plot(xtime,PTA-PDA,color='k')
ax3.set_ylabel('SOI Index')
plt.show()
Plot cumulative SOI: CSOI¶
In [271]:
SOI = PTA-PDA
cnegsoi = -SOI.cumsum()
In [272]:
plt.figure(figsize=(15,7))
plt.plot(xtime,cnegsoi,color='k')
plt.xlabel('Year')
plt.ylabel('Negative CSOI Index')
Out[272]:
In [273]:
ptada = np.append([PTA],[PDA], axis=0)
print(ptada.shape)
In [274]:
svdptd_u, svdptd_s, svdptd_vh = np.linalg.svd(ptada,full_matrices=False)
In [275]:
recontd = np.dot(np.dot(svdptd_u,np.diag(svdptd_s)),svdptd_vh)
In [276]:
print(recontd.shape)
print("Standard deviation of difference of arrays:",np.std(ptada-recontd))
Plot the WSOI 1
In [277]:
Dd = np.diag(svdptd_s)
Uu = svdptd_u
Vv = svdptd_vh
xtime = np.linspace(1951,2016,(65*12)+1)[0:780]
wsoi1 = Dd[0,0] * (Vv)[0,:]
In [278]:
print(wsoi1.shape)
In [279]:
plt.figure(figsize=(15,7))
plt.plot(xtime,wsoi1,color='k')
plt.xlabel('Year')
plt.ylabel('Weighted SOI 1')
plt.show()
Plot the WSOI 2
In [280]:
wsoi2 = Dd[1,1] * Vv[1,:]
print(wsoi2.shape)
In [281]:
fig, wsoi = plt.subplots(figsize=(15,7))
wsoi.plot(xtime,wsoi2,color='k')
wsoi.set_xlabel('Year')
wsoi.set_ylabel('SOI Index')
ax2 = wsoi.twiny()
ax2.plot(xtime,wsoi2,color='k')
ax2.set_xlabel('Year')
plt.show()
SVD Analysis for the Darwin-Tahiti Standardized SLP¶
In [282]:
# Remove the first column that is the year and create a vector
Pta = np.loadtxt("PSTANDtahiti.txt",skiprows=0)
PTA = np.reshape(Pta[:,1:],(780))
#Generate time ticks from Jan 1951 to Dec 2015
xtime = np.linspace(1951,2016,(65*12)+1)[0:780]
In [283]:
#Plot the Tahiti standardized SLP anomalies
plt.figure(figsize= (15.,7.))
plt.plot(xtime,PTA,color='r')
plt.xlabel('Year')
plt.ylabel('Pressure')
plt.title('Standardized Tahiti SLP Anomalies')
plt.show()
Do the same for Darwin SLP
In [284]:
Pda = np.loadtxt('PSTANDdarwin.txt',skiprows=0)
PDA = np.reshape(Pda[:,1:],(780))
plt.figure(figsize=(15.,7.))
plt.plot(xtime,PDA,color='b')
plt.xlabel('Year')
plt.ylabel('Pressure')
plt.title('Standardized Darwin SLP Anomalies')
plt.show()
Plot the SOI index
In [285]:
regr = linear_model.LinearRegression()
In [286]:
SOI = PTA - PDA
regr.fit(np.reshape(xtime,(780,1)), np.reshape(SOI,(780,1)))
ssoi = regr.predict(np.reshape(xtime,(780,1)))
In [287]:
fig, SOi = plt.subplots(figsize=(15,7))
SOi.plot(xtime,SOI,color='k')
SOi.set_xlabel('Year')
SOi.set_ylabel('SOI Index')
plt.plot(xtime,[0]*xtime.size,color='k')
ax2 = SOi.twiny()
ax2.plot(xtime,ssoi,color='r')
ax2.set_xlabel('Southern Oscillation Index')
ax3 = SOi.twinx()
ax3.plot(xtime,SOI,color='k')
plt.show()
Cumulative Negative SOI
In [288]:
CNegsoi = -SOI.cumsum()
In [289]:
plt.figure(figsize=(15,7))
plt.plot(xtime, CNegsoi,color='k')
plt.xlabel('Year')
plt.ylabel('Negative CSOI Index')
plt.title('Cumulative Southern Oscillation Index')
plt.show()
Space-time standardized SLP data matrix
In [290]:
ptada = np.append([PTA],[PDA], axis=0)
print(ptada.shape)
In [291]:
uu, ss, vvh = np.linalg.svd(ptada,full_matrices=False)
In [292]:
EOF = uu
print(EOF)
In [293]:
WSOI = np.dot(np.transpose(EOF),ptada)
print(WSOI.shape)
The two weighted SOI indices
In [294]:
WSOI1 = WSOI[0,:]
WSOI2 = WSOI[1,:]
In [295]:
fig, plott = plt.subplots(figsize=(15,7))
plott.plot(xtime,WSOI1,color='k',label='WSOI1')
plott.set_xlabel('Year')
plott.set_ylabel('WSOI Index')
ax2 = plott.twiny()
ax2.plot(xtime,WSOI2,color='b',label='WSOI2: 33% Variance')
ax2.set_xlabel('Weighted Southern Oscillation Index')
plott.legend(loc=(.05,.85))
ax2.legend(loc=(.05,.8))
plt.show()
Energy for each mode
In [296]:
Energy = ss
print(Energy)
Energy ratio
In [297]:
print('This is the Energy ratio:',(Energy**2.)/np.sum(Energy**2.))
In [298]:
PC = vvh
In [299]:
PC1 = PC[0,:]
PC2 = PC[1,:]
Plot the PCs
In [300]:
plt.figure(figsize=(15,7))
plt.plot(xtime,PC1,color='k',label='Principal Component 1')
plt.plot(xtime,PC2,color='b',label='Principal Component 2')
plt.xlabel('Year')
plt.ylabel('PC Index')
plt.title('Tahiti-Darwin Principal Components')
plt.legend(loc=(.05,.89))
plt.show()
Plot the Cumulative PCs
In [301]:
plt.figure(figsize=(15,7))
plt.plot(xtime,PC1.cumsum(),color='k',label='Principal Component 1')
plt.plot(xtime,PC2.cumsum(),color='b',label='Principal Component 2')
plt.xlabel('Year')
plt.ylabel('Cumulative PC Index')
plt.title('Cumulative Tahiti-Darwin Principal Components')
plt.legend(loc=(.75,.85))
plt.show()
Display the two ENSO modes on a world map
In [302]:
#choosing an arbitrary size for the image
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
wmap2 = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#draw coastlines, latitudes, and longitutes on the map
wmap2.drawcoastlines(color='black', linewidth=1)
wmap2.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
wmap2.drawmeridians(np.arange(-180.,181.,30.), labels= [1,0,0,1])
plt.scatter(131.,-12.,color='b',s=150)
plt.scatter(-129.,-30,color='r',s=150)
plt.text(100,-24, 'Darwin -0.79', color = 'b', fontsize = 22)
plt.text(-155,-42, 'Tahiti 0.61', color = 'r', fontsize = 22)
plt.text(-130,40, 'El Nino Southern Oscillation Model 1', color = 'm', fontsize = 30)
#plot
plt.show()
In [303]:
#choosing an arbitrary size for the image
plt.figure(figsize=(15., 7.))
#using the basemap package to put the cylindrical projection map into the image
wmap2 = Basemap(projection='cyl',llcrnrlat=-90.,urcrnrlat=90.,\
llcrnrlon=-180.,urcrnrlon=180.,resolution='c')
#draw coastlines, latitudes, and longitutes on the map
wmap2.drawcoastlines(color='black', linewidth=1)
wmap2.drawparallels(np.arange(-90.,91.,30.), labels = [1,0,0,1])
wmap2.drawmeridians(np.arange(-180.,181.,30.), labels= [1,0,0,1])
plt.scatter(131.,-12.,color='r',s=150)
plt.scatter(-129.,-30,color='r',s=150)
plt.text(105,-24, 'Darwin 0.61', color = 'r', fontsize = 22)
plt.text(-155,-42, 'Tahiti 0.79', color = 'r', fontsize = 22)
plt.text(-130,40, 'El Nino Southern Oscillation Model 2', color = 'm', fontsize = 30)
#plot
plt.show()
Plot principal components of the Darwin and Tahiti SLP
In [304]:
Pta = np.loadtxt("PSTANDtahiti.txt",skiprows=0)
PTA = np.reshape(Pta[:,1:],(780))
Pda = np.loadtxt('PSTANDdarwin.txt',skiprows=0)
PDA = np.reshape(Pda[:,1:],(780))
ptada = np.append([PTA],[PDA], axis=0)
svdptd_u, svdptd_s, svdptd_vh = np.linalg.svd(ptada,full_matrices=False)
In [305]:
dat = pd.DataFrame(svdptd_u,columns=['EOF1','EOF2'],index=['Tahiti','Darwin'])
In [306]:
dat
Out[306]:
In [307]:
lon = [131,231]
lat = [-24,-18]
In [308]:
plt.scatter(lon,lat)
plt.xlabel('Lon')
plt.ylabel('Lat')
plt.axis([0, 360, -90, 90])
plt.grid()
plt.show()
In [309]:
plt.figure(figsize=(12,7))
plt.scatter(np.linspace(1951,2015,780),svdptd_vh[0,:],marker = "o",facecolors='none',edgecolors='r')
plt.title("PC1 as the weighted SOI")
plt.xlabel("Year")
plt.ylabel("WSOI 1 Values")
plt.show
Out[309]:
Use NOAAGlobalTemp monthly time series as an example¶
In [310]:
d1 = np.loadtxt("aravg.ann.land_ocean.90S.90N.v4.0.1.201805.asc.txt"
,skiprows=0)
print(d1.shape)
In [311]:
columns = ['V1','V2','V3','V4','V5','V6']
D1 = pd.DataFrame(d1,columns=columns)
D1
Out[311]:
Extract the temperature data from 1880-2015 and show the data.
In [312]:
tmean18 = D1.iloc[0:139,1]
tmean18
Out[312]:
Linear trend
In [313]:
# Create linear regression object
regr = linear_model.LinearRegression()
In [314]:
yrtime18 = np.array(range(1880,2019))
reg8018 = regr.fit(np.reshape(yrtime18,(139,1)),np.reshape(tmean18,(139,)))
print("This is the coefficient:",reg8018.coef_)
print("This is the intercept:",[(reg8018.intercept_)])
In [315]:
# calculate the trendline for the data
z = np.polyfit(yrtime18, tmean18, 1)
p = np.poly1d(z)
Plot the temperature time series and its trend line.
In [316]:
plt.figure(figsize=(12,7))
plt.plot(yrtime18,tmean18,"k")
plt.plot(yrtime18,p(yrtime18),"r-")
plt.text(1900, .21,'Linear Temperature trend 0.67 $^{o}C$ per century',
color ="r", fontsize =15)
plt.title("Global Annual Mean Land and Ocean Surface Temperature Anomalies 1880-2015")
plt.ylabel("Temperature $^{o}C$")
plt.xlabel("Year")
plt.show()
Histogram and its density fit.
In [317]:
h, mids = np.histogram(tmean18)
In [318]:
from scipy.stats import norm
xfit = np.linspace(np.min(tmean18),np.max(tmean18),30)
areat = (np.diff(mids[0:2])*np.size(tmean18))
yfit = areat*(norm.pdf(xfit,np.mean(tmean18),np.std(tmean18)))
# yfit = areat*(mlab.normpdf(xfit,np.mean(tmean18),np.std(tmean18)))
In [319]:
plt.figure(figsize=(10,8))
plt.hist(tmean18,facecolor='w',edgecolor="k")
plt.plot(xfit,yfit)
plt.title("Histogram of 1880-2018 Temperature Anomalies")
plt.xlabel("Temperature anomalies")
plt.ylabel("Frequency")
plt.show()
Scatter plot
In [320]:
ust = pd.read_csv('USJantemp1951-2016-nohead.csv',header=None)
soi = pd.read_csv('soi-data-nohead.csv',header=None)
soid = np.array(soi.iloc[:,1])
soim = np.reshape(soid,(66,12))
soij = soim[:,0]
ustj = np.array(ust.iloc[:,2])
plt.figure(figsize=(12,7))
plt.scatter(soij,ustj,color='k')
#Create a linear trendline for the data
lin = np.polyfit(soij, ustj, 1)
linreg = np.poly1d(lin)
plt.plot(soij,linreg(soij),'r')
plt.ylabel('US Temp $^{o}$F')
plt.xlabel('SOI[dimensionless]')
plt.title("Scatter plot between Jan SOI and US temp",fontsize=15)
plt.show()
