Python 中的空间插值

在地理数据处理和分析中，空间插值是一项非常有用的技术，特别是在我们面对不完整的地理数据时。无论是在气候数据、环境监测数据，还是人口统计数据中，缺失的数据点都会对分析结果造成影响。在这篇文章中，我们将探讨如何在 Python 中利用空间插值技术，特别是反距离加权法（IDW），来推断缺失的空间数据点，并且会以非洲国家的人口密度为示例进行演示。

插值方法？

反距离加权 (IDW) 是一种简单而有效的空间插值方法，其核心思想基于第一地理定律，即“越近越相似”。具体来说，它假设离已知数据点越近的地方，其未知值越接近这些已知点的值。IDW 方法通过对周围点的距离进行加权平均，从而估算出目标点的值。由于它的直观性和易于实现，IDW 在地理信息系统（GIS）领域中被广泛使用。

数据准备

为了演示 IDW 方法，我们使用 GeoPandas 的内置地图数据集“naturalearth_lowres”。该数据集由 Natural Earth 提供，包含了全球各国家的边界和相关统计数据。我们首先将数据集过滤到非洲大陆，并计算每个国家的人口密度。虽然精确的面积计算通常需要将几何数据转换为投影坐标系，但在本文中，我们将简化处理，直接使用现有的经纬度数据来计算。

# Library import
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt

# Load world countries dataset
gdf = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# Filter for Africa
gdf_africa = gdf[gdf['continent'] == 'Africa']

gdf_africa['area'] = gdf_africa.geometry.area
gdf_africa['pop_density'] = gdf_africa['pop_est'] / gdf_africa['area']
gdf_africa.head(5)

在这个数据集中，‘area’代表国家的面积，‘pop_est’是人口估计值，而‘pop_density’则是我们计算得出的人口密度。在接下来的步骤中，我们将模拟一些缺失数据，并使用 IDW 方法来推断这些缺失点的值。

模拟缺失数据

为了测试 IDW 插值的效果，我们将在非洲地图的 GeoDataFrame 中随机选择几个国家，然后人为地将这些国家的人口密度值设为缺失值。

# get a copy of the original data
gdf_africa_missing = gdf_africa.copy()


# Simulating missing data, with 4 randomly picked countries
# By introduce missing values
indices_to_replace = np.asarray([57, 78, 48, 65])  
gdf_africa_missing.loc[indices_to_replace, 'pop_density'] = np.nan

接下来，我们可以将这些数据的原始分布与缺失数据的分布进行可视化对比。

# Show missing data
vmin = gdf_africa_missing.pop_density.min()
vmax = gdf_africa_missing.pop_density.max()

f, ax = plt.subplots(1, 2, figsize=(10, 5))

gdf_africa.plot(column='pop_density', ax=ax[0], cmap='pink', 
                edgecolor='k', vmin=vmin, vmax=vmax, legend=True)

gdf_africa_missing.plot(color='none', ax=ax[1], edgecolor='k')

gdf_africa_missing.plot(column='pop_density', ax=ax[1], cmap='pink',
                edgecolor='k', vmin=vmin, vmax=vmax, legend=True)

ax[0].set_title("Population density - Original", 
                fontsize = 12, pad = 12)
ax[1].set_title("Population density - Missing values", 
                fontsize = 12, pad = 12)

for aax in ax:
    aax.axis('off')

通过这些可视化图，可以清晰地看到原始地图中的人口密度分布以及我们手动设置为缺失的数据点。接下来，我们将使用 IDW 方法对这些缺失数据进行插值。

空间插值

数据准备就绪后，我们可以使用 SciPy 库中的 cKDTree 进行 IDW 插值。cKDTree 是一个快速的 k 近邻搜索算法，用于高效查找距离最近的已知点并计算权重。

# Importing scipy
from scipy.spatial import cKDTree

