学习Python中的逆向距离加权插值

地理空间插值是一个过程，用于使用已知的值来估计地理区域内未知点的值。

反距离加权，简称IDW，是用于地理空间数据插值的最流行的方法之一。这篇文章将教我们如何在Python中进行IDW插值。

IDW插值方法假定较近的数值比，较远的数值更有关系。IDW通过使用已知的值来估计未知点的值，其权重是它们与要估计的点的距离。

我们还可以在IDW中使用另一个叫做 "功率 "的变量来控制已知值的影响。较高的功率将增加附近数值的影响。

下面是一个在Python中对IDW的简单DIY实现。

Python

import numpy as np
def idw_custom(dist,val,power):
    numerator=np.sum(val/(dist**power))
    weights=np.sum(1/(dist**power))
    return numerator/weights
dist=np.array([1.38,3.04,7.36,7.59,9.67])
val=np.array([17.7,15.7,34.3,6.3,23.3])
print(idw_custom(dist,val,1))
分子包含已知值的总和除以距离，分母包括距离的倒数之和。功率可以用来控制已知值的影响。
使用内置库来进行IDW插值
在上面的例子中，我们使用了预先计算的距离，但在大多数实际使用情况下，我们必须自己计算距离。我们可以使用哈维线距离来计算，但是当我们有很多点的时候，这可能是很麻烦的。这时，我们可以使用预先存在的库来计算距离并进行插值。
在我们的例子中，我们将对班加罗尔市的PM2.5值进行插值。
Python
grid_space = 0.01
grid_lon = np.arrange(np.amin(lons), np.amax(lons), grid_space)
grid_lat = np.arrange(np.amin(lats), np.amax(lats), grid_space)
让我们首先生成一个我们需要估计数值的点的网格。我们已经设置了大约1公里的网格空间。
"lons "包含一个经度列表，"lats "包含一个纬度列表。我们使用经度和纬度的最小值和最大值来生成网格。
Python
all_lats = np.meshgrid(grid_lon, grid_lat)[1].ravel()
all_lons = np.meshgrid(grid_lon, grid_lat)[0].ravel()
itrp=pd.DataFrame()
itrp['lat']=all_lats
itrp['lng']=all_lons
在上面的代码中，我们创建了一个数据框架，其中包含我们需要估计数值的经度和纬度对。我们也可以使用 "for循环 "来做同样的事情，如下图所示。
Python
lat=[]
lng=[]
for i in range(len(grid_lat)):
    for j in range(len(grid_lon)):
        lat.append(grid_lat[i])
        lng.append(grid_lon[j])
itrp=pd.DataFrame()
itrp['lat']=lat
itrp['lng']=lng
我们可以使用Sklearn的KNN实现来模拟IDW。下面给出的代码就是这样做的。
Python
x=sample[['lat','lng']]
y=sample['PM2.5']
from sklearn.neighbors import KNeighborsRegressor
model=KNeighborsRegressor(algorithm='kd_tree',n_neighbors=8,weights='distance').fit(x,y)
样本 "数据帧包含单个时间戳的站点空气质量数据。我们给出经纬度作为解释变量，PM2.5作为需要插值的变量。我们应该使用 "kd_tree "作为算法，并设置 "n_neighbors "作为站点的数量，在本例中是8个。我们还应该设置 "weights "作为执行IDW的距离。
Python
pred=model.predict(itrp[['lat','lng']])
我们将使用predict方法来估计我们的网格点的数值，这些网格点存储在itrp数据框中。
现在我们将加载一些shapefiles来帮助我们可视化插值结果。
Python
data=gpd.read_file('Taluk_Boundary.shp')
data=data[data['KGISDistri']==20].iloc[0:4]
itrp=gpd.GeoDataFrame(itrp,geometry=gpd.points_from_xy(itrp.lng, itrp.lat))
stn=gpd.GeoDataFrame(stn,geometry=gpd.points_from_xy(stn.lng, stn.lat))
sample=gpd.GeoDataFrame(sample,geometry=gpd.points_from_xy(sample.lng, sample.lat))
sample.crs={'init' :'epsg:4326'}
sample=sample.to_crs(data.crs)
stn.crs={'init' :'epsg:4326'}
stn=stn.to_crs(data.crs)
itrp.crs={'init' :'epsg:4326'}
itrp=itrp.to_crs(data.crs)
数据 "包含班加罗尔市的shapefile。
我们将 "itrp"、"sample "和 "stn "转换为GeoDataFrame来绘制点。
最后，我们为所有新创建的地理数据框架设置坐标参考系统（简称CRS）。这应该与我们的shapefile "数据 "相同。
下面是插值的结果。
Python
ax=data.plot(color='white', edgecolor='black',figsize=(25,30))
itrp.plot(ax=ax, column='PM2.5', markersize=50,figsize=(25,30),legend=True)
sample.plot(ax=ax, marker='o', column='PM2.5', markersize=50,figsize=(25,30),label='SiteName')
for x, y, label in zip(sample.geometry.x, sample.geometry.y, sample.SiteName):
    ax.annotate(label.split(',')[0], xy=(x, y), xytext=(10, 10), textcoords="offset points")
我们也可以使用其他的IDW实现，下面是其中的一个。
Python
from photutils.utils import ShepardIDWInterpolator as idw
vals=idw(sample[['lat','lng']].to_numpy(),sample['PM2.5'].to_numpy())
itrp['PM2.5']=vals(itrp[['lat','lng']])
让我们来绘制并看看结果。
Python
ax=data.plot(color='white', edgecolor='black',figsize=(25,30))
itrp.plot(ax=ax, column='PM2.5', markersize=50,figsize=(25,30),legend=True)
sample.plot(ax=ax, marker='o', column='PM2.5', markersize=50,figsize=(25,30),label='SiteName')
for x, y, label in zip(sample.geometry.x, sample.geometry.y, sample.SiteName):
    ax.annotate(label.split(',')[0], xy=(x, y), xytext=(10, 10), textcoords="offset points")
我们可以看到，两种方法的结果都是一样的。
这就是本教程的结束。在本教程中，我们通过从头开始手动实现IDW和使用内置库来执行IDW，了解了IDW的基本原理。虽然有很多其他的方法来执行插值，但IDW是最容易理解和最强大的方法之一。因此，它是最受欢迎的方法之一。