本文的jupyter notebook已上傳至github，歡迎fork測(cè)試~

0.導(dǎo)言

我們都知道，氣候研究的時(shí)間跨度一般都較長，基本都在30年以上，這就意味著對(duì)應(yīng)的數(shù)據(jù)集十分龐大，既不能簡(jiǎn)單地對(duì)數(shù)據(jù)進(jìn)行描述，也無法輕易地從數(shù)據(jù)中提取特征。那么面對(duì)如此龐大的數(shù)據(jù)集，我們?nèi)绾尾拍軓闹刑崛〕鲎钅苊枋鏊闹饕卣髂?

在這種場(chǎng)景下，EOF分解就顯示出它的強(qiáng)大優(yōu)勢(shì)了。它可以把隨著時(shí)間變化的氣象要素場(chǎng)，分解為空間函數(shù)部分和時(shí)間函數(shù)部分，從而便于我們開展分析和研究，讓我們能夠從龐大的氣候數(shù)據(jù)中抓住他們的主要特征。本文將介紹EOF的基本原理，以及在python中如何實(shí)現(xiàn)EOF分解。如有不妥之處，還請(qǐng)大家不吝賜教!

1.什么是經(jīng)驗(yàn)正交分解

經(jīng)驗(yàn)正交分解(Empirical orthogonal function，EOF，以下簡(jiǎn)稱“EOF分解”)又被稱為主成分分析(Principal component analysis，PCA)，最早由Pearson提出，20世紀(jì)50年代Lorenz將該方法引入大氣科學(xué)領(lǐng)域，隨后被廣泛應(yīng)用至今。

EOF分解是一種分析矩陣數(shù)據(jù)中的結(jié)構(gòu)特征、提取主要數(shù)據(jù)特征量的方法。

EOF分解的分析對(duì)象，就是我們?cè)跉庀罂蒲泻蜆I(yè)務(wù)工作中，經(jīng)常要分析各種氣象要素場(chǎng)，如海表溫度場(chǎng)、氣壓場(chǎng)、降水場(chǎng)等。這些氣象要素場(chǎng)大多由不規(guī)則的網(wǎng)格點(diǎn)所組成。如果我們抽取這些場(chǎng)的某一段歷史時(shí)期的資料，就構(gòu)成一組以網(wǎng)格點(diǎn)為空間點(diǎn)(可看作多個(gè)變量)的隨時(shí)間變化的樣本。

EOF分解所做的，就是把隨時(shí)間變化的氣象要素場(chǎng)分解為空間函數(shù)部分和時(shí)間函數(shù)部分，且這兩部分相互正交。

其中，空間函數(shù)部分可在一定程度上概括要素場(chǎng)的空間分布特點(diǎn)，這部分是不隨時(shí)間變化的；而時(shí)間函數(shù)部分則由空間點(diǎn)(即多個(gè)變量)的線性組合所構(gòu)成，稱為主分量。

由于這些主分量中前幾個(gè)就可能占有原空間點(diǎn)(即多個(gè)變量)總方差的很大部分，即通過EOF分解可以很容易地將原始要素場(chǎng)的變化信息濃縮在前幾個(gè)主分量及其對(duì)應(yīng)的空間函數(shù)上，所以，EOF分解常被用于氣象要素場(chǎng)時(shí)空變化特征規(guī)律的提取和研究。

2.基本原理

在氣象科研和業(yè)務(wù)工作中，我們經(jīng)常會(huì)分析某一氣象要素的時(shí)空數(shù)據(jù)集，它們大多是3維的，包括2維的空間場(chǎng)以及1維的時(shí)間序列。其中2維空間場(chǎng)對(duì)應(yīng)著一個(gè)個(gè)的空間點(diǎn)，1維的時(shí)間序列對(duì)應(yīng)著觀測(cè)的時(shí)間點(diǎn)。

