如何在行政单元多边形内提取NetCDF温度数据

核心思路 #

放弃 bounding box 这种粗筛选方式，直接通过**空间掩码（Spatial Masking）**精准匹配每个行政单元多边形内的NetCDF像素，再完成逐日统计。推荐使用Python工具链：xarray（处理NetCDF时序数据）+ geopandas（读写shapefile）+ rioxarray（给xarray数据添加空间参考，支持空间裁剪）。

1. 环境依赖安装 #

先安装必备库：

pip install xarray geopandas rioxarray netCDF4 rasterio

2. 数据读取与预处理 #

先统一空间参考（CRS），确保shapefile和NetCDF数据的坐标系一致：

import geopandas as gpd
import rioxarray as rxr

# 读取行政单元shapefile
admin_shp = gpd.read_file("path/to/your/admin_units.shp")

# 读取NetCDF温度数据（假设变量名为`t2m`，对应2米气温）
# 自动加载空间参考，若NetCDF无CRS信息需手动指定（比如WGS84：EPSG:4326）
temp_data = rxr.open_dataset("path/to/your/global_temp.nc", variable=["t2m"])
if temp_data.rio.crs is None:
    temp_data = temp_data.rio.set_crs("EPSG:4326")

# 统一CRS：将shapefile转换为和NetCDF一致的坐标系
if admin_shp.crs != temp_data.rio.crs:
    admin_shp = admin_shp.to_crs(temp_data.rio.crs)

3. 遍历多边形提取并统计 #

对每个行政单元生成专属掩码，筛选内部像素后完成温度区间天数统计：

import pandas as pd

# 定义目标温度区间（示例：0℃~25℃，注意NetCDF温度单位，若为开尔文需转℃：t2m - 273.15）
temp_low = 0
temp_high = 25

# 存储统计结果的列表
results = []

# 遍历每个行政单元
for idx, row in admin_shp.iterrows():
    admin_name = row["admin_name"]  # 替换为你的shapefile中行政单元名称的字段
    admin_poly = row["geometry"]
    
    # 用多边形裁剪NetCDF数据，仅保留内部像素
    masked_temp = temp_data.rio.clip([admin_poly], drop=False)
    
    # 计算逐日平均温度（也可直接统计每个像素的达标天数后聚合，按需调整）
    daily_mean = masked_temp["t2m"].mean(dim=["x", "y"]) - 273.15
    
    # 统计温度落在目标区间内的天数
    target_days = ((daily_mean >= temp_low) & (daily_mean <= temp_high)).sum().item()
    
    # 存入结果
    results.append({
        "行政单元名称": admin_name,
        f"{temp_low}℃~{temp_high}℃天数": target_days,
        "总统计天数": len(daily_mean)
    })

# 导出结果为CSV
result_df = pd.DataFrame(results)
result_df.to_csv("admin_temp_stats.csv", index=False)

4. 性能优化建议 #

• 并行处理：行政单元数量极多时，用joblib并行加速：

  from joblib import Parallel, delayed
  
  def process_admin(row):
      admin_name = row["admin_name"]
      admin_poly = row["geometry"]
      masked_temp = temp_data.rio.clip([admin_poly], drop=False)
      daily_mean = masked_temp["t2m"].mean(dim=["x", "y"]) - 273.15
      target_days = ((daily_mean >= temp_low) & (daily_mean <= temp_high)).sum().item()
      return {"行政单元名称": admin_name, "目标区间天数": target_days}
  
  # 用所有CPU核心并行处理
  results = Parallel(n_jobs=-1)(delayed(process_admin)(row) for idx, row in admin_shp.iterrows())
  

• 内存优化：NetCDF数据过大时，用dask分块读取避免内存溢出：

  temp_data = rxr.open_dataset("path/to/global_temp.nc", variable=["t2m"], chunks={"time": 30})
  

• 加权统计：若需要考虑像素面积权重（大像素占比更高），可通过rasterio.features.rasterize生成带面积权重的掩码，再用xarray的weighted方法计算。

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