如何基于区域质心对NetCDF数据进行插值？

问题背景 #

我有一个分辨率为0.25x0.25、包含t_ave变量的NetCDF文件，想要将t_ave值分配给各行政区，但部分小型行政区未被现有网格覆盖。尝试重采样至0.01x0.01分辨率虽可行，但耗时久、输出文件过大，并非高效方案。现使用Python 3.9，希望直接基于行政区质心（来自带标准坐标的CSV文件）对原0.25分辨率的NetCDF网格（或轻度重采样后）进行插值，寻求可行方案。

核心思路 #

跳过全局重采样步骤，直接利用行政区质心的坐标点，从NetCDF原始网格数据中通过空间插值计算对应点的t_ave值，大幅提升处理效率。

实现步骤 #

1. 安装依赖库 #

确保所需库已安装：

pip install pandas geopandas xarray scipy netCDF4

2. 加载数据 #

import pandas as pd
import geopandas as gpd
import xarray as xr
import numpy as np
from scipy.interpolate import griddata

# 加载NetCDF气象数据
nc_dataset = xr.open_dataset("myfile.nc")
# 提取网格的经纬度序列和t_ave变量
lons = nc_dataset.lon.values
lats = nc_dataset.lat.values
t_ave_values = nc_dataset.t_ave.values

# 加载行政区质心CSV
centroids_df = pd.read_csv(
    "/Users/sandor/Documents/pythonProject/gis/centroids_xy_2015.csv",
    delimiter=";"
)
# 转换为GeoDataFrame（可选，用于空间校验）
centroids_gdf = gpd.GeoDataFrame(
    centroids_df,
    geometry=gpd.points_from_xy(centroids_df.x, centroids_df.y),
    crs="EPSG:4326"
)

3. 执行空间插值 #

根据需求选择插值方法，以下是三种常用方案：

# 将NetCDF网格转换为插值所需的点集合
grid_points = [(lon, lat) for lon in lons for lat in lats]
t_ave_flat = t_ave_values.flatten()

# 提取质心的坐标对
centroid_coords = list(zip(centroids_gdf.x, centroids_gdf.y))

# 方法1：最近邻插值（最快，直接取最近网格点的值）
centroids_df["t_ave_nearest"] = griddata(
    grid_points, t_ave_flat, centroid_coords, method="nearest"
)

# 方法2：线性插值（平衡速度与精度，基于邻近点计算）
centroids_df["t_ave_linear"] = griddata(
    grid_points, t_ave_flat, centroid_coords, method="linear"
)

# 方法3：立方插值（精度最高，耗时较长，适合平滑结果）
centroids_df["t_ave_cubic"] = griddata(
    grid_points, t_ave_flat, centroid_coords, method="cubic"
)

4. 关联插值结果与行政区数据 #

如果需要将插值结果绑定到行政区边界数据：

# 加载行政区边界GeoJSON
admin_areas = gpd.read_file("admin_areas.geojson")
admin_areas = admin_areas.explode()

# 通过行政区ID合并插值结果（需确保两边数据有匹配的ID字段）
admin_with_tave = admin_areas.merge(
    centroids_df,
    left_on="admin_id",  # 替换为行政区数据中的唯一ID字段名
    right_on="admin_id",  # 替换为质心CSV中的唯一ID字段名
    how="left"
)

# 导出最终结果
admin_with_tave.to_csv("admin_tave_results.csv", index=False)

可选优化：轻度重采样提升精度 #

如果原始0.25分辨率的插值精度不足，可以先轻度重采样（如0.1x0.1）再插值，平衡效率与精度：

# 生成更密的经纬度序列（示例：0.1分辨率）
new_lons = xr.DataArray(np.arange(6, 19, 0.1), dims="lon")
new_lats = xr.DataArray(np.arange(36, 48, 0.1), dims="lat")

# 轻度重采样NetCDF数据
nc_resampled = nc_dataset.interp(lon=new_lons, lat=new_lats, method="linear")

# 后续插值步骤同前，使用重采样后的lons、lats和t_ave_values即可

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