EOF分析

沉默的狙击手

数据资料介绍

NCEP 1948-2020年1~12月的17层月平均高度场资料，资料范围为90S-90N，0-360，分辨率2.5*2.5，纬向格点数144，经向格点数73.

代码：

import numpy as np

import xarray as xr

import matplotlib.pyplot as plt

from matplotlib.pyplot import MultipleLocator

import cartopy.crs as ccrs

import cartopy.feature as cft

from cartopy.util import add_cyclic_point

from eofs.standard import Eof

plt.rc('font',family="Times New Roman")

font1={'family':'Times New Roman','weight':'bold','size':10}

data = xr.open_dataset('E:\\duanqi\\task1\\hgt.mon.mean.nc')

#print(data.variables.keys())

#print(data.variables.items())

# for i in data.variables.keys():

#     print(i)

z=data.hgt.squeeze()/10

z=z[:,5,:,:].squeeze()#500hPa位势高度

# print(z.shape)

lat = data['lat'].loc[70:20]

lon = data['lon'].loc[40:140]

t1=[]

tt='-01-01T00:00:00.000000000'

for i in range(1979,2020):

    m=str(i)+tt

    t1.append(m)

z=z.loc['1979-01-01T00:00:00.000000000':'2019-12-01T00:00:00.000000000'].loc[:,70:20,40:140]

z1=[]

for i in t1:

    z1.append(z.loc)

z=np.array(z)

z1=np.array(z1)

date = []

for i in range(41):

    date.append(i)

eof = Eof(z1)

PCs = eof.pcs(npcs = 3)

EOFs = eof.eofs(neofs = 3)*100

EOF_var = eof.varianceFraction(neigs = 3)

fig=plt.figure(figsize=(8,8))

ax1=plt.subplot(321,projection=ccrs.PlateCarree(central_longitude=90))

ax2=plt.subplot(322)

ax3=plt.subplot(323,projection=ccrs.PlateCarree(central_longitude=90))

ax4=plt.subplot(324)

ax5=plt.subplot(325,projection=ccrs.PlateCarree(central_longitude=90))

ax6=plt.subplot(326)

plt.subplots_adjust(left=0.06, right=0.94,

                    bottom=0.1, top=0.95,

                      wspace=0.2, hspace=0.3)

##########################

ax1 = plt.subplot(321,projection =

                    ccrs.PlateCarree(central_longitude= 0))

ax1.coastlines()

ax1.add_feature(cft.OCEAN, color = 'lightblue')

ax1.add_feature(cft.LAKES, facecolor = 'none',

                edgecolor = 'gray', linewidth = 0.5)

### 添加网格和坐标 =========================================

ax1.gridlines()

gl1 = ax1.gridlines(draw_labels = True, linestyle = '--',

                   linewidth = 0.5, color ='gray',

                   x_inline=False, y_inline=False)

gl1.top_labels = False

gl1.right_labels = False

cycle_eof1,cycle_lon=add_cyclic_point(EOFs[0][:][:],coord=lon)

cycle_LON,cycle_LAT=np.meshgrid(cycle_lon,lat)

cf1 = ax1.contourf(cycle_LON,cycle_LAT,cycle_eof1,levels=np.arange(-2,9,1),cmap='YlOrBr',transform=ccrs.PlateCarree(),

                  vmax=9,vmin=-2,extend='both')

cb1 = fig.colorbar(cf1, orientation='horizontal',ticks = np.arange(-2, 9, 1),

                  fraction=0.08, pad = 0.12)

#设置坐标参数

cb1.set_label('Geopotential Height (gpm)',fontdict=font1)

cb1.ax.tick_params(labelsize=10)

plt.tick_params(labelsize=12)

ax1.set_title('(a) EOF1',font1,pad=7,loc = 'left')

ax1.set_title('%.2f%%' % (EOF_var[0]*100),font1,pad=7,loc = 'right')

############################

ax2.plot(date,PCs[:,0],c='k',linewidth = 0.5)

# ax2.set_xticks([1979,1984,1989,1994,1999,2004,2009,2014,2019])

x_major_locator = MultipleLocator(5)

ax2.xaxis.set_major_locator(x_major_locator)

#ax2.xlim(-0.5,492)

ax2.set_xticklabels(['','1979','1984','1989','1994','1999','2004','2009','2014','2019'])

ax2.set_title('(b) PC1',font1,pad=9,loc = 'left')

#############################

