Python降雨点数据IDW插值：如何按Shapefile边界裁剪结果

尝试用Python绘制降雨点数据的插值图，已导入研究区Shapefile作为边界，但当前代码生成的插值图是包含Shapefile边界的固定矩形。需要修改代码，实现类似ArcMap中仅显示Shapefile边界范围内插值结果的效果。

当前输出的插值图为矩形边界，曾尝试用min_x, min_y, max_x, max_y = study_area.total_bounds限制范围，但无效。

使用的插值绘图代码如下：

import geopandas as gpd
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

# Function to calculate IDW interpolation
def idw_interpolation(obs_points, values, new_points, power=2):
    tree = cKDTree(obs_points)
    distances, indices = tree.query(new_points, k=3)
    weights = 1.0 / distances**power
    weights /= weights.sum(axis=1)[:, np.newaxis]
    interpolated_values = np.sum(weights * values[indices], axis=1)
    return interpolated_values

# Read gauge data from Excel file
gauge_data = pd.read_excel('gauge_data.xlsx')

# Read shapefile of the study area
study_area = gpd.read_file('MountainousRegion.shp')

# Calculate extent of the shapefile
min_x, min_y, max_x, max_y = study_area.total_bounds

# Generate grid points within the study area
num_points = 100  # Number of grid points in each dimension
x = np.linspace(min_x, max_x, num_points)
y = np.linspace(min_y, max_y, num_points)
xx, yy = np.meshgrid(x, y)
grid_points = np.column_stack((xx.ravel(), yy.ravel()))

# Rainfall columns to interpolate
rainfall_columns = ['Gauge', 'Rainfall1', 'Rainfall2']

# Create subplots
fig, axes = plt.subplots(1, len(rainfall_columns), figsize=(7.5, 4.5))

# Collect maximum rainfall value from all columns
max_rainfall = gauge_data[rainfall_columns].max().max()

for i, rainfall_col in enumerate(rainfall_columns):
    # Perform IDW interpolation for current rainfall column
    interpolated_values = idw_interpolation(gauge_data[['Longitude', 'Latitude']].values, gauge_data[rainfall_col].values, grid_points)
    interpolated_grid = interpolated_values.reshape(xx.shape)
    
    # Plot interpolated rainfall within the extent of the shapefile
    im = axes[i].imshow(interpolated_grid, extent=[min_x, max_x, min_y, max_y], origin='lower', cmap='jet', aspect='auto', vmax=max_rainfall)
    
    # Plot shapefile
    study_area.plot(ax=axes[i], facecolor='none', edgecolor='black', linewidth=1)
    
    # Set title and axis labels with Times New Roman font
    axes[i].set_title(f'{rainfall_col}', fontname='Times New Roman', fontsize=12)
    axes[i].set_xlabel('Longitude', fontname='Times New Roman', fontsize=12)
    axes[i].set_ylabel('Latitude', fontname='Times New Roman', fontsize=12)

# Adjust space between subplots
plt.subplots_adjust(wspace=0)

# Plot colorbar
divider = make_axes_locatable(axes[-1])
cax = divider.append_axes("right", size="9%", pad=0.2)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label('Rainfall (mm/day)', fontname='Times New Roman', fontsize=12)

# Save the figure in high resolution (300 DPI)
plt.savefig('interpolation_maps.png', dpi=300)

# Show plot
plt.tight_layout()
plt.show()

核心思路是给研究区外的网格点插值值设置为NaN，Matplotlib绘制时会自动忽略NaN区域，只显示边界内的插值结果。具体修改步骤如下：

1. 将生成的网格点转换为GeoDataFrame，利用GeoPandas的within()方法判断每个点是否在研究区内
2. 对研究区外的点，将插值值替换为NaN
3. 重新生成带NaN的插值网格，再进行绘图

修改后的完整代码：

import geopandas as gpd
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from shapely.geometry import Point

