4、EOF分析方法

​

EOF简介

EOF(经验正交函数分解方法）是大气科学以及海洋科学中常用的分析方法，本文以海温数据为例讲解如何在Python中快速实现资料的EOF处理，关于EOF原理及分析方法不再赘述，感兴趣的小伙伴可以参考南京信息工程大学李丽平老师的课件EOF分析原理。

案例讲解

表征海洋的多个变量之间是相互影响的，这些变量之间的关系是复杂的。利用EOF（即经验正交函数Empirical Orthogonal Function）分析方法，可以实现用少数变量反映原来多变量体现的信息，从而有效地提取物理场的主要信息，并排除随机的干扰信息。EOF方法能够比较客观地反映要素场的主要特征，而且不受数据的空间站点和区域范围等条件的限制，适合中、大尺度分析。
本文将对来自于HadISST的多年海表温度（SST）数据进行EOF分析练习，从而掌握NetCDF文件的读取、相关的数据预处理以及EOF分析的基本原理和实现方法。另外，在分析结果的可视化方面，除了使用Matplotlib绘图以外，还会使用Basemap或Cartopy进行地图相关的绘制。

海温数据信息查看

0.1导入软件包

import xarray as xr
import numpy as np
import datetime as dt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.font_manager import FontProperties   
from eofs.standard import Eof
import numpy.ma as ma
import warnings
warnings.filterwarnings("ignore")

0.2读取海温数据

#读取海温数据并转换为np.array
with xr.open_dataset("./HadISST_npac.nc") as f:
    sst=f.sst
#查看数据信息
sst

0.3查看缺测值

num_nan=np.isnan(sst).sum()
print(num_nan) #可以看到sst数据中存在7135488个缺测数据
print(num_nan/1824) #总缺测数据个数/时间个数=缺测格点数，可以看到共有3912个缺测格点

0.4查看异常值

除了缺测值外，数据中的异常值也会对分析结果产生影响，因此数据处理前查看是否存在异常值也是十分重要的

sst_arr=np.array(sst)#将sst转化为ndarray
sst_arr=sst_arr.reshape(1,-1)#将sst_arr转化成1*26356800的二维数组
arr_u=np.unique(sst_arr)
print("去重后的数组为:\n",arr_u)#可以看到海温数据最小值为-1000，最大值为32.8，其中-1000为异常值，所以在分析时需要进行处理

1 修补异常值并去除缺测值后进行EOF分析

经过上述分析可以知道，数据中存在缺测值与异常值，所以接下来将对异常值进行修补并去除缺测值后进行EOF分解，查看效果

1.1 读取数据

import xarray as xr
import numpy as np
import datetime as dt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.font_manager import FontProperties   
from eofs.standard import Eof
import numpy.ma as ma
import warnings
warnings.filterwarnings("ignore")

f=xr.open_dataset("./HadISST_npac.nc")
f

sst=np.array(f.sst)
sst.shape

1.2异常数据修补（利用周围格点均值代替）

处理思路(以区域1为例）：
（1）首先，判断sst[time,lat,lon]格点是否存在异常值，如果存在则进行第（2）步，如果为正常值则直接跳过。
（2）sst[time,lat,lon]为异常值时，将其作为左上角点，判断其周围边长为i个格点的数据范围sst[time,lat:lat+i,lon:lon+i]是否存在正常值，如存在，则进行第（3）步,否则扩大数据范围，继续进行判断。
（3）使用sst[time,lat:lat+i,lon:lon+i]的平均值修补sst[time,lat,lon]。
说明：将全球范围分成4个区域是为了便于对异常值进行修补，其他区域修补方法见代码，思想都是一样

#区域1数据修补
num=43
for time in np.arange(1824):
    for lat in np.arange(0,43):
        for lon in np.arange(0,85):
            if sst[time,lat,lon]==-1000.0:
                for i in np.arange(2,num):
                    data=sst[time,lat:lat+i,lon:lon+i]
                    temp=data
                    a=np.isnan(data).sum()
                    b=np.sum(data==-1000)
                    if (a+b)==i**2:
                        pass
                    else:
                        data=data[~np.isnan(data)]
                        data=data[data!=-1000].mean()
                        sst[time,lat,lon]=data
                        break

#区域2数据修补
for time in np.arange(1824):
    for lat in np.arange(0,43):
        for lon in np.arange(85,170):
            if sst[time,lat,lon]==-1000:
                for i in np.arange(2,num):
                    data=sst[time,lat:lat+i,lon-i:lon]
                    a=np.isnan(data).sum()
                    b=np.sum(data==-1000)
                    if (a+b)==i**2:
                        pass
                    else:
                        data=data[~np.isnan(data)]
                        data=data[data!=-1000].mean()
                        sst[time,lat,lon]=data
                        break
#区域3数据修补
for time in np.arange(1824):
    for lat in np.arange(43,85):
        for lon in np.arange(0,85):
            if sst[time,lat,lon]==-1000:
                for i in np.arange(2,num):
                    data=sst[time,lat-i:lat,lon:lon+i]
                    a=np.isnan(data).sum()
                    b=np.sum(data==-1000)
                    if (a+b)==i**2:
                        pass
                    else:
                        data=data[~np.isnan(data)]
                        data=data[data!=-1000].mean()
                        sst[time,lat,lon]=data
                        break

#区域4数据修补
for time in np.arange(1824):
    for lat in np.arange(43,85):
        for lon in np.arange(85,170):
            if sst[time,lat,lon]==-1000:
                for i in np.arange(2,num):
                    data=sst[time,lat-i:lat,lon-i:lon]
                    a=np.isnan(data).sum()
                    b=np.sum(data==-1000)
                    if (a+b)==i**2:
                        pass
                    else:
                        data=data[~np.isnan(data)]
                        data=data[data!=-1000].mean()
                        sst[time,lat,lon]=data
                        break

1.3 去除季节趋势

#数据标准化之前需要去除季节变化对数据的影响
sst=sst.reshape(152,12,85,170)
climate_sst=np.empty((12,85,170))
for month in np.arange(12):
    for lat in np.arange(85):
        for lon in np.arange(170):
            climate_sst[month,lat,lon]=sst[:,month,lat,lon].mean()
for month in np.arange(12):
    for lat in np.arange(85):
        for lon in np.arange(170):
            sst[:,month,lat,lon]=sst[:,month,lat,lon]-climate_sst[month,lat,lon]
sst=sst.reshape(152*12,85,170)

1.4 EOF分解（不去除缺测值）

从eofs.standard参考文档中可以看到针对规则的缺测值，即由于陆地造成的缺测，eof.standard在分解时可以自动进行处理，无需手工处理。

lat = f.latitude
lon = f.longitude
lat = np.array(lat)
coslat = np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
solver1 = Eof(sst,center=True,weights=wgts)
eof1 = solver1.eofsAsCorrelation(neofs=3)
pc1 = solver1.pcs(npcs=3, pcscaling=1)
var1 = solver1.varianceFraction()
print("var1:\n",var1[:3])

1.5 EOF分解（去除缺测值）

虽然eof.standard模块可以对规则的缺测值进行处理，但同时本文也对缺测数据进行了删除处理，并比较两种不同数据处理方法结果是否存在较大差别

sst=sst[~np.isnan(sst)]#去除缺测值
sst=sst.reshape(1824,-1)#将sst转化为（1824，10538）的二维数组
sst.shape

lat = f.latitude
lon = f.longitude
lat = np.array(lat)
coslat = np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
solver2 = Eof(sst,center=True,weights=wgts)
eof2 = solver2.eofsAsCorrelation(neofs=3)
pc2 = solver2.pcs(npcs=3, pcscaling=1)
var2 = solver2.varianceFraction()
print("var2:\n",var2[:3])

1.6 结果分析

对比方差贡献var1(不去除缺测值时方差变量）与var2(去除缺测值时方差变量）可以看到两者差别并不大，说明缺测值的去除与否并不会对EOF分析结果产生大的影响。
同时对比参考答案，可以看到虽然本文结果与参考答案并非完全相同，但相差并不大，这可能是由于异常数据处理方法不同所导致的

print("题中参考答案:",[0.3257,0.0941,0.0824])
print("var1:",list(var1[:3]))
print("var2:",list(var2[:3]))

2 可视化与相关分析

2.1 EOF结果分析

def mapart(ax):
    '''添加地图元素'''
    ax.coastlines(color="k",lw=0.5)
    ax.add_feature(cfeature.LAND,facecolor="white")
    #设置地图范围
    ax.set_extent([105.5,275.5,-15,70],crs=ccrs.PlateCarree())#[leftlon,rightlon,lowerlat,upperlat]
    #设置经纬度标签
    ax.set_xticks([105,135,165,195,225,255],crs=ccrs.PlateCarree())
    #标注的位置是根据cartopy系统内部位置而定的，与数据本身无关
    ax.set_yticks([-15,-5,5,15,25,35,45,55,65],crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)


def eof_contourf(eofs,pcs,pers):
    
    plt.close
    fig=plt.figure(figsize=(18,10))
    projection=ccrs.PlateCarree(central_longitude=-180)
    year=range(10538)
    ax1=fig.add_subplot(3,2,1,projection=projection)
    mapart(ax1)
    lon=np.arange(105.5,275.5)
    lat=np.arange(69.5,-15.5,-1)
    p=ax1.contourf(lon,lat,eofs[0,:,:],transform=ccrs.PlateCarree())
    ax1.set_title('mode1 (%s'%(round(pers[0],2))+"%)",loc ='left')
    time=np.arange(1824)
    ax2=fig.add_subplot(3,2,2)
    ax2.plot(time,pcs[:,0],color="k",linewidth=2,linestyle="--")
    ax2.set_xticks(list(np.arange(120,1824,240)),[1880,1900,1920,1940,1960,1980,2000,2020])
    ax2.set_title('pc1 (%s'%(round(pers[0],2))+"%)",loc ='left')
    
    ax3 = fig.add_subplot(3,2, 3, projection=projection)
    mapart(ax3)
    pp = ax3.contourf(lon,lat,eofs[1,:,:],transform=ccrs.PlateCarree())
    ax3.set_title('mode2 (%s'%(round(pers[1],2))+"%)",loc ='left')
    
    ax4 = fig.add_subplot(3,2, 4)
    ax4.plot(time,pcs[:,1] ,color='k',linewidth=1.2,linestyle='--')
    ax4.set_xticks(list(np.arange(120,1824,240)),[1880,1900,1920,1940,1960,1980,2000,2020])
    ax4.set_title('pc2 (%s'%(round(pers[1],2))+"%)",loc ='left')
    
    ax5 = fig.add_subplot(3,2, 5, projection=projection)
    mapart(ax5)
    ppp = ax5.contourf(lon,lat,eofs[2,:,:],transform=ccrs.PlateCarree())
    ax5.set_title('mode3 (%s'%(round(pers[2],2))+"%)",loc ='left')
    
    ax6 = fig.add_subplot(3,2, 6)
    ax6.plot(time,pcs[:,2])
    ax6.set_xticks(list(np.arange(120,1824,240)),[1880,1900,1920,1940,1960,1980,2000,2020])
    ax6.set_title('pc3 (%s'%(round(pers[2],2))+"%)",loc ='left')
    
    #添加0线
    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='-')
    
    #在图下边留白边放colorbar        
    fig.subplots_adjust(bottom=0.1)
    #colorbar位置： 左 下 宽 高 
    l = 0.25
    b = 0.04
    w = 0.6
    h = 0.015
    #对应 l,b,w,h；设置colorbar位置；
    rect = [l,b,w,h] 
    cbar_ax = fig.add_axes(rect) 
    
    c=plt.colorbar(pp, cax=cbar_ax,orientation='horizontal', aspect=20, pad=0.1)
    c.ax.tick_params(labelsize=6)
    # c.set_label('%s'%(labelname),fontsize=20)#手动设置色标
    # c.set_ticks(np.arange(1,6,1))
    plt.subplots_adjust( wspace=-0.12,hspace=0.3)

    
    plt.savefig('eof.tif',dpi=600,bbox_inches = 'tight',transparent=True, pad_inches = 0)
    plt.show()    

eof_contourf(eof1,pc1,var1[:3]*100)

2.2 相关分析

#由于本人对dat格式的数据并不熟练，所以首先将数据导出成了excel文件，然后利用excel软件进行处理,最后将文件（PDO.xlsx)上传，利用程序读取
import pandas as pd
df=pd.read_excel("./PDO_index.xlsx")
df.head()
df=df.iloc[16:-1,:]
df.drop(["Year"],axis=1,inplace=True)
pdo_index=np.array(df).reshape(152*12,-1)
pdo_index.shape

#计算相关系数
corr = np.corrcoef(pc1[:,0],pdo_index[:].reshape(1824))
corr

import matplotlib.pyplot as plt
time=np.arange(1824)
fig,ax=plt.subplots(figsize=(12,6))
ax.plot(time,pc1[:,0],color="red",linewidth=2,linestyle="--")
ax.plot(time,pdo_index[:],color="blue",linewidth=2,linestyle="-")
ax.set_xticks(list(np.arange(120,1824,240)),[1880,1900,1920,1940,1960,1980,2000,2020])
ax.set_title("R=0.548",loc ='left')

数据链接

百度网盘数据链接