第python對站點數據做EOF且做插值繪制填色圖目錄前言讀取存儲的數據插值EOF處理定義繪圖函數并繪圖展示結果總結

前言

讀取站點資料數據對站點數據進行插值，插值到規則網格上繪制EOF第一模態和第二模態的空間分布圖繪制PC序列

關于插值，這里主要提供了兩個插值函數，一個是一般常用的規則網格插值:

griddata

另一個是metpy中的:

inverse_distance_to_grid()

本來只是測驗一下不同插值方法，但是發現兩種插值方法的結果差別很大，由于對于站點數據處理較少，所以不太清楚具體原因。如果有知道的朋友可以告知一下，不甚感謝！

本次數據存儲的文件格式為.txt，讀取的方法是通過pandas.read_csv()

同時，之前一直嘗試使用proplot進行繪圖，較長時間不用matplotlib.pyplot繪圖了，也當做是復習一下繪圖過程。

繪圖中的代碼主要通過封裝函數，這樣做的好處是大大減少了代碼量。

導入必要的庫：

fromeofs.standardimportEof

importmatplotlib.pyplotasplt

importnumpyasnp

fromerpolateimportgriddata

importpandasaspd

importmatplotlib.pyplotasplt

importcartopy.crsasccrs

importcartopy.featureascfeature

fromcartopy.mpl.tickerimportLongitudeFormatter,LatitudeFormatter

fromerpolateimportinverse_distance_to_grid

出現找不到庫的報錯，這里使用condainstallpackagename安裝一下就好

讀取存儲的數據

#####################readstationdata##########################################

path=r'D:/data.txt'

file=pd.read_csv(path,sep="\t",

header=None,

names=['station','lat','lon','year','data'],

na_values=-99.90)

data=file['data'].to_numpy()

lon=file['lon'].to_numpy()

lat=file['lat'].to_numpy()

year=file['year'].to_numpy()

station=file['station'].to_numpy()

year_max=np.max(year)

year_min=np.min(year)

year_range=np.arange(year_min,year_max+1,1)

data_all=data.reshape(70,53)

lon_all=lon.reshape(70,53)/100

lat_all=lat.reshape(70,53)/100

lon_real=lon_all[:,0]

lat_real=lat_all[:,0]

這里將讀取的數據全部轉為array格式，方便查看以及后續處理。本來存儲的文件中是沒有相關的經度、緯度、站點、時間的名稱的，這里我是自己加在上面方面數據讀取的。

本次處理的數據包含70個站點，53年

插值

#####################interpdata##########################################

###interpdatatotargetgrid

###settargetgrid

lon_target=np.arange(115,135.5,0.5)

lat_target=np.arange(38,55.5,0.5)

x_t,y_t=np.meshgrid(lon_target,lat_target)

z=np.zeros((len(year_range),lat_target.shape[0],lon_target.shape[0]))

foriinrange(len(year_range)):

print(i)

#z[i]=inverse_distance_to_grid(lon_real,lat_real,

#data_all[:,i],

#x_t,y_t,r=15,min_neighbors=3)

z[i]=griddata((lon_real,lat_real),

data_all[:,i],

(x_t,y_t),method='cubic')

這里顯示了使用griddata()的插值過程，metpy的過程注釋掉了，需要測試的同學之間取消注釋即可。

注意點：插值過程需要先設置目標的插值網格。

EOF處理

#計算緯度權重

lat_new=np.array(lat_target)

coslat=np.cos(np.deg2rad(lat_new))

wgts=np.sqrt(coslat)[...,np.newaxis]

#創建EOF分解器

solver=Eof(z,weights=wgts)

eof=solver.eofsAsCorrelation(neofs=2)

#此處的neofs的值是我們需要的空間模態數

pc=solver.pcs(npcs=2,pcscaling=1)#方差

var=solver.varianceFraction(neigs=2)

這里沒啥好說的，需要幾個模態自由選擇即可

