netcdf4 extract for subset of lat lon

Question

I would like to extract a spatial subset of a rather large netcdf file. From Loop through netcdf files and run calculations - Python or 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

How do I extract just the subset of netcdf file corresponding to a state (say Iowa). Iowa has following boundary lat lon:

Longitude: 89° 5' W to 96° 31' W

Latitude: 40° 36' N to 43° 30' N

Answer 1

Well this is pretty easy, you have to find the index for the upper and lower bound in latitude and longitude. You can do it by finding the value that is closest to the ones you're looking for.

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] ) )  

Then just subset the variable array.

# Air (time, latitude, longitude) 
airSubset = f.variables['air'][ : , latli:latui , lonli:lonui ] 

• Note, i'm assuming the longitude dimension variable is in degrees east, and the air variable has time, latitude, longitude dimensions.

Answer 2

Favo's answer works (I assume; haven't checked). A more direct and idiomatic way is to use numpy's where function to find the necessary indices.

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]

Answer 3

If you like pandas, then you should think about checking out 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))