GPM DPR不规则条带降水数据0.5°×0.5°保守重格网方法咨询及代码优化疑问

问题核心解答 #

• 你的判断完全正确：最近邻法转换到细网格会破坏降水总量守恒。原因是最近邻仅为每个细网格像素分配距离最近的条带像素值，完全忽略了条带像素与细网格像素的重叠面积权重，导致原始降水总量无法被保留，尤其是在条带边缘区域误差会更明显。
• 正确的思路是跳过中间细网格步骤，直接对原始条带数据应用守恒法重采样到目标0.5°网格，核心是要获取每个DPR条带像素的经纬度边界（而非仅中心坐标）。

关键前提：获取DPR条带像素的边界 #

GPM DPR的L2级产品通常自带每个像素的四个角点经纬度（变量名一般为Longitude_Bounds/Latitude_Bounds，具体取决于你使用的数据集版本），这些边界是守恒法重采样的必要输入。

修正后的实现代码 #

import numpy as np
import pandas as pd
import xarray as xr
from pyresample import SwathDefinition, AreaDefinition

# ══════════════════════════════════════════════════════════════════
# 1. CONFIGURATION
# ══════════════════════════════════════════════════════════════════
lon_min, lon_max = -22.7, 47.0
lat_min, lat_max = 23.75, 61.0
tgt_res = 0.5  # conservative target resolution (degrees)
target_date = "2025-07-15"

# ══════════════════════════════════════════════════════════════════
# 2. LOAD AND PREP DATA (注意提取边界变量)
# ══════════════════════════════════════════════════════════════════
# 假设你的ds里包含中心坐标和边界坐标
ds_day = (
    ds[['PrecipRateNearSurface', 'Latitude', 'Longitude', 
        'Latitude_Bounds', 'Longitude_Bounds', 'Time']]
    .compute()
    .to_dataframe()
    .reset_index()
)
ds_day['date'] = pd.to_datetime(ds_day['Time']).dt.floor('D')

# ══════════════════════════════════════════════════════════════════
# 3. BUILD TARGET COARSE GRID
# ══════════════════════════════════════════════════════════════════
lon_edges_c = np.arange(lon_min, lon_max + tgt_res, tgt_res)
lat_edges_c = np.arange(lat_min, lat_max + tgt_res, tgt_res)
ds_coarse = xr.Dataset(
    coords={
        "lon": (["lon"], (lon_edges_c[:-1] + lon_edges_c[1:]) / 2),
        "lat": (["lat"], (lat_edges_c[:-1] + lat_edges_c[1:]) / 2),
        "lon_b": (["lon_b"], lon_edges_c),
        "lat_b": (["lat_b"], lat_edges_c),
    }
)

# ══════════════════════════════════════════════════════════════════
# 4. LOOP OVER DATES — DIRECT CONSERVATIVE REGRID FROM SWATH
# ══════════════════════════════════════════════════════════════════
daily_grids = {}
for date, group in ds_day.groupby('date'):
    if group.empty:
        continue
    
    # 构建带边界的SwathDefinition
    swath_def = SwathDefinition(
        lons=group['Longitude'].to_numpy(),
        lats=group['Latitude'].to_numpy(),
        lon_bounds=group['Longitude_Bounds'].to_numpy().reshape(-1, 4),
        lat_bounds=group['Latitude_Bounds'].to_numpy().reshape(-1, 4)
    )
    
    # 构建目标区域的AreaDefinition
    n_cols_tgt = int(round((lon_max - lon_min) / tgt_res))
    n_rows_tgt = int(round((lat_max - lat_min) / tgt_res))
    tgt_area_def = AreaDefinition.from_extent(
        'target_area',
        'EPSG:4326',
        (n_rows_tgt, n_cols_tgt),
        [lon_min, lat_min, lon_max, lat_max]
    )
    
    # 直接用pyresample的conservative重采样
    precip_coarse = swath_def.resample(
        tgt_area_def,
        group['PrecipRateNearSurface'].values,
        method='conservative',
        fill_value=0.0
    )
    
    # 转成xarray DataArray以便后续拼接
    da_coarse = xr.DataArray(
        precip_coarse,
        dims=['lat', 'lon'],
        coords={
            'lat': ds_coarse['lat'].values,
            'lon': ds_coarse['lon'].values
        }
    )
    daily_grids[date] = da_coarse

# ══════════════════════════════════════════════════════════════════
# 5. STACK INTO SINGLE DATAARRAY
# ══════════════════════════════════════════════════════════════════
precip_all = xr.concat(
    list(daily_grids.values()),
    dim=pd.DatetimeIndex(list(daily_grids.keys()), name="time")
)

额外说明 #

• 如果你的数据集没有自带边界变量，可以通过DPR的扫描几何参数（比如扫描角、像素大小等）计算每个像素的边界，但这会复杂很多，优先推荐使用自带边界的产品。
• 用pyresample的conservative方法时，会自动计算每个条带像素与目标网格的重叠面积权重，从而严格保证降水总量守恒。
• 如果你更习惯用xesmf，也可以将条带数据转换成带边界的xarray Dataset，然后直接用xe.Regridder的conservative方法，核心同样是提供像素边界信息。

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