Нарезка переменной NetCDF с использованием подмножества ограничивающего прямоугольника - PullRequest
Новая
       25

Нарезка переменной NetCDF с использованием подмножества ограничивающего прямоугольника

0 голосов
/ 12 мая 2019

Фон

Я пытаюсь нарезать файл NetCDF, используя ограничивающую рамку в широтах / долях.Соответствующая информация об этом файле приведена ниже (переменные, форма, размеры):

enter image description here

Для большинства ответов здесь и в стандартных уроках это должно быть оченьЭто просто, и моя интерпретация заключается в том, что вы просто находите индексы широты / долготы и разрезаете массив переменных по этим индексам.

Попытка / Код 

 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

Я полагаю, что есть что-то незначительное, что я упускаю / не понимаю в разрезании многомерных массивов, и буду признателен за любую помощь.Я не заинтересован в каких-либо решениях, которые приносят какие-либо другие пакеты или работают вне Python (пожалуйста, никаких ответов на CDO или NCKS!).Спасибо за вашу помощь.

1 Ответ

2 голосов
/ 14 мая 2019

В 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, автоматически нарезаются.

В качестве альтернативы, если вы хотите выбрать поле вручную (используя 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]

...