Python气象数据处理与绘图：两个变量场的相关分析-SVD分解

1、前言

奇异值分解（Singular Value Decomposition，SVD），是一种提取信息的方法。
SVD和主成分分析一样，也是告诉我们数据中重要特征，奇异值是数据矩阵乘以该矩阵的转置的特征值的平方根。

2、一个简单的例子

 1import numpy as np
 2data = np.array([[0,0,0,2,2],
 3                 [0,0,0,3,3],
 4                 [0,0,0,1,1],
 5                 [1,1,1,0,0],
 6                 [2,2,2,0,0],
 7                 [5,5,5,0,0],
 8                 [1,1,1,0,0]])
 9u, sigma, vt = np.linalg.svd(data)
10sigma

1array([9.64365076e+00, 5.29150262e+00, 7.07756104e-16, 9.32863923e-17,
2       2.04123199e-33])

sigma是除了对角元素不为0，其他元素都为0，所以返回的时候，作为一维矩阵返回。

3、SVD在气象中的应用

1!pip install git+https://github.com/Yefee/xMCA.git
2Collecting git+https://github.com/Yefee/xMCA.git
3  Cloning https://github.com/Yefee/xMCA.git to /tmp/pip-req-build-bl2v33uq

1.导入模块

1from xMCA import xMCA
2import xarray as xr
3import matplotlib.pyplot as plt
4%matplotlib inline

2.读取数据
（1）读取美国2m气温距平场

1usta = xr.open_dataarray('/home/kesci/work/data/USTA.nc').transpose(*['time', 'lat', 'lon'])
2usta.name = 'USTA'
3usta

（2）太平洋海表温度距平场

1sstpc = xr.open_dataarray('/home/kesci/work/data/SSTPac.nc').transpose(*['time', 'lat', 'lon'])
2sstpc.name = 'SSTPC'
3sstpc

4、气温和海温的SVD分解

分解并得到前两个分量

1sst_ts = xMCA(sstpc, usta)
2sst_ts.solver()
3lp, rp = sst_ts.patterns(n=2)
4le, re = sst_ts.expansionCoefs(n=2)
5frac = sst_ts.covFracs(n=2)
6print(frac)

绘图：

1fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(18, 12))
2lp[0].plot(ax=ax1[0])
3le[0].plot(ax=ax1[1])
4
5rp[0].plot(ax=ax2[0])
6re[0].plot(ax=ax2[1])

气温和海温的协同回归

1lh, rh = sst_ts.homogeneousPatterns(n=1)
2le, re = sst_ts.heterogeneousPatterns(n=1)

绘图：

1fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(18, 12))
2lh[0].plot(ax=ax1[0])
3rh[0].plot(ax=ax1[1])
4le[0].plot(ax=ax2[0])
5re[0].plot(ax=ax2[1])

计算显著性水平（学生t检验）

1le, re, lphet, rphet = sst_ts.heterogeneousPatterns(n=1, statistical_test=True)

绘图：

1fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(18, 12))
2le[0].plot(ax=ax1[0])
3re[0].plot(ax=ax1[1])
4lphet[0].where(lphet[0]<0.01).plot(ax=ax2[0])
5rphet[0].where(rphet[0]<0.01).plot(ax=ax2[1])