Chapter 6: Covariance Matrices, EOFs, and PCs


The Python code was developed based on the R code written by
Samuel Shen, Distinguished Professor, San Diego State University, California, USA

Kaelia Okamura, Briana Ramirez, and Danielle Lafarga contributed to the Python code development.
Kaelia Okamura composed and finalized the code.

In [1]:
import pandas as pd
import numpy as np
import math
import statistics as stats
import scipy.stats as scistats
import matplotlib.pyplot as plt
import matplotlib.colors as color
import os
import sys
import statsmodels.api as sm
from pandas import read_table
from netCDF4 import Dataset as ds
import shapely
import warnings
warnings.filterwarnings("ignore")

from matplotlib import cm as cm1
from numpy.linalg import multi_dot
import netCDF4 as nc
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic
from matplotlib.patches import Polygon
from scipy.ndimage import gaussian_filter
In [2]:
# Style Dictionary to standarize plotting 
# scheme between different python scripts 
styledict = {'xtick.labelsize':21,
             'xtick.major.size':9,
             'xtick.major.width':1,
             'ytick.labelsize':21,
             'ytick.major.size':9,
             'ytick.major.width':1,
             'legend.framealpha':0.0,
             'legend.fontsize':15,
             'axes.labelsize':22,
             'axes.titlesize':22,
             'axes.linewidth':2,
             'figure.figsize':(12,8),
             'savefig.format':'jpg'}
plt.rcParams.update(**styledict)

# define font style for text
font = {'family': 'serif',
        'color':  'blue',
        'weight': 'normal',
        'size': 12,
        }
In [3]:
#Go to the working directory
import os
os.chdir("/Users/sshen/climstats")

Figure 6.1a: SAT covariance plot along a zonal band¶

In [4]:
import netCDF4 as nc
#import .nc data
nc = nc.Dataset("data/air.mon.anom.nc")
# extrat data for variables
lon = nc.variables['lon'][:]
lat = nc.variables['lat'][:]
time = nc.variables['time'][:]
sat = nc.variables['air'][:] 
#sat = surface air temp anomalies: mon, 5deg grid
# 1674 months from Jan 1880-June 2019: 139yrs 6mons
print(np.shape(sat)) #(1674, 36, 72): mon, lat, lon
#print(nc) 
Dec = list(np.arange(11, 1674, 12))#Dec index
Decsat = sat[Dec, :, :] #extract Dec data
# space-time data matrix STsat
N = 72 * 36
P = len(Dec)
STsat = np.zeros((N,P))
for k in range(P):
    STsat[:,k] = Decsat[k,:,:].flatten()
print(STsat.shape)#(2592, 139) 139 Decembers
#print(STsat[0:4, 0:4])

#Space-time dataframe with lat, lon, time
# repeat each value in lat 72 times
LAT = np.repeat(lat,72).reshape(2592,1) 
# repeat lon 36 times
LON = np.tile(lon,36).reshape(2592,1)
STanom0 = np.hstack((LAT, LON, STsat))
yrlist = list(np.arange(1880, 2019))
cnames = ["LAT", "LON"] + yrlist
STanom = pd.DataFrame(STanom0, columns=cnames)

#Select the zonal band
n1 = np.where((STanom["LAT"] > -4) & (STanom["LAT"] < 0)\
          & (STanom["LON"] > 0) & (STanom["LON"] < 360))
P1 = 81 #Dec 1961 in STsat space-time data matrix
P2 = 130 #Dec 2010
P = P2 - P1 + 1 #50 years
print(STsat.shape) #The space-time data
da0 = STsat[n1,P1:(P2 + 1)]
da1 = da0[0,:,:]

(1674, 36, 72)
(2592, 139)
(2592, 139)
In [5]:
STanom #Display the space-time data frame
Out[5]:
LAT LON 1880 1881 1882 1883 1884 1885 1886 1887 ... 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018
0 -87.5 2.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
1 -87.5 7.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
2 -87.5 12.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
3 -87.5 17.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
4 -87.5 22.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2587 87.5 337.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
2588 87.5 342.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
2589 87.5 347.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
2590 87.5 352.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36
2591 87.5 357.5 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 ... -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36 -9.969210e+36

2592 rows × 141 columns

In [6]:
# define AreaFrac as a float for future manipulation
AreaFrac = np.sqrt(np.cos(LAT[n1]*np.pi/180))
#AreaFrac = float(np.sqrt(np.cos(Lat1[0]*np.pi/180)))

# subtract row means from dat1 to get math anomalies
dat2 = da1 - da1.mean(axis=1, keepdims=True)
# multiply each value by AreaFrac
Banddat = AreaFrac * dat2
#print(Banddat.shape)
covBand = np.dot(Banddat, Banddat.transpose())
# check max and min values
print(np.max(covBand))
print(np.min(covBand))
covBand62 = covBand #To be used for Fig 6.2 later

90.67207862062392
-6.18168468483208
In [7]:
#plot the figure: Fig. 6.1a
int = np.linspace(-7,92,81) # define levels for plot
fig, ax = plt.subplots()
# define color bar
myColMap = LinearSegmentedColormap.from_list(name='my_list', 
            colors=['navy','darkgreen','lime','yellow',
                 'pink','red','maroon'], N=100)
# plot covariance plot 
Lon1 = np.arange(2.5, 360, 5)
colormap = plt.contourf(Lon1,Lon1,covBand, 
                        levels = int, cmap = myColMap)
