奇异值分解(SVD)

Function

本文将直接介绍该库的安装及使用，关于SVD的原理不做介绍。

一、安装

1
pip install xMCA

二、使用介绍

首先import

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from xMCA import xMCA
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import cartopy.mpl.ticker as cticker

def contour_map(fig,img_extent,spec1,spec2):
    fig.set_extent(img_extent, crs=ccrs.PlateCarree())
    fig.add_feature(cfeature.COASTLINE.with_scale('50m')) 
    fig.add_feature(cfeature.LAKES, alpha=0.5)
    fig.set_xticks(np.arange(leftlon,rightlon+spec1,spec1), crs=ccrs.PlateCarree())
    fig.set_yticks(np.arange(lowerlat,upperlat+spec2,spec2), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    fig.xaxis.set_major_formatter(lon_formatter)
    fig.yaxis.set_major_formatter(lat_formatter)

该库有几个基本函数是必须掌握的，我们一一介绍。

1
2
3
4
svd = xMCA(a, b).solver()
lp, rp = svd.patterns(n=2)
le, re = svd.expansionCoefs(n=2)
frac = svd.covFracs(n=2)

svd = xMCA(a, b).solver()建立一个SVD分解器，a，b为要进行分解的两个变量，

svd.patterns(n=2)，svd.expansionCoefs(n=2)，svd.covFracs(n=2)分别取出前两个模态的空间模态，时间系数和方差。lx为左场，rx为右场。

三、示例

我们以北大西洋海温和东亚区域夏季降水的分解为例，展示SVD分析完整的Python实现。

数据为GPCP逐月降水数据和哈德来中心的逐月海温数据。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#读取海温数据（sst格式为时间，纬度，经度）
f_sst = xr.open_dataset('sst.nc')
sst = f_sst['sst'].loc[f_sst.time.dt.month.isin([6,7,8])].loc['1979-06-01':'2018-08-31',60:-30,-100:-0]
sst_lat = sst.coords['latitude']
sst_lon = sst.coords['longitude']
#读取降水数据（pre格式为时间，纬度，经度）
f_pre = xr.open_dataset('precip.mon.total.v2020.nc')
pre = f_pre['precip'].loc[f_pre.time.dt.month.isin([6,7,8])].loc['1979-06-01':'2018-08-31',55:20,70:140]
pre_lat = pre.coords['lat']
pre_lon = pre.coords['lon']
#数据标准化，添补缺测值，并处理为DataArray格式
time = np.arange(1979,2019,1)
sst_summer = np.array(sst).reshape((40,3,sst_lat.shape[0],sst_lon.shape[0])).mean((1))
pre_summer = np.array(pre).reshape((40,3,pre_lat.shape[0],pre_lon.shape[0])).mean((1))
pre_summer[np.isnan(pre_summer)] = 0
sst_summer = (sst_summer - np.mean(sst_summer,axis = 0)) / np.std(sst_summer,axis = 0)
pre_summer = (pre_summer - np.mean(pre_summer,axis = 0)) / np.std(pre_summer,axis = 0)
pre_summer[np.isnan(pre_summer)] = 0
sst_summer = xr.DataArray(sst_summer, coords=[time,sst_lat,sst_lon], dims=['time', 'lat', 'lon'])
pre_summer = xr.DataArray(pre_summer, coords=[time,pre_lat,pre_lon], dims=['time', 'lat', 'lon'])
#SVD分解
svd = xMCA(sst_summer, pre_summer)
svd.solver()
lp, rp = svd.patterns(n=2)
le, re = svd.expansionCoefs(n=2)
frac = svd.covFracs(n=2)

以下是绘图部分，我们先给出其中不重复的部分，文章最末会给出完整代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
fig1 = plt.figure(figsize=(12,8))

f1_ax1 = fig1.add_axes([0.1, 0.6, 0.5, 0.5],projection = ccrs.PlateCarree(central_longitude=-50))
leftlon, rightlon, lowerlat, upperlat = (-100,0,-30,60)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
contour_map(f1_ax1,img_extent,25,30)
f1_ax1.set_title('1st pattern of left',loc='center',fontsize=18)
f1_ax1.set_title('(a)',loc='left',fontsize=18)
c1 = f1_ax1.contourf(sst_lon,sst_lat, lp[0], zorder=0, extend = 'both', levels =np.arange(-1,1.1,0.1),transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)

f1_ax2 = fig1.add_axes([0.15, 0.1, 0.4, 0.4],projection = ccrs.PlateCarree(central_longitude=105))
leftlon, rightlon, lowerlat, upperlat = (70,140,20,55)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
b = contour_map(f1_ax2,img_extent,35,7)
f1_ax2.set_title('1st pattern of right',loc='center',fontsize=18)
f1_ax2.set_title('(c)',loc='left',fontsize=18)
c2 = f1_ax2.contourf(pre_lon,pre_lat, rp[0], zorder=0, extend = 'both', levels =np.arange(-1,1.1,0.1), transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)

position=fig1.add_axes([0.16, 0.085, 0.35, 0.025])
fig1.colorbar(c2,cax=position,orientation='horizontal',format='%.1f',)

f1_ax3 = fig1.add_axes([0.6, 0.6, 0.4, 0.5])
f1_ax3.plot(time,le[0])
f1_ax3.set_title('(b)',loc='left',fontsize=18)
f1_ax3.set_title('1st time series of left',loc='center',fontsize=18)


f1_ax4 = fig1.add_axes([0.6, 0.15, 0.4, 0.3])
f1_ax4.plot(time,re[0])
f1_ax4.set_title('(d)',loc='left',fontsize=18)
f1_ax4.set_title('1st time series of right',loc='center',fontsize=18)