如何用Python Xarray高效访问OPeNDAP服务器上的稀疏格式网格子集

背景 #

我正在访问存储在OPeNDAP服务器上的GridRad数据，每个时间戳对应独立URL。常规使用Xarray访问OPeNDAP网格时，可先通过切片筛选目标区域，再调用.compute()下载子集，例如高程数据的处理代码：

elevation_url  = 'https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/srtm30plus_v11_land'
elevation_data = xr.open_dataarray(elevation_url)
elevation_data = elevation_data.sel(lon = slice(-105, -100), lat = slice(35, 45)).compute()

GridRad的稀疏存储特性 #

GridRad数据采用稀疏格式存储，打开单时间戳数据集的结构如下：

url = 'https://rda.ucar.edu/thredds/dodsC/files/g/ds841.1/202006/nexrad_3d_v4_2_20200601T000000Z.nc'
radar_data = xr.open_dataset(url)
print(radar_data)

输出结构：

<xarray.Dataset>
Dimensions:         (Longitude: 2832, Latitude: 1248, Altitude: 28,
                     Sweep: 1902, Index: 4288002, time: 1)
Coordinates:
  * Longitude       (Longitude) float64 235.0 235.0 235.1 ... 293.9 294.0 294.0
  * Latitude        (Latitude) float64 24.01 24.03 24.05 ... 49.95 49.97 49.99
  * Altitude        (Altitude) float64 1.0 1.5 2.0 2.5 ... 19.0 20.0 21.0 22.0
  * time            (time) datetime64[ns] 2020-06-01
Dimensions without coordinates: Sweep, Index
Data variables:
    sweeps_merged   (Sweep) |S64 ...
    index           (Index) int32 ...
    Reflectivity    (Index) float32 ...
    wReflectivity   (Index) float32 ...
    SpectrumWidth   (Index) float32 ...
    wSpectrumWidth  (Index) float32 ...
    datehour        (time) int32 ...
    Nradobs         (Altitude, Latitude, Longitude) int8 ...
    Nradecho        (Altitude, Latitude, Longitude) int8 ...
...

其中Index维度记录了变量在Altitude x Latitude x Longitude全量网格中的扁平化索引位置。常规提取方法需要下载完整的变量数据集，效率极低：

refl   = np.zeros(radar_data['Nradobs'].size)*np.nan
refl[radar_data['index']]   = radar_data['Reflectivity']
refl  = refl.reshape(radar_data['Nradobs'].shape)

已完成的筛选工作 #

我已经初步筛选出目标地理区域对应的索引，但不清楚后续如何仅下载该子集：

min_lon = 240
max_lon = 241
min_lat = 37
max_lat = 38
min_alt = 1
max_alt = 2

alt_inds, lat_inds, lon_inds = np.unravel_index(radar_data['index'], radar_data['Nradobs'].shape)

good_lon_inds = radar_data['Longitude'][lon_inds].to_series().between(min_lon, max_lon).to_numpy()
good_lat_inds = radar_data['Latitude'][lat_inds].to_series().between(min_lat, max_lat).to_numpy()
good_alt_inds = radar_data['Altitude'][alt_inds].to_series().between(min_alt, max_alt).to_numpy()

good_inds = good_lon_inds & good_lat_inds & good_alt_inds

问题 #

• 能否仅下载该稀疏格式网格的指定地理子集，避免全量下载？
• 如果可行，如何扩展到多时间戳（多个URL）的批量处理？

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