plt.title("Covariance of SAT Anomalies on \
Zonal Band 2.5$\degree$S",
          size = 21, pad = 15)
plt.xlabel("Longitude", labelpad = 15)
plt.ylabel("Longitude", labelpad = 15)
# plot color bar
cbar = plt.colorbar(colormap, shrink=0.9, pad = 0.05)
plt.text(372, 348, '[$\degree$C]$^2$', size = 20)
cbar.set_ticks([0,20,40,60,80])
plt.savefig("fig0601a.eps")# save figure

Figure 6.1b: SAT covariance plot along a meridional band¶

In [8]:
#Covert missing values into NaN
STanom[STanom < -999] =  pd.NA
STanom.iloc[0:2, 0:5]
Out[8]:
LAT LON 1880 1881 1882
0 -87.5 2.5 NaN NaN NaN
1 -87.5 7.5 NaN NaN NaN
In [9]:
plt.figure(figsize=(9, 7))

n1 = np.where((STanom.iloc[:, 0] > -90) & (STanom.iloc[:, 0] < 90) 
              & (STanom.iloc[:, 1] > 0) & (STanom.iloc[:, 1] < 4))
n1 = np.concatenate(n1)

P1 = 83  # Time index for Dec 1961
P2 = 133  # Dec 2010
P = P2 - P1   # =50 Decembers
dat1 = STanom.iloc[n1, P1:P2]  # 1961-2010
Lat1 = STanom.iloc[n1, 0]
Lon1 = STanom.iloc[n1, 1]

AreaFac = np.sqrt(np.cos(Lat1 * np.pi / 180))
dat2 = dat1.transpose() - np.array(np.mean(dat1, axis=1))
Banddat = AreaFac[:, np.newaxis] * dat2.transpose()
covBand = np.dot(Banddat, Banddat.T)

#print(dat2)
# Display max and min of covariance matrix
print("Max:", np.nanmax(covBand))
print("Min:", np.nanmin(covBand))

# Creating the filled contour plot
intervals = np.linspace(-35, 120, 81)
myColMap = LinearSegmentedColormap.from_list(name='my_list', 
            colors=['navy','darkgreen','lime','yellow',
                 'pink','red','maroon'], N=100)

ticks = np.arange(-90, 91, 30)
plt.contourf(Lat1, Lat1, covBand, levels=intervals, cmap=myColMap)

clb = plt.colorbar(label='', ticks=np.arange(-40, 120, 20),
            pad = 0.05, shrink=0.9)
clb.ax.set_title('[°C]²')
plt.title("SAT Meridional Covariance: 2.5°E")
plt.xlabel("Latitude")
plt.ylabel("Latitude")
plt.xticks(ticks)
plt.yticks(ticks)
plt.grid()
plt.show()

Max: 116.8833093892827
Min: -34.57828535281515

Figure 6.1c: Positions of the zonal band in Fig. 6.1a and meridional band in Fig. 6.1b indicated by the thick dashed black lines¶

In [10]:
sat.shape
mapmat = sat[1416, :, :]
sat[1410:1419, 0, 0]
Out[10]:
masked_array(data=[--, --, --, --, --, --, --, --, --],
             mask=[ True,  True,  True,  True,  True,  True,  True,  True,
                    True],
       fill_value=-9.96921e+36,
            dtype=float32)
In [11]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as mcolors
import netCDF4 as nc

plt.figure(figsize=(12, 6))

#import .nc data
nc = nc.Dataset("data/air.mon.anom.nc")
# extrat data for variables
lon = nc.variables['lon'][:]
lat = nc.variables['lat'][:]
time = nc.variables['time'][:]
sat = nc.variables['air'][:] #The temp anomalies data

# sat is a 3D array, and lon, lat are longitude and latitude 
mapmat = sat[1415, :, :] #Dec 1997
mapmat = np.clip(mapmat, -5, 5)

# Create the color palette
rgb_palette = mcolors.LinearSegmentedColormap.from_list(
    "custom", ['black', 'blue', 'darkgreen', 'green', 'white', 
               'yellow', 'pink', 'red', 'maroon'])

# Initialize the Basemap
m = Basemap(projection='cyl', llcrnrlat=-90, 
            urcrnrlat=90, llcrnrlon=0, urcrnrlon=360)

# Create meshgrid for plotting
lon, lat = np.meshgrid(lon, lat)
x, y = m(lon, lat)

# Plot the data
cf = m.contourf(x, y, mapmat, levels=np.linspace(-5, 5, 31), 
               cmap=myColMap, extend='both')

# Add map features
m.drawcoastlines()
m.drawparallels(np.arange(-90., 91., 30.), 
                labels=[1,0,0,0], fontsize = 15)
m.drawmeridians(np.arange(0., 361., 60.), 
                labels=[0,0,0,1], fontsize = 15)

# Add the two thick dashed lines
plt.plot([0,360], [-2.5, 2.5], 
         linestyle='dashed', color = 'k', linewidth = 5)
plt.plot([2.5, 2.5], [-90, 90], 
         linestyle='dashed', color = 'k', linewidth = 5)

# Title and labels
plt.title("December 1997 SAT Anomalies", fontsize=20)
clb = plt.colorbar(cf, label='', ticks=np.arange(-5, 5.1, 1),
            pad = 0.03, shrink=0.8)
clb.ax.set_title('°C')

# Save the figure
plt.savefig("fig0601c.eps", format='eps')
plt.show()