# Defining e function to perform the IDW
def idw_interpolation(xi, yi, zi, xi_interp, yi_interp, power=2):
    tree = cKDTree(np.c_[xi, yi])
    # k nearest neighbors
    distances, idx = tree.query(np.c_[xi_interp, yi_interp], k=8)  
    weights = 1 / distances**power
    weights /= weights.sum(axis=1)[:, None]
    zi_interp = np.sum(weights * zi[idx], axis=1)
    return zi_interp

# Prepare data for interpolation
gdf_africa_interpol = gdf_africa_missing.copy()
known = gdf_africa_interpol[gdf_africa_interpol['pop_density'].notna()]
unknown = gdf_africa_interpol[gdf_africa_interpol['pop_density'].isna()]

xi = known.geometry.centroid.x.values
yi = known.geometry.centroid.y.values
zi = known['pop_density'].values

xi_interp = unknown.geometry.centroid.x.values
yi_interp = unknown.geometry.centroid.y.values

# Perform IDW interpolation
zi_interp = idw_interpolation(xi, yi, zi, xi_interp, yi_interp)

# Assign interpolated values back to the GeoDataFrame
gdf_africa_interpol.loc[gdf_africa_interpol['pop_density'].isna(),
                        'pop_density'] = zi_interp

此时，插值后的数据已经生成。我们可以将插值后的地图与原始地图并排显示，来观察 IDW 方法的效果。

# Plot the results
f, ax = plt.subplots(1, 2, figsize=(10, 5))

gdf_africa.plot(column='pop_density', ax=ax[0], 
                cmap='pink', edgecolor='k',
                vmin=vmin, vmax=vmax, legend=True)

gdf_africa_interpol.plot(color='none', 
                ax=ax[1], edgecolor='k')

gdf_africa_interpol.plot(column='pop_density', ax=ax[1], 
                cmap='pink', edgecolor='k', 
                         vmin=vmin, vmax=vmax, legend=True)

ax[0].set_title("Population density - Original", \
                fontsize = 12, pad = 12)
ax[1].set_title("Population density - Interpolated", \
                fontsize = 12, pad = 12)

for aax in ax:
    aax.axis('off')

从可视化结果中可以看出，插值地图与原始地图的分布情况非常接近，这表明 IDW 方法在推断缺失数据时表现良好。为了进一步验证这一点，我们将使用相关性分析来定量评估插值的准确性。

相关性分析

最后，为了评估 IDW 插值的效果，我们将插值值与真实的原始值进行比较，并计算两者之间的相关性系数。

# Compare original and interpolated values for missing countries
missing_countries = \
        set(gdf_africa_missing[gdf_africa_missing['pop_density'].isna(\
                                                                )].name)

missing_original = \
        gdf_africa[gdf_africa.name.isin(missing_countries)][['name', \
                                         'pop_density']].set_index('name')

missing_interpol = \
        gdf_africa_interpol[gdf_africa_interpol.name.isin(missing_countries)]\
                [['name', 'pop_density']].rename(columns={'pop_density': \
                'pop_density_interpol'}).set_index('name')

interpol_comparison = missing_original.merge(missing_interpol, 
                                             left_index=True, 
                                             right_index=True)

print(interpol_comparison.corr())

plt.plot(interpol_comparison.pop_density, 
         interpol_comparison.pop_density_interpol, 'o')
plt.xlabel("Original Population Density")
plt.ylabel("Interpolated Population Density")
plt.title("Comparison of Original and Interpolated Values")
plt.show()

结果表明，相关性系数（r）达到了0.94，表明 IDW 插值所得的人口密度与原始数据高度相关。这说明即使使用 IDW 这种相对简单的方法，也能够在缺失数据的情况下生成高质量的空间数据估计。

结论

在本文中，我们展示了如何在 Python 中使用反距离加权法（IDW）进行空间插值，以推断缺失的地理数据。通过非洲的人口密度示例，我们验证了 IDW 方法的有效性。无论是在科研还是实际应用中，当面临缺失的空间数据时，IDW 都是一种强大且易于实现的工具。掌握这种技术，可以帮助我们在数据不完整的情况下，仍然能够获得可靠的分析结果。

求求你点个在看吧，这对我真的很重要