python對站點數(shù)據(jù)做EOF且做插值繪制填色圖

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

前言

讀取站點資料數(shù)據(jù)對站點數(shù)據(jù)進行插值，插值到規(guī)則網(wǎng)格上繪制EOF第一模態(tài)和第二模態(tài)的空間分布圖繪制PC序列

關(guān)于插值，這里主要提供了兩個插值函數(shù)，一個是一般常用的規(guī)則網(wǎng)格插值:

griddata

另一個是metpy中的:

inverse_distance_to_grid()

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

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

同時，之前一直嘗試使用proplot進行繪圖，較長時間不用matplotlib.pyplot繪圖了，也當(dāng)做是復(fù)習(xí)一下繪圖過程。

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

導(dǎo)入必要的庫：

fromeofs.standardimportEof

importmatplotlib.pyplotasplt

importnumpyasnp

fromerpolateimportgriddata

importpandasaspd

importmatplotlib.pyplotasplt

importcartopy.crsasccrs

importcartopy.featureascfeature

fromcartopy.mpl.tickerimportLongitudeFormatter,LatitudeFormatter

fromerpolateimportinverse_distance_to_grid

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

讀取存儲的數(shù)據(jù)

#####################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]

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

本次處理的數(shù)據(jù)包含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的過程注釋掉了，需要測試的同學(xué)之間取消注釋即可。

注意點：插值過程需要先設(shè)置目標(biāo)的插值網(wǎng)格。

EOF處理

#計算緯度權(quán)重

lat_new=np.array(lat_target)

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

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

#創(chuàng)建EOF分解器

solver=Eof(z,weights=wgts)

eof=solver.eofsAsCorrelation(neofs=2)

#此處的neofs的值是我們需要的空間模態(tài)數(shù)

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

var=solver.varianceFraction(neigs=2)

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

定義繪圖函數(shù)并繪圖

#####################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()

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

展示結(jié)果

使用griddata的繪圖結(jié)果：

使用metpt插值函數(shù)的結(jié)果：

附上全部的繪圖代碼：

#-*-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

#####