# Function to calculate IDW interpolation
def idw_interpolation(obs_points, values, new_points, power=2):
    tree = cKDTree(obs_points)
    distances, indices = tree.query(new_points, k=3)
    weights = 1.0 / distances**power
    weights /= weights.sum(axis=1)[:, np.newaxis]
    interpolated_values = np.sum(weights * values[indices], axis=1)
    return interpolated_values

# Read gauge data from Excel file
gauge_data = pd.read_excel('gauge_data.xlsx')

# Read shapefile of the study area
study_area = gpd.read_file('MountainousRegion.shp')
# 获取研究区的几何对象（若为多面数据，替换为下方注释的合并操作）
study_polygon = study_area.geometry.iloc[0]
# study_union = study_area.geometry.unary_union  # 多面数据时使用

# Calculate extent of the shapefile
min_x, min_y, max_x, max_y = study_area.total_bounds

# Generate grid points within the study area
num_points = 100  # Number of grid points in each dimension
x = np.linspace(min_x, max_x, num_points)
y = np.linspace(min_y, max_y, num_points)
xx, yy = np.meshgrid(x, y)
grid_points = np.column_stack((xx.ravel(), yy.ravel()))

# 将网格点转为GeoDataFrame，判断是否在研究区内
grid_gdf = gpd.GeoDataFrame(
    geometry=[Point(xy) for xy in grid_points],
    crs=study_area.crs  # 确保坐标系与研究区一致
)
# 标记每个点是否在研究区内，多面数据时替换为grid_gdf.geometry.within(study_union)
grid_gdf['in_study_area'] = grid_gdf.geometry.within(study_polygon)

# Rainfall columns to interpolate
rainfall_columns = ['Gauge', 'Rainfall1', 'Rainfall2']

# Create subplots
fig, axes = plt.subplots(1, len(rainfall_columns), figsize=(7.5, 4.5))

# Collect maximum rainfall value from all columns
max_rainfall = gauge_data[rainfall_columns].max().max()

for i, rainfall_col in enumerate(rainfall_columns):
    # Perform IDW interpolation for current rainfall column
    interpolated_values = idw_interpolation(gauge_data[['Longitude', 'Latitude']].values, gauge_data[rainfall_col].values, grid_points)
    
    # 给研究区外的点设置为NaN
    interpolated_values[~grid_gdf['in_study_area']] = np.nan
    interpolated_grid = interpolated_values.reshape(xx.shape)
    
    # Plot interpolated rainfall within the extent of the shapefile
    im = axes[i].imshow(interpolated_grid, extent=[min_x, max_x, min_y, max_y], origin='lower', cmap='jet', aspect='auto', vmax=max_rainfall)
    
    # Plot shapefile
    study_area.plot(ax=axes[i], facecolor='none', edgecolor='black', linewidth=1)
    
    # Set title and axis labels with Times New Roman font
    axes[i].set_title(f'{rainfall_col}', fontname='Times New Roman', fontsize=12)
    axes[i].set_xlabel('Longitude', fontname='Times New Roman', fontsize=12)
    axes[i].set_ylabel('Latitude', fontname='Times New Roman', fontsize=12)

# Adjust space between subplots
plt.subplots_adjust(wspace=0)

# Plot colorbar
divider = make_axes_locatable(axes[-1])
cax = divider.append_axes("right", size="9%", pad=0.2)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label('Rainfall (mm/day)', fontname='Times New Roman', fontsize=12)

# Save the figure in high resolution (300 DPI)
plt.savefig('interpolation_maps.png', dpi=300)

# Show plot
plt.tight_layout()
plt.show()

关键修改说明 #

• 导入shapely.geometry.Point用于创建网格点的几何对象
• 将网格点转为GeoDataFrame，并匹配研究区的坐标系，确保空间判断准确
• 用within()方法筛选出研究区内的点，将外部点的插值值设为NaN
• 若研究区Shapefile包含多个多边形，需使用study_union = study_area.geometry.unary_union合并所有面，并将空间判断语句改为grid_gdf.geometry.within(study_union)

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