如何基于经纬度与时间匹配从3D数组提取shp点值？

不需要转换数据集维度或把shapefile转栅格，直接利用xarray的维度索引/插值功能就能完成提取，步骤如下：

1. 准备数据 #

首先导入依赖库并读取数据：

import geopandas as gpd
import xarray as xr
import pandas as pd

# 读取点shapefile为GeoDataFrame
gdf = gpd.read_file("your_points.shp")
# 转换日期列为datetime格式，确保与xarray的时间维度格式一致
gdf["Date"] = pd.to_datetime(gdf["Date"])

# 读取你的xarray三维数据集（假设维度为time, lat, lon）
ds = xr.open_dataset("your_3d_dataset.nc")

2. 提取对应值 #

根据你的点是否正好落在数据集的网格点上，选择两种方法：

方法A：精确匹配或取最近网格点 #

如果你的点经纬度与数据集的网格点完全重合，或者允许取最近的网格点，使用sel方法：

# 提取指定变量的值（替换为你实际的变量名，比如'temperature'）
extracted = ds["variable_name"].sel(
    time=gdf["Date"],
    lat=gdf["geometry"].y,
    lon=gdf["geometry"].x,
    method="nearest"  # 用"exact"则只匹配完全重合的点，不匹配会报错
)

# 将结果转为数组或合并到GeoDataFrame
gdf["extracted_value"] = extracted.values

方法B：空间插值（点不在网格上时） #

如果点不在数据集的网格点上，需要进行空间插值，使用interp方法：

# 对时间、纬度、经度进行插值
extracted = ds["variable_name"].interp(
    time=gdf["Date"],
    lat=gdf["geometry"].y,
    lon=gdf["geometry"].x
)

# 结果存入GeoDataFrame
gdf["extracted_value"] = extracted.values

3. 输出结果 #

• 若要保存为数组：extracted_array = extracted.values
• 若要保存为新的数据集：extracted_ds = extracted.to_dataset(name="extracted_value")
• 若要保存为带提取值的shapefile：gdf.to_file("points_with_values.shp")

关键说明 #

• 无需转换3D数据集为2D：xarray原生支持多维度的索引操作，直接通过时间、经纬度三个维度定位即可。
• 无需将shapefile转为栅格：GeoDataFrame可以直接提供点的经纬度和时间序列，供xarray进行索引或插值。

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