現(xiàn)在假設(shè)我們要對(duì)某一個(gè)氣象要素的時(shí)空數(shù)據(jù)集進(jìn)行研究，它共有 m 個(gè)空間點(diǎn) ( $x_1,x_2,···,x_m$ ), n 個(gè)時(shí)間點(diǎn) ( $t_1,t_2,···,t_n$ )。那對(duì)于這樣一個(gè)時(shí)空數(shù)據(jù)集，我們就可以用一個(gè) $m×n$ 的矩陣 $F$ 來表示。

其中 $F$ 的行對(duì)應(yīng)著某個(gè)空間點(diǎn)在觀測(cè)期內(nèi)所有時(shí)間點(diǎn)上的值， $F$ 的列對(duì)應(yīng)著某個(gè)時(shí)間點(diǎn)在觀測(cè)區(qū)域內(nèi)所有空間點(diǎn)上的值。

$F= \begin{pmatrix} f_{11} & f_{12} & ··· & f_{1n}\\ f_{21} & f_{22} & ··· & f_{2n}\\ ··· & ··· & ··· & ···\\ f_{m1} & f_{m2} & ··· & f_{mn}\\ \end{pmatrix}$

而我們的EOF分解方法，就是將時(shí)空數(shù)據(jù)集 $F$ 分解為空間函數(shù) $V$ 和時(shí)間函數(shù) $Z$ 兩部分，即

$F=VZ$

$\begin{pmatrix} f_{11} & f_{12} & ··· & f_{1n}\\ f_{21} & f_{22} & ··· & f_{2n}\\ ··· & ··· & ··· & ···\\ f_{m1} & f_{m2} & ··· & f_{mn}\\ \end{pmatrix} = \begin{pmatrix} v_{11} & v_{12} & ··· & v_{1m}\\ v_{21} & v_{22} & ··· & v_{2m}\\ ··· & ··· & ··· & ···\\ v_{m1} & v_{m2} & ··· & v_{mm}\\ \end{pmatrix} \begin{pmatrix} z_{11} & z_{12} & ··· & z_{1n}\\ z_{21} & z_{22} & ··· & z_{2n}\\ ··· & ··· & ··· & ···\\ z_{m1} & z_{m2} & ··· & z_{mn}\\ \end{pmatrix}$

也可表示為

$F_{it}=\sum_{k=1}^{m}v_{ik}z_{kt}=v_{i1}z_{1t}+v_{i2}z_{2t}+···+v_{im}z_{mt}$

其中下標(biāo) $i=1,2,···,m$ 表示空間點(diǎn)，下標(biāo) $t=1,2,···,n$ 表示時(shí)間點(diǎn)。

因此，該分解式的含義為：氣象要素場(chǎng)中第 $i$ 個(gè)空間點(diǎn)上的第 $t$ 次觀測(cè)值，可以看作是 m 個(gè)空間函數(shù) $v_{i1},v_{i2},···,v_{im}$ 和 m 個(gè)時(shí)間函數(shù) $z_{1t},z_{2t},···,z_{mt}$ 的線性組合。