'''

ax3 = plt.subplot(323,projection =

                    ccrs.PlateCarree(central_longitude= 0))

ax3.coastlines()

ax3.add_feature(cft.OCEAN, color = 'lightblue')

ax3.add_feature(cft.LAKES, facecolor = 'none',

                edgecolor = 'gray', linewidth = 0.5)

### 添加网格和坐标 =========================================

ax3.gridlines()

gl3 = ax3.gridlines(draw_labels = True, linestyle = '--',

                   linewidth = 0.5, color ='gray',

                   x_inline=False, y_inline=False)

gl3.top_labels = False

gl3.right_labels = False

cycle_eof2,cycle_lon=add_cyclic_point(EOFs[1][:][:],coord=lon)

cycle_LON,cycle_LAT=np.meshgrid(cycle_lon,lat)

cf3 = ax3.contourf(cycle_LON,cycle_LAT,cycle_eof2,levels=np.arange(-7,7,1),cmap='RdBu_r',transform=ccrs.PlateCarree(),

                  vmax=7,vmin=-7,extend='both')

cb3 = fig.colorbar(cf3, orientation='horizontal',ticks = np.arange(-7, 7, 1),

                  fraction=0.08, pad = 0.12)

#设置坐标参数

cb3.set_label('Geopotential Height (gpm)',fontdict=font1)

cb3.ax.tick_params(labelsize=10)

plt.tick_params(labelsize=12)

ax3.set_title('(c) EOF2',font1,pad=7,loc = 'left')

ax3.set_title('%.2f%%' % (EOF_var[1]*100),font1,pad=7,loc = 'right')

#############################

ax4.plot(date,PCs[:,1],c='k',linewidth = 0.5)

# ax2.set_xticks([1979,1984,1989,1994,1999,2004,2009,2014,2019])

x_major_locator = MultipleLocator(5)

ax4.xaxis.set_major_locator(x_major_locator)

#ax2.xlim(-0.5,492)

ax4.set_xticklabels(['','1979','1984','1989','1994','1999','2004','2009','2014','2019'])

ax4.set_title('(d) PC2',font1,pad=9,loc = 'left')

##############################

ax5 = plt.subplot(325,projection = ccrs.PlateCarree(central_longitude= 0))

ax5.coastlines()

ax5.add_feature(cft.OCEAN, color = 'lightblue')

ax5.add_feature(cft.LAKES, facecolor = 'none',

                edgecolor = 'gray', linewidth = 0.5)

### 添加网格和坐标 =========================================

ax5.gridlines()

gl5 = ax5.gridlines(draw_labels = True, linestyle = '--',

                   linewidth = 0.5, color ='gray',

                   x_inline=False, y_inline=False)

gl5.top_labels = False

gl5.right_labels = False

cycle_eof3,cycle_lon=add_cyclic_point(EOFs[2][:][:],coord=lon)

cycle_LON,cycle_LAT=np.meshgrid(cycle_lon,lat)

cf5 = ax5.contourf(cycle_LON,cycle_LAT,cycle_eof3,levels=np.arange(-6,6,1),cmap='RdBu_r',transform=ccrs.PlateCarree(),

                  vmax=6,vmin=-6,extend='both')

cb5 = fig.colorbar(cf5, orientation='horizontal',ticks = np.arange(-6, 6, 1),

                  fraction=0.08, pad = 0.12)

#设置坐标参数

cb5.set_label('Geopotential Height (gpm)',fontdict=font1)

cb5.ax.tick_params(labelsize=10)

plt.tick_params(labelsize=12)

ax5.set_title('(e) EOF3',font1,pad=7,loc = 'left')

ax5.set_title('%.2f%%' % (EOF_var[2]*100),font1,pad=7,loc = 'right')

#############################

ax6.plot(date,PCs[:,2],c='k',linewidth = 0.5)

# ax2.set_xticks([1979,1984,1989,1994,1999,2004,2009,2014,2019])

x_major_locator = MultipleLocator(5)

ax6.xaxis.set_major_locator(x_major_locator)

#ax2.xlim(-0.5,492)

ax6.set_xticklabels(['','1979','1984','1989','1994','1999','2004','2009','2014','2019'])

ax6.set_title('(f) PC3',font1,pad=9,loc = 'left')

plt.savefig('E:\\duanqi\\task2\\picture1.eps')

'''

乌衣年少

请问为什么要乘100？

西红柿毛毛雨

乌衣年少 发表于 2023-10-17 10:17
请问为什么要乘100？

将解释方差转为百分数

fuyixin

楼主，用你的代码重新运行一遍，遇到错误，Z1是一维的？
File f:\data\untitled0.py:38
    eof = Eof(z1)

  File ~\anaconda3\Lib\site-packages\eofs\standard.py:112 in __init__
    raise ValueError('the input data set must be '

ValueError: the input data set must be at least two dimensional

fuyixin

for i in t1:
    z1.append(z.loc)

这个循环里，没有体现 i 呀？