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
warnings.filterwarnings("ignore")
from scipy import optimize as opt
import sympy as sm
import math as m
import os
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
# import cartopy.crs as ccrs
# import cartopy.feature as cfeature
from scipy.ndimage import gaussian_filter
from sklearn.linear_model import LinearRegression
Starting with Figure 11.4, due to the fact that Figure 11.1-11.3 are snapshots of the NCDC NOAA data access tool. Please follow link below for more information regarding these plots. www.ncdc.noaa.gov/data-access/marineocean-data/noaa-global-surface-temperature-noaaglobaltemp
# open .asc file
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)
4267130
print(da1[0 : 3])
# create time variables
tm1 = np.arange(0, 4267129, 2594)
tm2 = np.arange(1, 4267130, 2594)
print(len(tm1))
print(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])))
print(mm1[0 : 6])
print(yy1[0 : 6])
print(len(mm1))
print(len(yy1))
print(tm1[0 : 6])
print(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]
print(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
[1.0, 1880.0, -999.9] 1645 1645 [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] [1880.0, 1880.0, 1880.0, 1880.0, 1880.0, 1880.0] 1645 1645 [ 0 2594 5188 7782 10376 12970] [ 1 2595 5189 7783 10377 12971] 1645.0
# 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
#Fig 11.4
# 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 = 15)
dmap.drawmeridians(np.arange(50, 350, 50.), labels = [0, 0, 0, 1],
fontsize = 15)
# 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%", ticks=[-6,-4,-2,0,2,4,6])
cbar.ax.set_yticklabels(['-6','-4','-2','0','2','4','6'])
plt.text(362, 90, '[$\degree$C]', size=20) # add colorbar title string
# add plot title
plt.title('NOAAGlobalTemp Anomalies [$\degree$C]: Dec. 2015')
# label x and y
plt.xlabel('Longitude', labelpad = 30)
plt.ylabel('Latitude', labelpad = 40)
# display on screen
plt.show()
(36, 72) (36, 72) (36, 72)
# 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]
(2592, 1647) 263520
#Fig 11.5
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.drawcountries()
# create and draw meridians and parallels grid lines
smap.drawparallels(np.arange( -40, 60, 20.),labels=[1, 0, 0, 0],fontsize = 15)
smap.drawmeridians(np.arange(50, 350, 50.),labels=[0, 0, 0, 1],fontsize = 15)
# convert latitude/longitude values to plot x/y values
x, y = smap(*np.meshgrid(Lon, Lat))
# contour levels
clevs = np.arange(-4, 4.2, 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%", ticks=[-4,-2,0,2,4])
# add colorbar title string
plt.text(302, 42, '[$\degree$C]', size=20)
# add plot title
plt.title('Tropic Pacific SAT Anomalies [$\degree$C]: Dec. 1997')
# label x and y
plt.xlabel('Longitude', labelpad = 30)
plt.ylabel('Latitude', labelpad = 40)
# display on screen
plt.show()
# prepare the data
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]
# 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])
#Fig 11.6
plt.close()
# create time variable
timemo = np.linspace(1880, 2017, 1645)
# create plot and plot labels
plt.plot(timemo, avev, '-', color="k")
plt.title('Area-Weighted Global Average of Monthly \n SAT Anomalies: Jan. 1880-Jan. 2017')
plt.xlabel('Year')
plt.ylabel('Temperature Anomaly [$\degree$C]')
plt.yticks([-1,-0.5,0,0.5,1])
# create linear regression model
lm = LinearRegression()
lm.fit(timemo.reshape(-1, 1), avev)
predictions = lm.predict(timemo.reshape(-1, 1))
# plot linear regression model
plt.plot(timemo, predictions, '-', color="blue",linewidth=2)
# add text indicating the trend
plt.text(1885, .55, 'Linear Trend: 0.69 [$\degree$C] per Century', fontsize=16, color="blue")
plt.grid()
plt.show()