在空間函數(shù) $V$ 中，它的每一列(如第k列，可表示為 $\mathbf{v_k}=(v_{1k},v_{2k},···,v_{mk})$ 都可以看作為一個(gè)典型場(chǎng)，稱為一個(gè)空間模態(tài)(EOF)，它只跟空間點(diǎn)有關(guān)，而與時(shí)間無關(guān)。

在時(shí)間函數(shù) $Z$ 中，它的每一行(如第k行，可表示為 $\mathbf{z_k}=(z_{k1},z_{k2},···,z_{kn})$ 都可以看作是一個(gè)時(shí)間序列，稱為時(shí)間系數(shù)(PC)。

在此基礎(chǔ)上，時(shí)空數(shù)據(jù)集 $F$ 中第 t 時(shí)刻的空間場(chǎng)(共有m個(gè)空間點(diǎn))，就可以表示為

$$
\begin{pmatrix}
f_{1t} \
f_{2t} \
··· \
f_{mt} \
\end{pmatrix}
=
\begin{pmatrix}
v_{11} \
v_{21} \
··· \
v_{m1} \
\end{pmatrix}
z_{1t}

\begin{pmatrix}
v_{12} \
v_{22} \
··· \
v_{m2} \
\end{pmatrix}
z_{2t}

···

\begin{pmatrix}
v_{1m} \
v_{2m} \
··· \
v_{mm} \
\end{pmatrix}
z_{mt}
$$

該式表明，時(shí)空數(shù)據(jù)集中第 t 時(shí)刻的空間場(chǎng)可以表示為 m 個(gè)空間模態(tài)( $\mathbf{v_1},\mathbf{v_2},···,\mathbf{v_m}$ )，按照不同的權(quán)重線性 ( $z_{1t},z_{2t},···,z_{mt}$ ) 疊加而成。

對(duì)于任意一個(gè)空間模態(tài) $\mathbf{v_k}$ 而言，其在不同時(shí)刻對(duì)應(yīng)的權(quán)重有所不同，若將其在全部觀測(cè)時(shí)間點(diǎn)所對(duì)應(yīng)的權(quán)重 $(z_{k1},z_{k2},···,z_{kn})$ 列在一起，就是我們上面提到的時(shí)間系數(shù) $\mathbf{z_k}$ 。

每一個(gè)空間模態(tài)都對(duì)應(yīng)著一個(gè)自己的時(shí)間系數(shù) $\mathbf{z_k}$ ，如果 $\mathbf{z_k}$ 中的正值的絕對(duì)值越大，表明該時(shí)間點(diǎn)所對(duì)應(yīng)的數(shù)據(jù)集與該時(shí)間系數(shù)所對(duì)應(yīng)的空間模態(tài)相似度越高；如果 $\mathbf{z_k}$ 中的負(fù)值的絕對(duì)值越大，則表明該時(shí)間點(diǎn)所對(duì)應(yīng)的數(shù)據(jù)集越傾向于和該時(shí)間系數(shù)所對(duì)應(yīng)的空間模態(tài)呈相反的狀態(tài)。

除此以外，我們還可以根據(jù)每一個(gè)空間模態(tài)(或時(shí)間系數(shù))計(jì)算該空間模態(tài)(或時(shí)間系數(shù))對(duì)應(yīng)的方差貢獻(xiàn)率，方差貢獻(xiàn)率越高，表明該空間模態(tài)(或時(shí)間系數(shù))所包含的關(guān)于原始數(shù)據(jù)的變化信息越多。

我們通常根據(jù)方差貢獻(xiàn)率的大小對(duì)空間模態(tài)進(jìn)行排序，方差貢獻(xiàn)率最大的稱為第一空間模態(tài)，方差貢獻(xiàn)率第二大的稱為第二空間模態(tài)，以此類推。

因?yàn)樗锌臻g模態(tài)(或時(shí)間系數(shù))所對(duì)應(yīng)的方差貢獻(xiàn)率加起來等于 1，所以，當(dāng)方差貢獻(xiàn)率最高的幾個(gè)空間模態(tài)(或時(shí)間系數(shù))加起來得到的累計(jì)方差貢獻(xiàn)率較大時(shí)(比如達(dá)到80%)，我們就會(huì)用它們幾個(gè)來表示原始數(shù)據(jù)中所包含的主要特征。

但需要強(qiáng)調(diào)的是，EOF分解只是一種統(tǒng)計(jì)方法，其得到的分解結(jié)果本身沒有任何的物理意義，需要我們結(jié)合實(shí)際及專業(yè)知識(shí)去對(duì)結(jié)果進(jìn)行理解和解釋的。

3.在Python實(shí)現(xiàn)EOF分解

import struct
import numpy as np
import xarray as xr
from eofs.xarray import Eof

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
import cartopy.crs as ccrs       
import cartopy.feature as cfeature 
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter, LatitudeLocator, LongitudeLocator)
import cmaps
import warnings
warnings.filterwarnings('ignore')

3.1進(jìn)行EOF分析

