
######################################################################################################################
#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


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

Python Code Version 3.0 for¶

Climate Mathematics: Theory and Applications - Chapter 11¶

Chapter 11: Python Analysis of Incomplete Climate Data¶

11.2 Read NOAAGlobalTemp and From the Space-Time Data Matrix¶

11.3 Spatial Averages and Their Trends¶