使用Xarray从Shapefile区域提取NetCDF温度数据及点位值

依赖库安装 #

先安装所需工具包：

pip install xarray geopandas shapely numpy netCDF4

步骤1：加载数据 #

import xarray as xr
import geopandas as gpd
from shapely.geometry import Point
import numpy as np

# 加载NetCDF数据集
ds = xr.open_dataset("your_netcdf_file.nc")
# 加载Shapefile（自动识别EPSG:4326投影）
gdf = gpd.read_file("your_shapefile.shp")

# 选择目标区域多边形（示例取第一个，可根据名称/ID筛选）
target_polygon = gdf.geometry.iloc[0]
# 如果需要匹配所有多边形，合并为单一几何对象：
# target_polygon = gdf.geometry.unary_union

步骤2：构建网格点几何集 #

NetCDF的lon/lat是二维数组（Y,X），需转换为可用于空间判断的点对象：

# 展平二维经纬度数组
lon_flat = ds.lon.values.flatten()
lat_flat = ds.lat.values.flatten()

# 生成所有网格点的Shapely Point对象
points = [Point(lon, lat) for lon, lat in zip(lon_flat, lat_flat)]

# 创建GeoDataFrame，绑定CRS为EPSG:4326
grid_gdf = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326")

# 添加X/Y维度的索引列（注意：Y维度优先展平，索引顺序要对应）
grid_gdf["X_idx"] = np.tile(np.arange(ds.X.size), ds.Y.size)  # 每个Y对应完整X序列
grid_gdf["Y_idx"] = np.repeat(np.arange(ds.Y.size), ds.X.size)  # 每个X重复对应Y索引

步骤3：筛选区域内的网格点 #

通过空间关系判断，筛选出目标多边形内的网格点：

# 判断每个网格点是否在目标多边形内
in_region_mask = grid_gdf.geometry.within(target_polygon)
# 筛选符合条件的点
filtered_grid = grid_gdf[in_region_mask]

# 获取唯一的X/Y索引（用于后续提取数据）
x_indices = filtered_grid["X_idx"].unique()
y_indices = filtered_grid["Y_idx"].unique()

步骤4：提取区域温度数据 #

利用筛选出的索引，提取区域内的温度数组（保留时间维度）：

# 提取区域内的温度数据（xarray Dataset格式，保留所有维度）
region_temp = ds.temperature.isel(X=x_indices, Y=y_indices)

# 转换为NumPy数组（如果需要）
region_temp_np = region_temp.values
# 输出维度：(time, 筛选后的Y数量, 筛选后的X数量)

步骤5：获取特定点位的温度值 #

方法1：匹配最近网格点（适合快速查询）

# 目标点位的经纬度（EPSG:4326）
target_lon = 11.9
target_lat = 47.95

# 计算所有网格点与目标点的距离（简化为经纬度差的绝对值和）
lon_diff = np.abs(ds.lon.values - target_lon)
lat_diff = np.abs(ds.lat.values - target_lat)
distance = lon_diff + lat_diff

# 找到距离最近的网格点索引
min_pos = np.unravel_index(np.argmin(distance), ds.lon.shape)
y_idx, x_idx = min_pos

# 获取该点所有时间步的温度值
point_temp = ds.temperature.isel(Y=y_idx, X=x_idx).values

方法2：线性插值（适合非网格点的精确查询）

# 对经纬度进行线性插值，获取目标点的温度
interp_temp = ds.temperature.interp(lon=target_lon, lat=target_lat, method="linear")

# 获取所有时间步的插值结果
interp_temp_values = interp_temp.values

注意事项 #

• 如果NetCDF未提供lon/lat二维数组，需通过Lambert投影参数将X/Y投影坐标转换为EPSG:4326经纬度，可使用pyproj库完成转换。
• 若Shapefile的CRS不是EPSG:4326，需先通过gdf.to_crs("EPSG:4326")转换坐标。

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