Chapter 5: Matrices for Climate Data


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

Kaelia Okamura, Joseph Diaz, Liu Yang, and Briana Ramirez contributed to the Python code development.
Kaelia Okamura composed and finalized the code.

In [1]:
import numpy as np
import math
import statistics as stats
import scipy.stats as scistats
import matplotlib.pyplot as plt
import os
from pandas import read_table, read_csv
from netCDF4 import Dataset as ds
import matplotlib.ticker as mticker
import matplotlib.patches as patches

plt.rcParams['figure.figsize'] = (12, 12)
In [2]:
#Style Dictionary to standarize plotting scheme between different python scripts 
styledict = {'xtick.labelsize':25,
             'xtick.major.size':9,
             'xtick.major.width':1,
             'ytick.labelsize':25,
             '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)
In [3]:
import os
os.chdir("/Users/sshen/climstats")

Computational examples of matrices¶

In [4]:
# Matrix multiplication
A = [[1, 0],  [0, 4],[3, 2]]
B = [[0,-1],[1,2]]
C = np.matmul(A,B) #Or C = np.dot(A,B)
print('C=', C)
#C= [[ 0 -1]
# [ 4  8]
# [ 2  1]]
print('Transpose matrix of C =', C.transpose())
#Transpose matrix of C = [[ 0  4  2]
# [-1  8  1]]

C= [[ 0 -1]
 [ 4  8]
 [ 2  1]]
Transpose matrix of C = [[ 0  4  2]
 [-1  8  1]]
In [5]:
#matrix inversion
A = [[1,-1],[1,2]]
np.linalg.inv(A)# compute the inverse of A
#array([[ 0.66666667,  0.33333333],
#       [-0.33333333,  0.33333333]])
Out[5]:
array([[ 0.66666667,  0.33333333],
       [-0.33333333,  0.33333333]])
In [6]:
#Solve a system of linear equations
A = [[30, 40],[1, 1]]
b = [[1000],[30]]
x = np.linalg.solve(A,b)
print('x=',x)
#x= [[20.]
# [10.]]

x= [[20.]
 [10.]]
In [7]:
#Compute determinant of the previous matrix A
A = [[30, 40],[1, 1]]
print('Determinant det(A)=', np.linalg.det(A))
#Determinant det(A)= -10.000000000000002

Determinant det(A)= -10.000000000000002
In [8]:
#An orthogonal matrix
p = np.sqrt(2)/2
Q = [[p,p],[-p,p]]
print('Orthogonal matrix Q=', np.round(Q,2))
T = np.transpose(Q)
print('Q times transpose of Q = ', np.matmul(Q,T))
print('Determinant of Q =', np.linalg.det(Q))
#Orthogonal matrix Q= [[ 0.71  0.71]
# [-0.71  0.71]]
#Q times transpose of Q =  [[1. 0.]
# [0. 1.]]
#Determinant of Q = 1.0 

Orthogonal matrix Q= [[ 0.71  0.71]
 [-0.71  0.71]]
Q times transpose of Q =  [[1. 0.]
 [0. 1.]]
Determinant of Q = 1.0

Compute eigenvectors and eigenvalues¶

In [9]:
#Compute eigenvectors and eigenvalues
A = [[1,2], [2,1]] 
np.linalg.eig(A)
#(array([ 3., -1.]), array([[ 0.70710678, -0.70710678],
#        [ 0.70710678,  0.70710678]]))
Out[9]:
(array([ 3., -1.]),
 array([[ 0.70710678, -0.70710678],
        [ 0.70710678,  0.70710678]]))

Fig. 5.2: An eigenvector vs a non-eigenvector¶

In [10]:
import matplotlib.patches as patches
fig = plt.figure(figsize = (12,12))
x = [0]
y = [0]
plt.plot(x,y, 'ko-')
plt.xlim([-0.1,3.1])
plt.ylim([-0.1,3.1])

style = "Simple, tail_width=0.5,head_width=8,head_length=28"
kw1 = dict(arrowstyle=style, color="blue")
kw2 = dict(arrowstyle=style, color="red")

a1 = patches.FancyArrowPatch((0,0), (1,0), **kw1, 
                             linewidth =5)
a2 = patches.FancyArrowPatch((0, 0), (1,2),**kw1)
a3 = patches.FancyArrowPatch((0, 0), (1,1), **kw2, 
                             linewidth = 5)
a4 = patches.FancyArrowPatch((0, 0), (3,3),**kw2)
for a in [a1, a2, a3, a4]:
    plt.gca().add_patch(a)
plt.title('An eigenvector v vs a non-eigenvector u')
plt.xlabel(r'$x_1$', fontsize = 25)
plt.ylabel(r'$x_2$', fontsize = 25)
plt.text(0.6, 0.1, 'Non-eigenvector u',
          color = 'blue', fontsize =25)
plt.text(0.8, 2.0, 'Au',
          color = 'blue', fontsize =25)
plt.text(1.03,0.85, 'Eigenvector v',
          color = 'red', fontsize =25)
plt.text(2.7, 2.9, 'Av',
          color = 'red', fontsize =25)
plt.show()

Diagonalization and matrix decomposition¶

In [11]:
# Verify diagonalization and decomposition
C = [[2,1],[1,2]]
valC, Q = np.linalg.eig(C)
print('eigenvalues of C =', valC)
#eigenvalues of C = [3. 1.]
print('eigenvectors of C =', Q)
#eigenvectors of C = [[ 0.70710678 -0.70710678]
# [ 0.70710678  0.70710678]]
D = Q.transpose(1,0).dot(C).dot(Q)
print('D = ', D)
#D =  [[3. 0.]
# [0. 1.]]
# Matrix C is decomposed into three matrices: C = Q D Q'
Q.dot(D).dot(Q.transpose(1,0))
#array([[2., 1.],
#       [1., 2.]]) 

eigenvalues of C = [3. 1.]
eigenvectors of C = [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
D =  [[3. 0.]
 [0. 1.]]
Out[11]:
array([[2., 1.],
       [1., 2.]])

Decompose a matrix into vectors¶

In [12]:
D[0][0]*np.outer(Q[:][0],Q.transpose()[:][0]) + \
D[1][1]*np.outer(Q[:][1],Q.transpose()[:][1]) 
#array([[ 1.,  2.],
#       [-2., -1.]]) #matrix C is recovered from vectors
Out[12]:
array([[ 1.,  2.],
       [-2., -1.]])

SVD example for a 2-by-3 matrix¶

In [13]:
#SVD example for a 2-by-3 matrix: Python code
A = [[-1,0, -2],[1,2,3]]
UsvdA, DsvdA, VsvdA = np.linalg.svd(A)
print('Singular values= ', np.round(DsvdA,2))
#Singular values=  [4.22 1.09]
print('Spatial singular vectors= ', np.round(UsvdA,2))
#Spatial singular vectors=  [[-0.48  0.88]
# [ 0.88  0.48]]
print('Temporal singular vectors= ', np.round(VsvdA,2))
#Temporal singular vectors=  [[ 0.32  0.42  0.85]
# [-0.37  0.88 -0.29]
# [-0.87 -0.22  0.44]]

Singular values=  [4.22 1.09]
Spatial singular vectors=  [[ 0.48  0.88]
 [-0.88  0.48]]
Temporal singular vectors=  [[-0.32 -0.42 -0.85]
 [-0.37  0.88 -0.29]
 [-0.87 -0.22  0.44]]
In [14]:
B = np.array(A)
C = np.matmul(B, B.T) #B times B transpose
valC, vecC = np.linalg.eig(C)
np.sqrt(valC)
#array([1.0855144 , 4.22157062])
Out[14]:
array([1.0855144 , 4.22157062])
In [15]:
#Data reconstruction by singular vectors: Python code
np.round(DsvdA[0]*np.outer(UsvdA[:][0], VsvdA[:][0]),1)
#array([[-0.7, -0.8, -1.7],
#       [ 1.2,  1.5,  3.2]])
Out[15]:
array([[-0.7, -0.8, -1.7],
       [-1.2, -1.5, -3.2]])
In [16]:
A1 = DsvdA[0]*np.outer(UsvdA[:][0], VsvdA[:][0])
A2 = DsvdA[1]*np.outer(UsvdA[:][1], VsvdA[:][1])
np.round(A1 + A2, 2) 
#array([[-1.,  0., -2.],
#       [ 1.,  2.,  3.]])
Out[16]:
array([[-0.3 , -1.68, -1.44],
       [-1.38, -1.08, -3.3 ]])

Fig. 5.3: Schematic diagram of SVD¶

In [17]:
import matplotlib.pyplot as plt
import matplotlib.patches as patches

def draw_matrix(ax, lower_left, upper_right, color, 
                n_label, m_label, label_color, matrix_label):
    # Draw the rectangle
    ax.add_patch(patches.Rectangle(lower_left, 
                                   upper_right[0] - lower_left[0], 
                                   upper_right[1] - lower_left[1], 
                                   edgecolor=color, 
                                   facecolor='none', lw=3))
    
    # Add n and m labels
    ax.text(lower_left[0] - 0.8, (lower_left[1] + upper_right[1]) / 2, 
            n_label, ha='right', va='center', 
            fontsize=20, color=label_color, rotation=90)
    ax.text((lower_left[0] + upper_right[0]) / 2, 
            upper_right[1] + 0.8, m_label, ha='center', va='bottom', 
            fontsize=20, color=label_color)

    # Add matrix label
    ax.text((lower_left[0] + upper_right[0]) / 2, 
            lower_left[1] - 1.5, matrix_label, 
            ha='center', va='top', fontsize=20)

def draw_diagonal(ax, start, end, color):
    ax.plot([start[0], end[0]], [start[1], end[1]], 
            color=color, lw=1.3, linestyle='--')

def main():
    # Set up the plot
    fig, ax = plt.subplots(figsize=(15, 10))
    ax.set_xlim(-3, 28)
    ax.set_ylim(-3, 16)
    ax.axis('off')  # Hide axis

    # SVD Text
    ax.text(13, 15.5, 
            r"SVD: $A = UDV^T$ when n > m (top panel) or n < m (bottom panel)", 
            ha='center', fontsize=22)
    
    #Draw the large equality sign =
    plt.text(4, 10, '=', color = 'k',
           fontsize =40, rotation=0)

    # Draw matrices for n > m
    draw_matrix(ax, (0, 7), (3, 13), 'blue', 'n', 'm', 'blue', 'A[n×m]')
    draw_matrix(ax, (7, 7), (10, 13), 'blue', 'n', 'm', 'blue', 'U[n×m]')
    draw_matrix(ax, (12, 10), (15, 13), 'brown', 'm', 'm', 'red', 'D[m×m]')
    draw_diagonal(ax, (12, 13), (15, 10), 'brown')
    draw_matrix(ax, (17, 10), (20, 13), 'red', 'm', 'm', 'red', "(Vᵀ)[m×m]")

    # Draw matrices for n < m
    draw_matrix(ax, (0, 0), (6, 3), 'blue', 'n', 'm', 'blue', 'B[n×m]')
    draw_matrix(ax, (11, 0), (14, 3), 'blue', 'n', 'n', 'blue', 'U[n×n]')
    draw_matrix(ax, (16, 0), (19, 3), 'brown', 'n', 'n', 'blue', 'D[n×n]')
    draw_diagonal(ax, (16, 3), (19, 0), 'brown')
    draw_matrix(ax, (21, 0), (27, 3), 'red', 'n', 'm', 'red', "(Vᵀ)[n×m]")
    
     #Draw the large equality sign =
    plt.text(7, 1.1, '=', color = 'k',
           fontsize =40, rotation=0)
    
#Draw the dots for matrix columns and rows for the upper panel
    plt.plot([0.5,0.5], [7, 13], 
         linestyle='dotted', color = 'k')
    plt.plot([1, 1], [7, 13], 
        linestyle='dotted', color = 'k')
    plt.text(1.7, 10, r'$\cdots$', color = 'k',
           fontsize =25, rotation=0)
    
    plt.plot([0.5+7,0.5+7], [7, 13], 
         linestyle='dotted', color = 'b')
    plt.plot([1+7, 1+7], [7, 13], 
        linestyle='dotted', color = 'b')
    plt.text(1.7+7, 10, r'$\cdots$', color = 'blue',
           fontsize =25, rotation=0)
    
    
    plt.text(1.7+7+5, 10+2, '0', color = 'maroon',
           fontsize =25, rotation=0)
    plt.text(1.7+7+3.7, 10+0.8, '0', color = 'maroon',
           fontsize =25, rotation=0)
    
    plt.plot([17,20], [12, 12], 
         linestyle='dotted', color = 'r')
    plt.plot([17,20], [12.5, 12.5], 
         linestyle='dotted', color = 'r')
    plt.text(18, 11, r'$\vdots$', color = 'red',
           fontsize =25, rotation=0)
    
#Draw the dots for matrix columns and rows for the lower panel
    plt.plot([0, 6], [2.5, 2.5], 
         linestyle='dotted', color = 'k')
    plt.plot([0, 6], [2.0, 2.0], 
         linestyle='dotted', color = 'k')
    plt.text(2.7, 1.0, r'$\vdots$', color = 'k',
           fontsize =25, rotation=0)
    
    plt.plot([12, 12], [0, 3], 
         linestyle='dotted', color = 'b')
    plt.plot([12-0.5, 12-0.5], [0, 3], 
         linestyle='dotted', color = 'b')
    plt.text(1.7+11, 1.2, r'$\cdots$', color = 'blue',
           fontsize =25, rotation=0)
    
    
    plt.text(1.7+16.1, 1.6, '0', color = 'maroon',
           fontsize =25, rotation=0)
    plt.text(1.7+14.9, 0.7, '0', color = 'maroon',
           fontsize =25, rotation=0)
    
    plt.plot([21, 27], [2, 2], 
         linestyle='dotted', color = 'r')
    plt.plot([21, 27], [2.5, 2.5], 
         linestyle='dotted', color = 'r')
    plt.text(23.5, 1, r'$\vdots$', color = 'red',
           fontsize =25, rotation=0)

    # Show the plot
    plt.show()

if __name__ == "__main__":
    main()

SVD analysis of Darwin and Tahiti SLP¶

In [18]:
#Python SVD analysis of Darwin and Tahiti SLP
import os
os.chdir("/Users/sshen/climstats")
PDA = np.array(read_table("data/PSTANDdarwin.txt", \
                          header = None, delimiter = "\s+"))
PTA = np.array(read_table("data/PSTANDtahiti.txt", \
                          header = None, delimiter = "\s+"))
pdata = np.stack([PDA[58:65,12], PTA[58:65,12]], axis=0) 
print('The Darwin and Tahiti SLP data 2009-2015 =', pdata)
#The Darwin and Tahiti Standardized SLP anomalies = 
#[[ 0.5 -2.3 -2.2  0.3  0.3  0.1 -0.4]
# [-0.7  2.5  1.9 -0.7  0.4 -0.8 -1.3]]
u, d, v = np.linalg.svd(pdata)
print('Spatial singular vectors EOFs U =', np.round(u,2))
#Spatial singular vectors EOFs U = [[-0.66  0.75]
# [ 0.75  0.66]]
print('Diagonal matrix D =', np.round(np.diag(d),2))
#Diagonal matrix D = [[4.7  0.  ]
# [0.   1.42]]
print('Temporal singular vectors PCs V=', np.round(v,2))
#Temporal singular vectors PCs V= 
#[[-0.18  0.72  0.61 -0.15  0.02 -0.14 -0.15]
# [-0.06 -0.06 -0.28 -0.17  0.34 -0.32 -0.82] ...]

The Darwin and Tahiti SLP data 2009-2015 = [[ 0.5 -2.3 -2.2  0.3  0.3  0.1 -0.4]
 [-0.7  2.5  1.9 -0.7  0.4 -0.8 -1.3]]
Spatial singular vectors EOFs U = [[ 0.66  0.75]
 [-0.75  0.66]]
Diagonal matrix D = [[4.7  0.  ]
 [0.   1.42]]
Temporal singular vectors PCs V= [[ 0.18 -0.72 -0.61  0.15 -0.02  0.14  0.15]
 [-0.06 -0.06 -0.28 -0.17  0.34 -0.32 -0.82]
 [ 0.74 -0.32  0.55  0.03  0.1  -0.02 -0.18]
 [ 0.01  0.19 -0.03  0.97  0.06 -0.06 -0.15]
 [-0.27 -0.18  0.2   0.04  0.87  0.1   0.29]
 [ 0.15  0.25 -0.13 -0.05  0.12  0.9  -0.27]
 [ 0.56  0.49 -0.43 -0.1   0.3  -0.24  0.32]]