定義繪圖函數并繪圖

#####################defplotfunction##########################################

defcontourf(ax,i,level,cmap):

extents=[115,135,35,55]

ax.set_extent(extents,crs=proj)

ax.add_feature(cfeature.LAND,edgecolor='black',facecolor='silver',

ax.add_feature(cfeature.LAKES,edgecolor='black',facecolor='w',

ax.add_feature(cfeature.BORDERS)

xtick=np.arange(extents[0],extents[1],5)

ytick=np.arange(extents[2],extents[3],5)

ax.coastlines()

tick_proj=ccrs.PlateCarree()

c=ax.contourf(lon_target,lat_target,eof[i],

transform=ccrs.PlateCarree(),

levels=level,

extend='both',

cmap=cmap)

ax.set_xticks(xtick,crs=tick_proj)

ax.set_xticks(xtick,crs=tick_proj)

ax.set_yticks(ytick,crs=tick_proj)

ax.set_yticks(ytick,crs=tick_proj)

ax.xaxis.set_major_formatter(LongitudeFormatter())

ax.yaxis.set_major_formatter(LatitudeFormatter())

plt.yticks(fontproperties='TimesNewRoman',size=10)

plt.xticks(fontproperties='TimesNewRoman',size=10)

ax.tick_params(which='major',direction='out',

length=4,width=0.5,

pad=5,bottom=True,left=True,right=True,top=True)

ax.tick_params(which='minor',direction='out',

bottom=True,left=True,right=True,top=True)

ax.set_title('EOF'+str(i),loc='left',fontsize=15)

returnc

defp_line(ax,i):

ax.set_title('pc'+str(i)+'',loc='left',fontsize=15)

#ax.set_ylim(-3.5,3.5)

ax.axhline(0,line)

ax.plot(year_range,pc[:,i],color='blue')

ax.set_ylim(-3,3)

#####################plot##########################################

fig=plt.figure(figsize=(8,6),dpi=200)

proj=ccrs.PlateCarree()

contourf_1=fig.add_axes([0.02,0.63,0.5,0.3],projection=proj)

c1=contourf(contourf_1,0,np.arange(0.7,1,0.05),plt.cm.bwr)

plot_1=fig.add_axes([0.45,0.63,0.5,0.3])

p_line(plot_1,0)

contourf_2=fig.add_axes([0.02,0.15,0.5,0.3],projection=proj)

c2=contourf(contourf_2,1,np.arange(-0.5,0.6,0.1),plt.cm.bwr)

plot_2=fig.add_axes([0.45,0.15,0.5,0.3],)

p_line(plot_2,1)

cbposition1=fig.add_axes([0.16,0.55,0.22,0.02])

cb=fig.colorbar(c1,cax=cbposition1,

orientation='horizontal',format='%.1f')

cb.ax.tick_params(which='both',direction='in')

cbposition2=fig.add_axes([0.16,0.08,0.22,0.02])

cb2=fig.colorbar(c2,cax=cbposition2,

orientation='horizontal',format='%.1f')

cb2.ax.tick_params(which='both',direction='in')

plt.show()

這里將大部分重復的繪圖代碼，進行了封裝，通過封裝好的函數進行調用，大大簡潔了代碼量。相關的封裝過程之前也有講過，可以翻看之前的記錄。

展示結果

使用griddata的繪圖結果：

使用metpt插值函數的結果：

附上全部的繪圖代碼：

#-*-coding:utf-8-*-

CreatedonFriSep2317:46:422025

@author:Administrator

fromeofs.standardimportEof

importmatplotlib.pyplotasplt

importnumpyasnp

fromerpolateimportgriddata

importpandasaspd

importmatplotlib.pyplotasplt

importcartopy.crsasccrs

importcartopy.featureascfeature

fromcartopy.mpl.tickerimportLongitudeFormatter,LatitudeFormatter

fromerpolateimportinverse_distance_to_grid

#####