地图插值叠加异常求助：采样村庄数据插值后无法匹配区域地图

我从包含至少105个村庄的区域采样了10个村庄的数据，预测后得到含经度、纬度及预测值的数据集。

我的需求是通过插值将采样数据覆盖到未采样村庄，第一步的插值等高线图生成正常，但将插值结果叠加到区域地图时完全异常，无法覆盖未采样村庄。

第一步正常运行的插值代码 #

from scipy.interpolate import griddata
import numpy as np
import matplotlib.pyplot as plt

# 提取经纬度和预测值列
interpolation_data = decoded_df[['longitude', 'latitude', 'prediction']]

# 删除缺失值行
interpolation_data = interpolation_data.dropna()

# 转换为numpy数组
points = interpolation_data[['longitude', 'latitude']].values
values = interpolation_data['prediction'].values

# 定义插值网格点
grid_points = np.vstack((grid_lon.flatten(), grid_lat.flatten())).T

# 执行线性插值（注：代码注释写了IDW但实际调用的是linear方法）
interpolated_values = griddata(points, values, grid_points, method='linear')

interpolated_values = interpolated_values.reshape(grid_lon.shape)

# 绘制等高线图
plt.contourf(grid_lon, grid_lat, interpolated_values)
plt.colorbar()
plt.scatter(decoded_df['longitude'], decoded_df['latitude'], c=decoded_df['prediction'], cmap='viridis', edgecolors='black')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.title('插值预测结果')
plt.show()

叠加地图的异常代码及问题 #

import geopandas as gpd
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import box

# 读取Babati村庄shapefile
shapefile_path = "Babati Villages/Babati_villages.shp"
gdf_babati = gpd.read_file(shapefile_path)

gdf_bti= gdf_babati[gdf_babati["District_N"] == "Babati"]

# 定义插值网格点
grid_points = np.vstack((grid_lon.flatten(), grid_lat.flatten())).T

# 执行线性插值
interpolated_values = griddata(points, values, grid_points, method='linear')

# 重塑插值结果匹配网格形状
interpolated_values = interpolated_values.reshape(grid_lon.shape)

# 创建Babati区域边界框
bbox = box(gdf_bti.total_bounds[0], gdf_bti.total_bounds[1],
           gdf_bti.total_bounds[2], gdf_bti.total_bounds[3])

# 裁剪插值结果到区域范围（此处存在核心错误）
interpolated_predictions = gpd.clip(interpolated_predictions, bbox)

# 创建子图
fig, ax = plt.subplots(figsize=(10, 10))

# 绘制村庄边界
gdf_bti.plot(ax=ax, facecolor='none', edgecolor='black')

# 绘制插值结果
interpolated_predictions.plot(ax=ax, column='prediction', cmap='viridis', markersize=30, legend=True)

# 添加颜色条
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
interpolated_predictions.plot(ax=cax, column='prediction', cmap='viridis', legend=True, cax=cax)

# 设置标题和标签
ax.set_title('Babati区域插值预测结果')
ax.set_xlabel('经度')
ax.set_ylabel('纬度')

plt.show()

运行上述代码后，插值结果叠加地图完全异常，无法覆盖未采样村庄。

问题分析与修复方案 #

1. 核心错误：未将插值网格转换为GeoDataFrame

你直接对interpolated_values（numpy数组）调用gpd.clip，但interpolated_predictions从未被正确定义为GeoDataFrame，这是导致异常的根本原因。

修复步骤：
将插值网格和对应的预测值转换为带地理信息的GeoDataFrame：

from shapely.geometry import Point

# 生成网格点的几何对象列表
geometry = [Point(xy) for xy in grid_points]
# 创建GeoDataFrame，必须指定坐标系与shapefile一致
interpolated_predictions = gpd.GeoDataFrame(
    {'prediction': interpolated_values.flatten(),
     'geometry': geometry},
    crs=gdf_bti.crs  # 同步shapefile的坐标系
)

2. 坐标系不匹配问题

如果采样数据的经纬度是WGS84（EPSG:4326），而shapefile使用了其他投影坐标系，直接叠加会导致位置错位。

修复：
确保插值数据与shapefile坐标系一致：

# 查看shapefile的坐标系
print(gdf_bti.crs)
# 若插值数据是WGS84，转换为shapefile的坐标系
interpolated_predictions = interpolated_predictions.to_crs(gdf_bti.crs)

3. 针对「覆盖未采样村庄」的优化方案

如果你目标是给每个未采样村庄赋值，而非覆盖整个网格区域，更高效的方式是直接对村庄中心点插值：

# 筛选未采样村庄（需替换为实际的村庄ID列名和采样ID列表）
unsampled_villages = gdf_bti[~gdf_bti['Village_ID'].isin(sampled_village_ids)]
# 获取村庄几何中心点
village_centers = unsampled_villages.geometry.centroid
# 提取中心点坐标
village_coords = np.array([[point.x, point.y] for point in village_centers])
# 对村庄中心点执行插值
village_predictions = griddata(points, values, village_coords, method='linear')
# 将预测值关联到村庄数据
unsampled_villages['prediction'] = village_predictions

# 可视化所有村庄
fig, ax = plt.subplots(figsize=(10,10))
gdf_bti.plot(ax=ax, facecolor='none', edgecolor='black')
# 绘制采样村庄（标记为红色边框）
sampled_villages = gdf_bti[gdf_bti['Village_ID'].isin(sampled_village_ids)]
sampled_villages.plot(ax=ax, column='prediction', cmap='viridis', markersize=50, edgecolor='red', label='采样村庄')
# 绘制未采样村庄
unsampled_villages.plot(ax=ax, column='prediction', cmap='viridis', markersize=30, label='未采样村庄')
plt.legend()
plt.colorbar(ax=ax)
plt.title('Babati区域村庄预测值分布')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.show()

4. 代码冗余清理

删除重复导入的geopandas和matplotlib.pyplot语句，简化代码结构。

内容的提问来源于stack exchange，提问作者LivingstoneM