一、反距離權重插值
假設：彼此距離較近的事物要比彼此距離較遠的事物更相似。當為任何未測量的位置預測值時，反距離權重法會采用預測位置周圍的測量值與距離預測位置較遠的測量值相比，距離預測位置最近的測量值對預測值的影響更大。反距離權重法假定每個測量點都有一種局部影響，而這種影響會隨著距離的增大而減小。由于這種方法為距離預測位置最近的點分配的權重較大，而權重卻作為距離的函數(shù)而減小，因此稱之為反距離權重法。
1.1 庫導入
import pandas as pd
import numpy as np
import geopandas as gpd
import netCDF4 as nc
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import matplotlib.pyplot as plt
from matplotlib.image import imread
from cartopy.io.shapereader import Reader
1.2 站點觀測數(shù)據(jù)
# 站點觀測
data = pd.read_excel('weather_station_data.xlsx', 'Sheet1', na_values=['NA'])
lon_sta = data['經度'][:].to_numpy()
lat_sta = data['緯度'][:].to_numpy()
t2m_sta = data['t2m_bilinear'][:].to_numpy()
1.3 半正矢公式（Haversine）計算球面兩點的距離
#給定經緯度，計算地球上任意兩點距離，單位m
import numpy as np
def haversine_dist(lon1,lat1,lon2,lat2):
lon1,lat1,lon2,lat2 = map(np.radians, (lon1,lat1,lon2,lat2))
radius = 6378.135E3 # radius of Earth, unit:m
dlat = lat2 - lat1
dlon = lon2 - lon1
arg = np.sin(dlat/2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2) ** 2
dist = 2 * radius * np.arcsin(np.sqrt(arg))
return dist
1.4 IDW插值
def interp_IDW(lon_sta, lat_sta, data_sta, lon2D, lat2D):
n_sta = len(lon_sta)
ny, nx = np.shape(lon2D)
data2D = np.zeros((ny, nx))
for j in range(ny):
for i in range(nx): #遍歷二維每一個格點
dist = [] # 格點至所有站點的距離
for s in range(n_sta):
d = haversine_dist(lon_sta[s], lat_sta[s], lon2D[j,i], lat2D[j,i])
d = np.max([1.0, d]) # aviod divide by zero
dist.append(d)
wgt = 1.0 / np.power(dist, 2)
wgt_sum = np.sum(wgt)
arg_sum = np.sum(np.array(wgt) * np.array(data_sta))
data2D[j,i] = arg_sum / wgt_sum
return data2D
1.5 插值結果
# 插值的結果
t2m_IDW = interp_IDW(lon_sta, lat_sta, t2m_sta, lon2D, lat2D)
1.6 繪制反距離權重插值結果
# 繪制IDW插值結果
plot_contour_scatter(t2m_IDW, lon2D, lat2D, lon_sta, lat_sta, 't2m_IDW.png')
二、克里金法（Kriging）
Kriging是依據(jù)協(xié)方差函數(shù)對隨機過程/隨機場進行空間建模和預測（插值）的回歸算法。在特定的隨機過程，例如固有平穩(wěn)過程中，克里金法能夠給出最優(yōu)線性無偏估計（Best Linear Unbiased Prediction, BLUP），因此在地統(tǒng)計學中也被稱為空間最優(yōu)無偏估計器(spatial BLUP)?？死锝穑╧riging）插值是在有限區(qū)域內對區(qū)域化變量進行無偏最優(yōu)估計的一種方法（用于估計在空間上有相關性的值，比如空氣質量，相隔很近的位置的數(shù)值接近）。無偏指的是估計值和實際值之差的期望等于零，最優(yōu)指的是估計值和實際值的方差最小。基于這一特點使得克里金插值的效果比其他插值方法要好很多。
普通克里金（Ordinary Kriging, OK） 泛克里金（Universal Kriging, UK） 協(xié)同克里金（Co-Kriging, CK） 析取克里金（Disjunctive Kriging, DK） 回歸克里金（regression-Kriging） 神經網絡克里金（neural Kriging） 貝葉斯克里金（Bayesian Kriging） 在Python里，有兩個常用的克里金插值包，pykrige和pykriging。
總體而言，pykrige相對方便好用，支持普通克里金、泛克里金和回歸克里金。
github網址：https://github.com/whdc/PyKrige
2.1 普通克里金用法，其他同理
1. 構造 ordinary kriging object
OK_obj = OrdinaryKriging(lon_sta, lat_sta, t2m_sta, variogram_model="linear")
輸入：
lon_sta: 站點經度，一維
lat_sta: 站點緯度，一維
t2m_sta：站點數(shù)據(jù)，一維，這里是溫度
參數(shù)：
variogram_model：是變差函數(shù)模型，pykrige提供 linear, power, gaussian, spherical,
exponential, hole-effect幾種選擇，默認的為linear模型。
輸出：
OK_obj : 一個對象，pykrige.ok.OrdinaryKriging
2. 輸出插值網格點的值
t2m_OK, ss = OK_obj.execute("grid", longitude, latitude) # 網格點處對應的值
輸入：網格的經度和緯度，一維
輸出：網格值和方差
2.2 普通克里金插值結果
from pykrige.ok import OrdinaryKriging
### 普通克里金（Ordinary Kriging, OK）
OK_obj = OrdinaryKriging(lon_sta, lat_sta, t2m_sta, variogram_model="linear")
# variogram_model是變差函數(shù)模型，
# pykrige提供 linear, power, gaussian, spherical, exponential, hole-effect幾種選擇，默認的為linear模型。
# 使用不同的variogram_model，插值效果是不一樣的，應該針對自己的任務多嘗試不同選項。
t2m_OK, ss = OK_obj.execute("grid", longitude, latitude) # 網格點處對應的值
2.3 泛克里金（Universal Kriging, UK）插值結果
# 泛克里金
from pykrige.uk import UniversalKriging
### 泛克里金（Universal Kriging, UK）
UK_obj = UniversalKriging(lon_sta, lat_sta, t2m_sta, variogram_model="linear")
t2m_UK, ss = UK_obj.execute("grid", longitude, latitude) # 給出網格點處對應的值
2.4 繪制普通克里金法插值結果
# 繪制克里金插值結果
plot_contour_scatter(t2m_OK, lon2D, lat2D, lon_sta, lat_sta, 't2m_OK.png')
三、徑向基函數(shù)（RBF）插值
3.1 RBF方法是一系列精確插值方法的組合；即表面必須通過每一個測得的采樣值。
有以下五種基函數(shù)：
- 薄板樣條函數(shù)
- 張力樣條函數(shù)
- 規(guī)則樣條函數(shù)
- 高次曲面函數(shù)
- 反高次曲面函數(shù)
3.2 RBF 插值結果
# cubic, gaussian, inverse_multiquadric, linear, multiquadric, quintic, thin_plate
from scipy import interpolate
func = interpolate.Rbf(lon_sta, lat_sta, t2m_sta, function='multiquadric')
t2m_RBF = func(lon2D,lat2D)#輸入輸出都是二維
3.3 繪制徑向基函數(shù)（RBF）插值結果
# 繪制徑向基函數(shù)插值結果
plot_contour_scatter(t2m_RBF, lon2D, lat2D, lon_sta, lat_sta, 't2m_RBF.png')
四、結果可視化和對比
4.1 導庫
import xarray as xr
import salem
import numpy as np
import geopandas as gpd
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.image import imread
from cartopy.io.shapereader import Reader
4.2 驗證數(shù)據(jù)
## ERA5 data，用于驗證插值效果，并提供待插值的網格 lon2D，lat2D
## 北京地區(qū)2020年1，4，7，10月份（月平均）的地面數(shù)據(jù)
dataset = nc.Dataset("../data/ERA5_single_level_2020_1-4-7-10_beijing.nc")
print(dataset.variables.keys())
# 經緯度
longitude = dataset.variables['longitude'][:].data
latitude = dataset.variables['latitude'][:].data
# 第2時次為7月份的數(shù)據(jù)
t2m_ERA5 = dataset.variables['t2m'][2].data
lon2D, lat2D = np.meshgrid(longitude,latitude)
4.3 定義通用方法
# 統(tǒng)一定制繪制方法
def plot_contour_scatter_subplot(data_list, data_name, lons2D, lats2D, lonSta, latSta, figname):
lonMin = np.min(lons2D)
lonMax = np.max(lons2D)
latMin = np.min(lats2D)
latMax = np.max(lats2D)
proj = ccrs.PlateCarree() # 創(chuàng)建坐標系
fig = plt.figure(figsize=(16, 18), dpi=400) # 創(chuàng)建頁面
axes = fig.subplots(2,2, subplot_kw={'projection': proj}) # 創(chuàng)建子圖
rows = [0, 0, 1, 1]
cols = [0, 1, 0, 1]
for i in range(len(data_list)):
ax = axes[rows[i], cols[i]]
data2D = data_list[i]
ax.set_extent([lonMin, lonMax, latMin, latMax])
# --設置地圖屬性
provinces = cfeat.ShapelyFeature(
gpd.read_file('china_shp/province.shp')['geometry'],
proj,
edgecolor='k',
facecolor='none'
)
line = cfeat.ShapelyFeature(
Reader('china_shp/nine_line.shp').geometries(),
ccrs.PlateCarree(),
edgecolor='k',
facecolor='none'
)
ax.add_feature(provinces, linewidth=0.6, zorder=1)
ax.add_feature(line, linewidth=0.6, zorder=1)
ax.add_feature(cfeat.RIVERS.with_scale('110m'), linewidth=0.5, zorder=1)
ax.add_feature(cfeat.LAKES.with_scale('110m'), linewidth=0.5, zorder=1)
gl = ax.gridlines(
crs = ccrs.PlateCarree(),
draw_labels = False,
linewidth = 0.9,
color = 'k',
alpha = 0.5,
linestyle = '--'
)
gl.top_labels = False # 關閉經緯度標簽
gl.right_labels =False
# --設置刻度
ax.set_xticks(np.linspace(int(lonMin), int(lonMax), num=5, endpoint=True))
ax.set_yticks(np.linspace(int(latMin), int(latMax), num=5, endpoint=True))
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.xaxis.set_minor_locator(plt.MultipleLocator(1))
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.yaxis.set_minor_locator(plt.MultipleLocator(1))
ax.tick_params(axis='both', labelsize=5, direction='out')
levels = np.linspace(290.0, 300.0, 11, endpoint=True)
ctr = ax.contourf(lons2D, lats2D, data2D,
levels=levels,
cmap=get_cmap("rainbow"),
extend='both')
cb = plt.colorbar(ctr, ax=ax, orientation="vertical", pad=.02, fraction=0.05)
cb.ax.tick_params(labelsize=5)
ax.set_title(data_name[i], {"fontsize" : 15})
ax.scatter(
lonSta,
latSta,
marker='*',
s=8 ,
color ="black"
)
plt.savefig(figname)
return None
4.4 三種插值結果和驗證數(shù)據(jù)對比
文章來源：http://www.zghlxwxcb.cn/news/detail-447420.html
4.5 三種插值結果對比差異
diff_IDW = t2m_IDW - t2m_ERA5
diff_OK = t2m_OK - t2m_ERA5
diff_RBF = t2m_RBF - t2m_ERA5
data_list = [diff_IDW, diff_OK, diff_RBF ]
data_name = ['IDW-ERA5', 'OK-ERA5', 'RBF-ERA5']
plot_contour_scatter_diff(data_list, data_name, lon2D, lat2D, lon_sta, lat_sta, 'diff_interp.png')
文章來源地址http://www.zghlxwxcb.cn/news/detail-447420.html
到了這里，關于Python 站點數(shù)據(jù)插值到格點 反距離權重插值 克里金法 徑向基函數(shù)（RBF）插值的文章就介紹完了。如果您還想了解更多內容，請在右上角搜索TOY模板網以前的文章或繼續(xù)瀏覽下面的相關文章，希望大家以后多多支持TOY模板網！