# 讀取數(shù)據(jù)

data = np.fromfile('ssta.dat', dtype='float32')
data = np.reshape(data,(56, 11, 336), order='F')
lon = np.arange(160, 272, 2) #經(jīng)度160E~90W（間隔2度)
lat = np.arange(-10, 12, 2)  #緯度10S~10N（間隔2度） 
time = xr.cftime_range(start='1979-01', end='2006-12', freq='MS') #時(shí)間1979.1~2006.12（逐月)
ds = xr.Dataset({'ssta': ([ 'lon', 'lat', 'time'], data)},
                coords={'lon': lon, 'lat': lat, 'time': time})

#檢驗(yàn)數(shù)據(jù)讀取
fig = plt.figure(figsize=(18, 12), dpi=300)  
projection = ccrs.PlateCarree(central_longitude=180) #指定投影為經(jīng)緯度投影，中心經(jīng)緯度為180°
ax1 = fig.add_subplot(1, 1, 1, projection=projection) #繪制第一個(gè)子圖
projection = ccrs.PlateCarree(central_longitude=180) #指定投影為經(jīng)緯度投影，并指定中心經(jīng)度為180°
# 設(shè)置地圖范圍，經(jīng)度為（160, 270)，緯度為（-10, 10）
ax1.set_extent([-20, 90, -10, 10], crs=ccrs.PlateCarree(central_longitude=180))
# 設(shè)置經(jīng)緯度標(biāo)簽
ax1.set_xticks([-20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90], crs=projection)
ax1.set_yticks([-10, -5, 0, 5, 10], crs=projection)
lon_formatter = LongitudeFormatter(zero_direction_label=True) #給標(biāo)簽添加對(duì)應(yīng)的N,S,E,W
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
ax1.contourf(ds.lon,ds.lat,ds.ssta[:,:,0].transpose('lat', 'lon'), np.linspace(-1,1,51),cmap=cmaps.BlueWhiteOrangeRed,transform=ccrs.PlateCarree())  #在第一個(gè)地圖上繪制第一空間模態(tài)的填充色圖，EOFs[0]表示從EOF中取出第一個(gè)空間模態(tài)
ax1.coastlines(color='k', lw=0.5)                #添加海岸線
ax1.add_feature(cfeature.LAND, facecolor='white') #添加陸地

# 轉(zhuǎn)置數(shù)據(jù)
ds = ds.transpose('time', 'lat', 'lon')

#計(jì)算網(wǎng)格點(diǎn)的權(quán)重
coslat = np.cos(np.deg2rad(ds.coords['lat'].values)) 
wgts = np.sqrt(coslat)[..., np.newaxis]

# 計(jì)算EOF分解
solver = Eof(ds['ssta'], weights=wgts, center=False)
eofs = solver.eofs(neofs=10, eofscaling=2)
pcs = solver.pcs(npcs=10, pcscaling=1)
variance = solver.varianceFraction(neigs=10)

# 繪制方差解釋率圖
plt.figure(figsize=(8, 6),dpi=200)
plt.plot(np.arange(1, 11), variance[:10]*100, 'o-')
plt.xlabel('EOF modes')
plt.ylabel('Variance explained (%)')
plt.title('Variance explained by EOF modes')
plt.grid()
# 在每個(gè)點(diǎn)上添加文本標(biāo)簽
for i, v in enumerate(variance.values[:10]*100):
    plt.text(i+1, v+0.9, f'{v:.2f}', fontsize=8, color='black')
plt.show()

#計(jì)算特征值和NorthTest的結(jié)果
eigenvalues = solver.eigenvalues(neigs=10)
errors = solver.northTest(neigs=10)
# 繪制Error Bar圖
plt.figure(figsize=(8, 6),dpi=200)
plt.errorbar(np.arange(1, 11), eigenvalues, yerr=errors, fmt='o-', capsize=5)
plt.xlabel('EOF modes')
plt.ylabel('Eigenvalues')
plt.title('NorthTest(with error bars)')
plt.grid()
plt.show()

