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

更新時間：2022年09月26日 11:46:22 作者：oceanography-Rookie

這篇文章主要介紹了python對站點數(shù)據(jù)做EOF且做插值繪制填色圖，文章圍繞主題展開詳細的內(nèi)容介紹，具有一定的參考價值，，需要的小伙伴可以參考一下

前言

讀取站點資料數(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繪圖了，也當做是復習一下繪圖過程。

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

導入必要的庫：

from eofs.standard import Eof
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from metpy.interpolate import inverse_distance_to_grid

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

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

##################### read station  data   ##########################################
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年

插值

#####################   interp data   ##########################################
### interp data to target grid
### set target grid
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]))
for i in range(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的過程注釋掉了，需要測試的同學之間取消注釋即可。
注意點：插值過程需要先設(shè)置目標的插值網(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ù)并繪圖

##################### def  plot function ##########################################
def contourf(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='Times New Roman',size=10)
    plt.xticks(fontproperties='Times New Roman',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)

    return c

def p_line(ax,i):

    ax.set_title('pc'+str(i)+'',loc='left',fontsize = 15)
    # ax.set_ylim(-3.5,3.5)
    ax.axhline(0,linestyle="--")
    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()

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

展示結(jié)果

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

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

附上全部的繪圖代碼：

# -*- coding: utf-8 -*-
"""
Created on Fri Sep 23 17:46:42 2022

@author: Administrator
"""
from eofs.standard import Eof
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from metpy.interpolate import inverse_distance_to_grid

##################### read station  data   ##########################################
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]


#####################   interp data   ##########################################
### interp data to target grid
### set target grid
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]))
for i in range(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')

#計算緯度權(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)
    
##################### def  plot function ##########################################
def contourf(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='Times New Roman',size=10)
    plt.xticks(fontproperties='Times New Roman',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)
    
    return c

def p_line(ax,i):
    
    ax.set_title('pc'+str(i)+'',loc='left',fontsize = 15)
    # ax.set_ylim(-3.5,3.5)
    ax.axhline(0,linestyle="--")
    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()