怎么使用Python對站點數(shù)據(jù)執(zhí)行徑向基函數(shù)插值可視化

這篇文章給大家介紹怎么使用Python對站點數(shù)據(jù)執(zhí)行徑向基函數(shù)插值可視化，內(nèi)容非常詳細(xì)，感興趣的小伙伴們可以參考借鑒，希望對大家能有所幫助。

成都創(chuàng)新互聯(lián)服務(wù)項目包括張掖網(wǎng)站建設(shè)、張掖網(wǎng)站制作、張掖網(wǎng)頁制作以及張掖網(wǎng)絡(luò)營銷策劃等。多年來，我們專注于互聯(lián)網(wǎng)行業(yè)，利用自身積累的技術(shù)優(yōu)勢、行業(yè)經(jīng)驗、深度合作伙伴關(guān)系等，向廣大中小型企業(yè)、政府機(jī)構(gòu)等提供互聯(lián)網(wǎng)行業(yè)的解決方案，張掖網(wǎng)站推廣取得了明顯的社會效益與經(jīng)濟(jì)效益。目前，我們服務(wù)的客戶以成都為中心已經(jīng)輻射到張掖省份的部分城市，未來相信會繼續(xù)擴(kuò)大服務(wù)區(qū)域并繼續(xù)獲得客戶的支持與信任！

以云南省2015年6月的29個氣象站點數(shù)據(jù)為例進(jìn)行徑向基函數(shù)(Rbf)插值。數(shù)據(jù)格式如下：

今天需要使用到cartopy庫來繪圖，因此需要先安裝好，據(jù)說安裝很煩人，可以去uci下載.whl文件來安裝，安裝好后先測試一下是否可以運行，如下簡單測試：

怎么使用Python對站點數(shù)據(jù)執(zhí)行徑向基函數(shù)插值可視化

首先，這是一個不成功的嘗試，因為沒能成功加載shp圖層導(dǎo)致最后的插值沒有落在特定的地理范圍內(nèi)。如果有伙伴知道這個問題的解決方法，希望不吝賜教。我相信只要是問題就一定存在解決方案。
那么現(xiàn)在就開始今天的測試，首先還是老套路導(dǎo)入所需要的庫或模塊，如下：

import numpy as np
import pandas as pd
from scipy.interpolate import Rbf
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.io.shapereader import Reader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER,LATITUDE_FORMATTER
import matplotlib.pyplot as plt

怎么使用Python對站點數(shù)據(jù)執(zhí)行徑向基函數(shù)插值可視化

然后設(shè)置能夠正常顯示中文和負(fù)號

plt.rcParams['font.sans-serif']=['SimHei']#正常顯示中文
plt.rcParams['axes.unicode_minus']=False#正常顯示負(fù)號

怎么使用Python對站點數(shù)據(jù)執(zhí)行徑向基函數(shù)插值可視化

接下來用pandas讀取整理好的數(shù)據(jù)，使用的數(shù)據(jù)是云南省2015年6月的降水?dāng)?shù)據(jù)和氣溫數(shù)據(jù)，文中只測試了降水的Rbf插值。

climate = pd.read_excel('yunnanclimate.xlsx',header=None,sep='\s+',names=['站號','lon','lat','降水','氣溫'])#print(climate)

然后提取加載的xlsx數(shù)據(jù)里面的特定要素，比如經(jīng)緯度，降水啥的

#讀取climate里面的具體數(shù)據(jù)
lon = climate['lon']
lat = climate['lat']
rain_data = climate['降水']

然后設(shè)置范圍和格網(wǎng)化

#設(shè)置經(jīng)緯度范圍
slon = np.linspace(97.85,105.21,36)
slat = np.linspace(21.45,28.50,36)
slon,slat = np.meshgrid(slon,slat)

接下來開始插值，function，我這里測試了3種，linear+cubic+gaussian,對于他們有啥區(qū)別哈哈，在最后我將分別展示這三種函數(shù)的模樣

#開始插值，使用徑向基函數(shù)Rbf
func = Rbf(lon,lat,rain_data,function='cubic')
rainX = func(slon,slat)

接下來加載云南省的shp圖層以及修改經(jīng)緯度啥的信息，就是這個位置我出錯了，如果有知道解決方案的還請不吝賜教。

shp_path = r'D:\GeoCAS\Python\yunnan.shp'#設(shè)置絕對路徑
proj = ccrs.PlateCarree()#設(shè)置投影
fig = plt.figure(figsize=(5,5),dpi=300)
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})
extent = [97.85,105.21,21.45,28.50]#限定經(jīng)緯度范圍
reader = Reader(shp_path)#讀取shp矢量
yunnan = cfeat.ShapelyFeature(reader.geometries(),proj,edgecolor='k',facecolor='none')#設(shè)定shp特征
ax.add_feature(yunnan,linewidth=0.5)
ax.set_extent(extent,crs=proj)
gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,linestyle='-')
gl.xlabels_top=False
gl.ylabels_right=False
gl.xformatter = LONGITUDE_FORMATTER#以經(jīng)緯度格式顯示xy軸
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style={'size':6}#設(shè)置xy字體大小
gl.ylabel_style={'size':6}
h = ax.contourf(slon,slat,rainX,cmap='viridis',extend='both')
plt.colorbar(h,orientation='vertical')

plt.show()