使用边界框子集切片 NetCDF 变量

Question

背景

我正在尝试使用 lat/lons 的边界框对 NetCDF 文件进行切片。该文件的相关信息如下（变量、形状、尺寸）：

根据此处的大多数答案和标准教程，这应该非常简单，我的解释是您只需找到 lat/lons 的索引并根据这些索引对变量数组进行切片。

Attempt/Code

 def netcdf_worker(nc_file, bbox):
    dataset = Dataset(nc_file)
    for variable in dataset.variables.keys():
        if (variable != 'lat') and (variable != 'lon'):
            var_name = variable
            break

    # Full extent of data
    lats = dataset.variables['lat'][:]
    lons = dataset.variables['lon'][:]

    if bbox:
        lat_bnds = [bbox[0], bbox[2]]  # min lat, max lat
        lon_bnds = [bbox[1], bbox[3]]  # min lon, max lon
        lat_inds = np.where((lats > lat_bnds[0]) & (lats < lat_bnds[1]))
        lon_inds = np.where((lons > lon_bnds[0]) & (lons < lon_bnds[1]))

        var_subset = dataset.variables[var_name][:, lat_inds[0], lon_inds[0]]

        # would also be great to slice the lats and lons too for visualization

问题

当尝试通过上述代码实施在 SO 上列出的其他答案中找到的解决方案时，我遇到了错误：

File "/Users/XXXXXX/Desktop/Viewer/viewer.py", line 41, in netcdf_worker
    var_subset = dataset.variables[var_name][:, lat_inds[0], lon_inds[0]]
  File "netCDF4/_netCDF4.pyx", line 4095, in netCDF4._netCDF4.Variable.__getitem__
  File "/Users/XXXXXX/Viewer/lib/python3.6/site-packages/netCDF4/utils.py", line 242, in _StartCountStride
    ea = np.where(ea < 0, ea + shape[i], ea)
IndexError: tuple index out of range

我相信我 missing/not 了解关于切片多维数组的一些小问题，希望得到任何帮助。我对任何带来任何其他包或在 python 外部运行的解决方案不感兴趣（请不要回答 CDO 或 NCKS！）。谢谢你的帮助。

Answer 1

在Python中，我认为最简单的解决方案是使用xarray。最小示例（使用一些 ERA5 数据）：

import xarray as xr

f = xr.open_dataset('model_fc.nc')

print(f['latitude'].values)  # [52.771 52.471 52.171 51.871 51.571 51.271 50.971]
print(f['longitude'].values) # [3.927 4.227 4.527 4.827 5.127 5.427 5.727]

f2 = f.sel(longitude=slice(4.5, 5.4), latitude=slice(52.45, 51.5))  

print(f2['latitude'].values)  # [52.171 51.871 51.571]
print(f2['longitude'].values) # [4.527 4.827 5.127]

作为示例，我只显示了 latitude 和 longitude 变量，但是 NetCDF 文件中具有 latitude 和 longitude 维度的所有变量都被切片了自动。

或者，如果您想 select 手动框（使用 NetCDF4）：

import netCDF4 as nc4
import numpy as np

f = nc4.Dataset('model_fc.nc')

lat = f.variables['latitude'][:]
lon = f.variables['longitude'][:]

# All indices in bounding box:
where_j = np.where((lon >= 4.5) & (lon <= 5.4))[0]
where_i = np.where((lat >= 51.5) & (lat <= 52.45))[0]

# Start and end+1 indices in each dimension:
i0 = where_i[0]
i1 = where_i[-1]+1

j0 = where_j[0]
j1 = where_j[-1]+1

print(lat[i0:i1])  # [52.171 51.871 51.571]
print(lon[j0:j1])  # [4.527 4.827 5.127]

当然，现在您必须手动对每个数据数组进行切片，例如使用var_slice = var[j0:j1, i0:i1]

使用边界框子集切片 NetCDF 变量

Slicing NetCDF variable using bounding box subset

python

numpy

multidimensional-array

netcdf

netcdf4