如何使用地理坐标裁剪WRADLIB中的RADOLAN雷达数据？

使用wradlib读取德国气象局（DWD）的RADOLAN雷达数据集，需要通过指定经纬度范围裁剪数据。现有代码读取数据后，DataArray仅包含x、y坐标，已尝试创建投影和获取网格，但不清楚如何完成投影转换与数据裁剪。

现有代码 #

import wradlib as wrl

rv_file = "DE1200_RV2201271935_000"
ds = wrl.io.open_radolan_dataset(rv_file)
data = ds.RV

print(data)

输出结果 #

<xarray.DataArray 'RV' (y: 1200, x: 1100)>
[1320000 values with dtype=float32]
Coordinates:
  * y        (y) float64 -4.809e+03 -4.808e+03 ... -3.611e+03 -3.61e+03
  * x        (x) float64 -543.5 -542.5 -541.5 -540.5 ... 552.5 553.5 554.5 555.5
Attributes:
    valid_min:      0
    valid_max:      4095
    standard_name:  rainfall_rate
    long_name:      RV
    unit:           mm h-1

已尝试的代码片段 #

proj_stereo = wrl.georef.create_osr("dwd-radolan")
proj_wgs = osr.SpatialReference()
proj_wgs.ImportFromEPSG(4326)

radolan_grid_xy = wrl.georef.get_radolan_grid(1200, 1100)

按以下步骤完成投影转换和数据裁剪：

1. 导入完整依赖模块 #

确保导入GDAL的osr模块（wradlib依赖该模块处理投影）：

import wradlib as wrl
from osgeo import osr

2. 转换RADOLAN网格为经纬度坐标 #

RADOLAN的x/y是以德国诺德韦克为原点的极射投影坐标（单位km），需转换为WGS84经纬度：

# 创建投影对象
proj_stereo = wrl.georef.create_osr("dwd-radolan")
proj_wgs = osr.SpatialReference()
proj_wgs.ImportFromEPSG(4326)

# 获取RADOLAN网格的扁平化x/y坐标
radolan_grid_xy = wrl.georef.get_radolan_grid(1200, 1100)
# 投影转换为经纬度（lon, lat）
lonlat = wrl.georef.reproject(radolan_grid_xy, projection_source=proj_stereo, projection_target=proj_wgs)
# 重塑为与数据匹配的二维网格
lon_grid = lonlat[:, 0].reshape(1200, 1100)
lat_grid = lonlat[:, 1].reshape(1200, 1100)

3. 为DataArray添加经纬度坐标 #

将生成的经纬度网格绑定到原DataArray，方便后续筛选：

data = data.assign_coords(lon=(["y", "x"], lon_grid), lat=(["y", "x"], lat_grid))

4. 按指定经纬度范围裁剪数据 #

自定义经纬度范围，通过布尔索引筛选目标区域：

# 示例裁剪范围：东经8-10度，北纬49-51度
lon_min, lon_max = 8.0, 10.0
lat_min, lat_max = 49.0, 51.0

# 生成掩码：标记需要保留的网格点
mask = (data.lon >= lon_min) & (data.lon <= lon_max) & (data.lat >= lat_min) & (data.lat <= lat_max)

# 应用掩码并移除空维度
cropped_data = data.where(mask, drop=True)

5. 验证裁剪结果 #

打印裁剪后的DataArray，确认坐标范围和数据有效性：

print(cropped_data)

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