3.2給出前4個(gè)模態(tài)的空間分布和時(shí)間序列（前4個(gè)modes的誤差棒無相重疊部分）

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
import cmaps

def mapart(ax):
    '''
    添加地圖元素
    '''
    # 指定投影為經(jīng)緯度投影，并指定中心經(jīng)度為180°
    projection = ccrs.PlateCarree(central_longitude=180)
    # 設(shè)置地圖范圍，經(jīng)度為（160, 270)，緯度為（-10, 10）
    ax.set_extent([-20, 90, -10, 10], crs=ccrs.PlateCarree(central_longitude=180))
    # 設(shè)置經(jīng)緯度標(biāo)簽
    ax.set_xticks([-20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90], crs=projection)
    ax.set_yticks([-10, -5, 0, 5, 10], crs=projection)
    # 給標(biāo)簽添加對(duì)應(yīng)的N,S,E,W
    lon_formatter = LongitudeFormatter(zero_direction_label=True)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    # 添加海岸線
    #ax.coastlines(color='k', lw=0.5)
    # 添加陸地
    #ax.add_feature(cfeature.LAND, facecolor='white')

def eof_contourf(EOFs, PCs, pers):
    '''
    繪制EOF填充圖
    '''
    plt.close
    # 將方差轉(zhuǎn)換為百分?jǐn)?shù)的形式，如0.55變?yōu)?5%
    pers=(pers*100).values
    # 指定畫布大小以及像素
    fig = plt.figure(figsize=(18, 12), dpi=300)
    # 指定投影為經(jīng)緯度投影，中心經(jīng)緯度為180°
    projection = ccrs.PlateCarree(central_longitude=180)
    # 將橫坐標(biāo)轉(zhuǎn)換為日期格式
    dates_num = [mdates.date2num(t) for t in time]
    
    # 繪制第一個(gè)子圖
    ax1 = fig.add_subplot(4, 2, 1, projection=projection)
    # 為第一個(gè)子圖添加地圖底圖
    mapart(ax1)
    # 在第一個(gè)地圖上繪制第一空間模態(tài)的填充色圖，EOFs[0]表示從EOF中取出第一個(gè)空間模態(tài)
    p = ax1.contourf(ds.lon,
    ds.lat,
    EOFs[0],
    np.linspace(-1, 1, 51),
    cmap=cmaps.BlueWhiteOrangeRed,
    transform=ccrs.PlateCarree(),
    extend='both')
    # 為第一個(gè)子圖添加標(biāo)題，其中，pers[0]表示從pers中取出第一個(gè)空間模態(tài)對(duì)應(yīng)的方差貢獻(xiàn)率
    ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left')
    ax1.set_aspect(1.6)

    # 繪制第二個(gè)子圖
    ax2 = fig.add_subplot(4, 2, 2) 
    ax2.set_xticks(dates_num[::108])  # 每9年的1月作為橫坐標(biāo)標(biāo)簽
    ax2.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
    ax2.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108]))  # 每9年的1月作為橫坐標(biāo)刻度線
    b = ax2.bar(dates_num, PCs[:, 0], width=50, color='r')
    # 對(duì)時(shí)間系數(shù)值小于0的柱子設(shè)置為藍(lán)色
    for bar, height in zip(b, PCs[:, 0]):
        if height < 0:
            bar.set(color='blue') 
    # 為第二個(gè)子圖添加標(biāo)題
    ax2.set_title('PC1' % (round(pers[0], 2)), loc='left') 
    


    # 后面ax3和ax5與ax1繪制方法相同，ax4和ax6與ax2繪制方法相同，不再注釋

    ax3 = fig.add_subplot(4, 2, 3, projection=projection) 
    mapart(ax3)
    pp = ax3.contourf(ds.lon,
                      ds.lat,
                      EOFs[1],
                      np.linspace(-1, 1, 51),
                      cmap=cmaps.BlueWhiteOrangeRed,
                      transform=ccrs.PlateCarree(),
                      extend='both')
    ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')
    ax3.set_aspect(1.6)

    ax4 = fig.add_subplot(4, 2, 4)
    ax4.set_xticks(dates_num[::108])  # 每9年的1月作為橫坐標(biāo)標(biāo)簽
    ax4.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
    ax4.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108]))  # 每9年的1月作為橫坐標(biāo)刻度線
    ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')
    bb = ax4.bar(dates_num, PCs[:, 1], width=50, color='r')
    for bar, height in zip(bb, PCs[:, 1]):
        if height < 0:
            bar.set(color='blue')

    ax5 = fig.add_subplot(4, 2, 5, projection=projection)
    mapart(ax5)
    ppp = ax5.contourf(ds.lon,
                       ds.lat,
                       EOFs[2],
                       np.linspace(-1, 1, 51),
                       cmap=cmaps.BlueWhiteOrangeRed,
                       transform=ccrs.PlateCarree(),
                       extend='both')
    ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')
    ax5.set_aspect(1.6)

    ax6 = fig.add_subplot(4, 2, 6)
    ax6.set_xticks(dates_num[::108])  # 每9年的1月作為橫坐標(biāo)標(biāo)簽
    ax6.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
    ax6.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108]))  # 每9年的1月作為橫坐標(biāo)刻度線
    ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')
    bbb = ax6.bar(dates_num, PCs[:, 2], width=50, color='r')
    for bar, height in zip(bbb, PCs[:, 2]):
        if height < 0:
            bar.set(color='blue')
            
    ax7 = fig.add_subplot(4, 2, 7, projection=projection)
    mapart(ax7)
    pppp = ax7.contourf(ds.lon,
    ds.lat,
    EOFs[3],
    np.linspace(-1, 1, 51),
    cmap=cmaps.BlueWhiteOrangeRed,
    transform=ccrs.PlateCarree(),
    extend='both')
    ax7.set_title('mode4 (%s' % (round(pers[3], 2))+"%)", loc='left')
    ax7.set_aspect(1.6)

    ax8 = fig.add_subplot(4, 2, 8) 
    ax8.set_xticks(dates_num[::108])  # 每9年的1月作為橫坐標(biāo)標(biāo)簽
    ax8.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
    ax8.xaxis.set_major_locator(ticker.FixedLocator(dates_num[::108]))  # 每9年的1月作為橫坐標(biāo)刻度線
    #ax8.set_xlim(dates_num[0], dates_num[-1])
    bbbb = ax8.bar(dates_num, PCs[:, 3], width=50, color='r') 
    for bar, height in zip(bbbb, PCs[:, 3]):
        if height < 0:
            bar.set(color='blue') 
    ax8.set_title('PC4' % (round(pers[3], 2)), loc='left') 



    # 為柱狀圖添加0標(biāo)準(zhǔn)線
    ax2.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax4.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax6.axhline(y=0,  linewidth=1, color='k', linestyle='-')
    ax8.axhline(y=0,  linewidth=1, color='k', linestyle='-')

    # 在圖下邊留白邊放colorbar
    fig.subplots_adjust(bottom=0.1)
    # colorbar位置： 左 下 寬 高
    l = 0.25
    b = 0.04
    w = 0.6
    h = 0.015
    # 對(duì)應(yīng) l,b,w,h；設(shè)置colorbar位置；
    rect = [l, b, w, h]
    cbar_ax = fig.add_axes(rect)

    # 繪制colorbar
    c = plt.colorbar(
        p,
        cax=cbar_ax,
        orientation='horizontal',
        aspect=20,
        pad=0.1) 

    # 設(shè)置colorbar的標(biāo)簽大小
    c.ax.tick_params(labelsize=14) 

    # 調(diào)整子圖之間的間距
    plt.subplots_adjust(wspace=0.1, hspace=0.3)
    plt.show()
    
eof_contourf(eofs[:4]*(-1),pcs[:,:4]*(-1),variance[:4])