python实现春季海温的EOF分解并进行north检验
作者：第八星系-欣妹儿
邮箱：3035245582@qq.com
备注：此代码计算了海温距平、海温标准化后的分解情况。此篇代码主要参考了CSDN博主森屿星球。
导入函数
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from eofs.standard import Eof
import xarray as xr
from cartopy.mpl import ticker
import pandas as pd
da = xr.open_dataset("D:\\PAPER\\date\\HadISST_sst.nc", engine="netcdf4")数据处理
# 确定选择的时间范围
spring_months =
# 加载数据到内存中
SST.load()
# 选择所需子集数据
SST_spring = SST.sel(latitude=slice(30, -30), longitude=slice(-180, 180), time=slice("1979", "2021")).sel(
time=SST['time.month'].isin(spring_months))
# 提取所需的数据变量
sst1 = SST_spring.sst[:]
sst2 = np.array(sst1)
lat = SST_spring.latitude[:]
lon = SST_spring.longitude[:]
#因为我们需要分析的是海温异常，就需要分析与平均值不同的异常，所以需要将原来的数据和平均值做一个差值，
#通过差值的大小来判断海温异常的趋势以及分布。
sst=np.array(sst1)
ano=sst1.groupby('time.month')-sst1.groupby('time.month').mean('time', skipna=True)##距平
ano1=np.array(ano)
##标准化 std=标准差
ano1_normol = (sst1.groupby('time.month') - sst1.groupby('time.month').mean('time', skipna=True)) / np.std(sst1,axis=0)
#ano1_normol = (sst1 - sst1_spring_mean) / np.std(sst1,axis=0)
ano2=np.array(ano1_normol)
print(sst.shape)
#计算纬度权重
lat=np.array(lat)
coslat=np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]创建EOF分解器
solver=Eof(sst,weights=wgts)###sst为原来的数据；ano1为距平；ano2为标准化
eof=solver.eofsAsCorrelation(neofs=3)
#此处的neofs的值是我们需要的空间模态数，比如这里我们打算画3个模态
pc=solver.pcs(npcs=3,pcscaling=1)#方差
var=solver.varianceFraction(neigs=3)画图-EOF
##绘图
fig=plt.figure(figsize=(15,15))
proj=ccrs.PlateCarree(central_longitude=180)
leftlon,rightlon,lowerlat,upperlat=(-180,180,-30,30)
lon_formatter=ticker.LongitudeFormatter()
lat_formatter=ticker.LatitudeFormatter()
# 绘制第一模态
fig_ax1=fig.add_axes(,projection=proj)
fig_ax1.set_extent(,crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax1.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax1.add_feature(cfeature.COASTLINE,lw=1)
fig_ax1.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax1.set_title('(a) EOF1(1979-2021)',loc='left',fontsize =15)
fig_ax1.set_title('sst',loc='center',fontsize =18)##给每组图加标题
fig_ax1.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(lon,lat, eof, levels=np.arange(-0.9,1.0,0.1 ), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
# 绘制第二模态
fig_ax2=fig.add_axes(,projection=proj)
fig_ax2.set_extent(,crs=ccrs.PlateCarree())
fig_ax2.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax2.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax2.add_feature(cfeature.COASTLINE,lw=1)
fig_ax2.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax2.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax2.xaxis.set_major_formatter(lon_formatter)
fig_ax2.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax2.set_title('(b) EOF2(1979-2021)',loc='left',fontsize =15)
fig_ax2.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c2=fig_ax2.contourf(lon,lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
# 绘制第三模态
fig_ax3=fig.add_axes(,projection=proj)
fig_ax3.set_extent(,crs=ccrs.PlateCarree())
fig_ax3.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax3.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax3.add_feature(cfeature.COASTLINE,lw=1)
fig_ax3.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax3.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax3.xaxis.set_major_formatter(lon_formatter)
fig_ax3.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax3.set_title('(c) EOF3(1979-2021)',loc='left',fontsize =15)
fig_ax3.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c3=fig_ax3.contourf(lon,lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
cbposition=fig.add_axes()
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f')
画图-时间序列
fig=plt.figure(figsize=(15,15))
fig_ax5=fig.add_axes()
fig_ax5.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax5.set_title('sst',loc='center',fontsize = 15)
fig_ax5.set_ylim(-3.5,3.5)
fig_ax5.axhline(0,linestyle="--")
fig_ax5.plot(np.arange(1979,2022,4/12),pc[:,0],color='blue')
fig_ax6 = fig.add_axes()
fig_ax6.set_title('(d) PC2',loc='left',fontsize = 15)
fig_ax6.set_ylim(-3.5,3.5)
fig_ax6.axhline(0,linestyle="--")
fig_ax6.plot(np.arange(1979,2022,4/12),pc[:,1],color='blue')
fig_ax7 = fig.add_axes()
fig_ax7.set_title('(f) PC3',loc='left',fontsize = 15)
fig_ax7.set_ylim(-3.5,3.5)
fig_ax7.axhline(0,linestyle="--")
fig_ax7.plot(np.arange(1979,2022,4/12),pc[:,2],color='blue')
画图-North检验
##进行north检验
lat=np.array(lat)
coslat=np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver=Eof(sst,weights=wgts)
eof=solver.eofsAsCorrelation(neofs=10)
#此处的neofs的值是我们需要的空间模态数，比如这里我们打算画四个模态
pc=solver.pcs(npcs=10,pcscaling=1)#方差
var=solver.varianceFraction(neigs=10)
##计算特征值和NorthTest结果
eigenvalues = solver.eigenvalues(neigs=10)
errors = solver.northTest(neigs=10)
##绘制NorthTest结果Error Bar
plt.figure()
plt.errorbar(np.arange(1,11),eigenvalues,yerr=errors, fmt='o-', capsize=5 ,color = 'k')
plt.xlabel('EOF modes')
plt.ylabel('Eigenvalues')
plt.title('NorthTest(with error bar)_sst')
plt.show()
完整代码
###导入函数
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from eofs.standard import Eof
import xarray as xr
from cartopy.mpl import ticker
import pandas as pd
da = xr.open_dataset("D:\\PAPER\\date\\HadISST_sst.nc", engine="netcdf4")
# 确定选择的时间范围
spring_months =
# 加载数据到内存中
SST.load()
# 选择所需子集数据
SST_spring = SST.sel(latitude=slice(30, -30), longitude=slice(-180, 180), time=slice("1979", "2021")).sel(
time=SST['time.month'].isin(spring_months))
# 提取所需的数据变量
sst1 = SST_spring.sst[:]
sst2 = np.array(sst1)
lat = SST_spring.latitude[:]
lon = SST_spring.longitude[:]
#因为我们需要分析的是海温异常，就需要分析与平均值不同的异常，所以需要将原来的数据和平均值做一个差值，
#通过差值的大小来判断海温异常的趋势以及分布。
sst=np.array(sst1)
ano=sst1.groupby('time.month')-sst1.groupby('time.month').mean('time', skipna=True)##距平
ano1=np.array(ano)
##标准化 std=标准差
ano1_normol = (sst1.groupby('time.month') - sst1.groupby('time.month').mean('time', skipna=True)) / np.std(sst1,axis=0)
#ano1_normol = (sst1 - sst1_spring_mean) / np.std(sst1,axis=0)
ano2=np.array(ano1_normol)
print(sst.shape)
#计算纬度权重
lat=np.array(lat)
coslat=np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver=Eof(sst,weights=wgts)###sst为原来的数据；ano1为距平；ano2为标准化
eof=solver.eofsAsCorrelation(neofs=3)
#此处的neofs的值是我们需要的空间模态数，比如这里我们打算画3个模态
pc=solver.pcs(npcs=3,pcscaling=1)#方差
var=solver.varianceFraction(neigs=3)
##绘图
fig=plt.figure(figsize=(15,15))
proj=ccrs.PlateCarree(central_longitude=180)
leftlon,rightlon,lowerlat,upperlat=(-180,180,-30,30)
lon_formatter=ticker.LongitudeFormatter()
lat_formatter=ticker.LatitudeFormatter()
# 绘制第一模态
fig_ax1=fig.add_axes(,projection=proj)
fig_ax1.set_extent(,crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax1.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax1.add_feature(cfeature.COASTLINE,lw=1)
fig_ax1.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax1.set_title('(a) EOF1(1979-2021)',loc='left',fontsize =15)
fig_ax1.set_title('sst',loc='center',fontsize =18)##给每组图加标题
fig_ax1.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(lon,lat, eof, levels=np.arange(-0.9,1.0,0.1 ), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
# 绘制第二模态
fig_ax2=fig.add_axes(,projection=proj)
fig_ax2.set_extent(,crs=ccrs.PlateCarree())
fig_ax2.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax2.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax2.add_feature(cfeature.COASTLINE,lw=1)
fig_ax2.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax2.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax2.xaxis.set_major_formatter(lon_formatter)
fig_ax2.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax2.set_title('(b) EOF2(1979-2021)',loc='left',fontsize =15)
fig_ax2.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c2=fig_ax2.contourf(lon,lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
# 绘制第三模态
fig_ax3=fig.add_axes(,projection=proj)
fig_ax3.set_extent(,crs=ccrs.PlateCarree())
fig_ax3.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax3.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax3.add_feature(cfeature.COASTLINE,lw=1)
fig_ax3.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax3.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax3.xaxis.set_major_formatter(lon_formatter)
fig_ax3.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax3.set_title('(c) EOF3(1979-2021)',loc='left',fontsize =15)
fig_ax3.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c3=fig_ax3.contourf(lon,lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
cbposition=fig.add_axes()
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f')
#绘制时间序列
fig=plt.figure(figsize=(15,15))
fig_ax5=fig.add_axes()
fig_ax5.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax5.set_title('sst',loc='center',fontsize = 15)
fig_ax5.set_ylim(-3.5,3.5)
fig_ax5.axhline(0,linestyle="--")
fig_ax5.plot(np.arange(1979,2022,4/12),pc[:,0],color='blue')
fig_ax6 = fig.add_axes()
fig_ax6.set_title('(d) PC2',loc='left',fontsize = 15)
fig_ax6.set_ylim(-3.5,3.5)
fig_ax6.axhline(0,linestyle="--")
fig_ax6.plot(np.arange(1979,2022,4/12),pc[:,1],color='blue')
fig_ax7 = fig.add_axes()
fig_ax7.set_title('(f) PC3',loc='left',fontsize = 15)
fig_ax7.set_ylim(-3.5,3.5)
fig_ax7.axhline(0,linestyle="--")
fig_ax7.plot(np.arange(1979,2022,4/12),pc[:,2],color='blue')
##进行north检验
lat=np.array(lat)
coslat=np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver=Eof(sst,weights=wgts)
eof=solver.eofsAsCorrelation(neofs=10)
#此处的neofs的值是我们需要的空间模态数，比如这里我们打算画四个模态
pc=solver.pcs(npcs=10,pcscaling=1)#方差
var=solver.varianceFraction(neigs=10)
##计算特征值和NorthTest结果
eigenvalues = solver.eigenvalues(neigs=10)
errors = solver.northTest(neigs=10)
##绘制NorthTest结果Error Bar
plt.figure()
plt.errorbar(np.arange(1,11),eigenvalues,yerr=errors, fmt='o-', capsize=5 ,color = 'k')
plt.xlabel('EOF modes')
plt.ylabel('Eigenvalues')
plt.title('NorthTest(with error bar)_sst')
plt.show()
微信搜索“第八星系人造大气理论爱好者”公众号，关注获取文章数据
页:
[1]