geopandas 空间插值
时间: 2025-02-19 18:08:32 浏览: 49
### 使用 GeoPandas 进行空间插值
在地理信息系统 (GIS) 应用中,空间插值是一种常用的技术,用于估计未知位置的数据值。GeoPandas 结合其他 Python 库可以实现多种空间插值方法。
#### 数据准备
首先,确保安装必要的库:
```bash
pip install geopandas pandas numpy scipy scikit-learn folium
```
加载所需的模块并读取数据:
```python
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
import pandas as pd
from sklearn.gaussian_process.kernels import RBF
from pykrige.ok import OrdinaryKriging
from scipy.interpolate import griddata
from sklearn.neighbors import KNeighborsRegressor
```
#### 处理缺失值
对于存在缺失值的情况,可以通过模拟的方式引入一些缺失值来测试插值效果[^2]:
```python
# 加载原始数据
gdf_africa = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
gdf_africa = gdf_africa[gdf_africa['continent'] == 'Africa']
# 模拟缺失人口密度数据
indices_to_replace = [57, 78, 48, 65]
gdf_africa_missing = gdf_africa.copy()
gdf_africa_missing.loc[indices_to_replace, 'pop_density'] = np.nan
```
#### 空间插值方法
##### 1. IDW 插值法
IDW(Inverse Distance Weighted)是最常用的插值算法之一。它基于距离权重计算目标点的属性值。
```python
def idw_interpolation(points, values, xi, yi, power=2):
dists = np.sqrt((xi[:, None] - points[:, 0])**2 + (yi[:, None] - points[:, 1])**2)**power
weights = 1 / dists
sum_weights = np.sum(weights, axis=1)
# 防止除零错误
sum_weights[sum_weights==0] = 1e-10
zi = np.dot(weights.T, values) / sum_weights
return zi.reshape(xi.shape)
points = [(row.geometry.x, row.geometry.y) for _, row in gdf_africa.iterrows()]
values = gdf_africa['pop_density'].dropna().tolist()
grid_x, grid_y = np.mgrid[gdf_africa.total_bounds[:2], gdf_africa.total_bounds[2:]]
interpolated_values_idw = idw_interpolation(np.array(points), np.array(values), grid_x.flatten(), grid_y.flatten())
```
##### 2. Kriging 法
克里金插值考虑了样本之间的自相关性,适用于具有较强空间结构的数据集。
```python
coords = list(zip(*[(p.x, p.y) for p in gdf_africa.geometry]))
z = gdf_africa.pop_density.dropna().to_numpy()
OK = OrdinaryKriging(
coords,
z,
variogram_model="linear",
)
grid_lon = np.linspace(min(x_coords), max(x_coords))
grid_lat = np.linspace(min(y_coords), max(y_coords))
zi_krige, _ = OK.execute("grid", grid_lon, grid_lat)
```
##### 3. GridData 方法
`scipy.griddata` 提供了一种简单的网格化方式来进行插值。
```python
xi = np.linspace(min(x_coords), max(x_coords))
yi = np.linspace(min(y_coords), max(y_coords))
XI, YI = np.meshgrid(xi, yi)
ZI = griddata(points, values, (XI, YI), method='cubic')
```
#### 可视化结果
利用 Folium 或 Matplotlib 来展示插值后的地图图像。
```python
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.contourf(XI, YI, ZI, alpha=.75, cmap='viridis')
plt.colorbar(label='Population Density')
plt.title('Interpolated Population Density Map using Cubic Method')
plt.show()
```
阅读全文
相关推荐



















