python對站點數據做EOF且做插值繪制填色圖
更新時間：2022年09月26日 11:46:22 作者：oceanography-Rookie
這篇文章主要介紹了python對站點數據做EOF且做插值繪制填色圖，文章圍繞主題展開詳細的內容介紹，具有一定的參考價值，，需要的小伙伴可以參考一下
前言
讀取站點資料數據對站點數據進行插值，插值到規(guī)則網格上繪制EOF第一模態(tài)和第二模態(tài)的空間分布圖繪制PC序列
關于插值，這里主要提供了兩個插值函數，一個是一般常用的規(guī)則網格插值:
griddata
另一個是metpy中的:
inverse_distance_to_grid()
本來只是測驗一下不同插值方法，但是發(fā)現兩種插值方法的結果差別很大，由于對于站點數據處理較少，所以不太清楚具體原因。如果有知道的朋友可以告知一下，不甚感謝！
本次數據存儲的文件格式為.txt，讀取的方法是通過pandas.read_csv()
同時，之前一直嘗試使用proplot進行繪圖，較長時間不用matplotlib.pyplot繪圖了，也當做是復習一下繪圖過程。
繪圖中的代碼主要通過封裝函數，這樣做的好處是大大減少了代碼量。
導入必要的庫：
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
出現找不到庫的報錯，這里使用conda install packagename 安裝一下就好
讀取存儲的數據
##################### 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]這里將讀取的數據全部轉為array格式，方便查看以及后續(xù)處理。本來存儲的文件中是沒有相關的經度、緯度、站點、時間的名稱的，這里我是自己加在上面方面數據讀取的。
本次處理的數據包含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的過程注釋掉了，需要測試的同學之間取消注釋即可。
注意點：插值過程需要先設置目標的插值網格。
EOF處理
#計算緯度權重 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)數 pc=solver.pcs(npcs=2,pcscaling=1)#方差 var=solver.varianceFraction(neigs=2)
這里沒啥好說的，需要幾個模態(tài)自由選擇即可
定義繪圖函數并繪圖
##################### 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()這里將大部分重復的繪圖代碼，進行了封裝，通過封裝好的函數進行調用，大大簡潔了代碼量。相關的封裝過程之前也有講過，可以翻看之前的記錄。
展示結果
使用griddata的繪圖結果：
使用metpt插值函數的結果：
附上全部的繪圖代碼：
# -*- 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')
#計算緯度權重
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)數
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()總結
metpy的插值函數好處在于可以自由填充整個繪圖區(qū)域，但是感覺griddata函數的插值結果更加符合預期，雖然也有點怪怪的。
到此這篇關于python對站點數據做EOF且做插值繪制填色圖的文章就介紹到這了,更多相關python EOF插值繪制內容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關文章希望大家以后多多支持腳本之家！
相關文章
Python處理字節(jié)串:struct.pack和struct.unpack使用
這篇文章主要介紹了Python處理字節(jié)串:struct.pack和struct.unpack使用方式,具有很好的參考價值,希望對大家有所幫助,如有錯誤或未考慮完全的地方,望不吝賜教2024-01-01
生活中除了給男神女神表白,還會遇到朋友過生日,所以小編今天就給大家?guī)砹艘粋€送給朋友的生日祝福程序,這篇文章主要給大家介紹了關于如何用python給朋友做生日快樂祝福的相關資料,需要的朋友可以參考下2023-11-11
這篇文章主要介紹了django 創(chuàng)建過濾器的實例詳解的相關資料,主要說明django 創(chuàng)建過濾器來統(tǒng)一處理字符串，需要的朋友可以參考下2017-08-08
詳解利用django中間件django.middleware.csrf.CsrfViewMiddleware防止csrf
這篇文章主要介紹了詳解利用django中間件django.middleware.csrf.CsrfViewMiddleware防止csrf攻擊，小編覺得挺不錯的，現在分享給大家，也給大家做個參考。一起跟隨小編過來看看吧2018-10-10
python?泛型函數--singledispatch的使用解讀
這篇文章主要介紹了python?泛型函數--singledispatch的使用解讀，具有很好的參考價值，希望對大家有所幫助。如有錯誤或未考慮完全的地方，望不吝賜教2022-09-09
在VS Code上搭建Python開發(fā)環(huán)境的方法
這篇文章主要介紹了在VS Code上搭建Python開發(fā)環(huán)境的方法,需要的朋友可以參考下2018-04-04
這篇文章主要為大家介紹了Python語法中實現格式化輸出的方法，本文通過幾個實例為大家進行了詳細的講解，感興趣的小伙伴可以了解一下2022-08-08
這篇文章主要介紹了利用一個簡單的例子窺探CPython內核的運行機制,作者通過一個簡單的輸出函數深入、介紹了CPython源碼C代碼中的一些函數，需要的朋友可以參考下2015-03-03
這篇文章主要介紹了基于PyInstaller各參數的含義說明，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧2021-03-03
本篇內容給大家詳細講述了Python3的TFTP文件傳輸的相關知識點，有需要的朋友可以參考下。2018-06-06