netcdf4 извлечение для подмножества лат-лон

Я хотел бы извлечь пространственное подмножество довольно большого файла netcdf. Из цикла через файлы netcdf и запускайте вычисления - Python или R

from pylab import *
import netCDF4

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.1989.nc')
# print variables
f.variables.keys()
atemp = f.variables['air'] # TODO: extract spatial subset

Как извлечь только подмножество файла netcdf, соответствующего состоянию (скажем, Айова). Айова имеет следующие границы:

Долгота: от 89° 5' W до 96° 31' W

Широта: 40° 36'с.ш. до 43° 30' с.ш.

Источник

6 ответов

Решение

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

latbounds = [ 40 , 43 ]
lonbounds = [ -96 , -89 ] # degrees east ? 
lats = f.variables['latitude'][:] 
lons = f.variables['longitude'][:]

# latitude lower and upper index
latli = np.argmin( np.abs( lats - latbounds[0] ) )
latui = np.argmin( np.abs( lats - latbounds[1] ) ) 

# longitude lower and upper index
lonli = np.argmin( np.abs( lons - lonbounds[0] ) )
lonui = np.argmin( np.abs( lons - lonbounds[1] ) )

Затем просто подмножество массива переменных.

# Air (time, latitude, longitude) 
airSubset = f.variables['air'][ : , latli:latui , lonli:lonui ]
  • Обратите внимание, я предполагаю, что переменная измерения долготы находится в градусах на восток, а воздушная переменная имеет измерения времени, широты и долготы.
Источник

Ответ Фаво работает (полагаю, не проверял). Более прямым и идиоматическим способом является использование функции numpy's where для поиска необходимых индексов.

lats = f.variables['latitude'][:] 
lons = f.variables['longitude'][:]
lat_bnds, lon_bnds = [40, 43], [-96, -89]

lat_inds = np.where((lats > lat_bnds[0]) & (lats < lat_bnds[1]))
lon_inds = np.where((lons > lon_bnds[0]) & (lons < lon_bnds[1]))

air_subset = f.variables['air'][:,lat_inds,lon_inds]
Источник

Если вы любите панд, вам стоит подумать о проверке xarray.

import xarray as xr

ds = xr.open_dataset('http://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/models/ncep/narr/air.2m.1980.nc',
                     decode_cf=False)
lat_bnds, lon_bnds = [40, 43], [-96, -89]
ds.sel(lat=slice(*lat_bnds), lon=slice(*lon_bnds))
Источник

Обратите внимание, что это может быть выполнено еще быстрее в командной строке, используя Ncks NCO.

ncks -v air -d latitude,40.,43. -d longitude,-89.,-96. infile.nc -O subset_infile.nc

Источник

Чтобы отразить ответ от N1B4, вы также можете сделать это в одной строке с операторами климатических данных (cdo):

cdo sellonlatbox,-96.5,-89,40,43 in.nc out.nc

Таким образом, чтобы перебрать набор файлов, я бы сделал это в скрипте BASH, используя cdo для обработки каждого файла, а затем вызвав ваш скрипт на python:

#!/bin/bash

# pick up a list of files (I'm presuming the loop is over the years)
files=`ls /usgs/data2/rsignell/models/ncep/narr/air.2m.*.nc`

for file in $files ; do 
   # extract the location, I haven't used your exact lat/lons
   cdo sellonlatbox,-96.5,-89,40,43 $file iowa.nc

   # Call your python or R script here to process file iowa.nc
   python script
done

Я всегда стараюсь выполнять обработку файлов в автономном режиме, так как считаю, что она менее подвержена ошибкам. cdo - альтернатива ncks, я не говорю, что это лучше, мне просто легче запомнить команды. nco в целом более мощный, и я прибегаю к нему, когда cdo не может выполнить задачу, которую я хочу выполнить.

Источник

Небольшое изменение необходимо внести в часть lonbounds (данные в градусах на восток), потому что значение долготы в данных колеблется от 0 до 359, поэтому отрицательные числа в этом случае работать не будут. Также необходимо переключать вычисления для latli и latui, потому что значение идет с севера на юг, от 89 до -89.

latbounds = [ 40 , 43 ]
lonbounds = [ 260 , 270 ] # degrees east
lats = f.variables['latitude'][:] 
lons = f.variables['longitude'][:]

# latitude lower and upper index
latli = np.argmin( np.abs( lats - latbounds[1] ) )
latui = np.argmin( np.abs( lats - latbounds[0] ) ) 

# longitude lower and upper index
lonli = np.argmin( np.abs( lons - lonbounds[0] ) )
lonui = np.argmin( np.abs( lons - lonbounds[1] ) )
Источник

Если вы работаете в Linux, с этим очень легко справиться с помощью nctoolkit ( https://nctoolkit.readthedocs.io/en/latest/):

import nctoolkit as nc
data = nc.open_data('/usgs/data2/rsignell/models/ncep/narr/air.2m.1989.nc')
data.clip(lon = [-(96+31/60), -(89+5/6)], lat = [40 + 36/60, 43 + 30/60])
Источник
Другие вопросы по тегам