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

Small change needs to be made to the lonbounds part (data are degrees east), because the longitude value ranges from 0 to 359 in the data, so negative numbers will not work in this case. Also the calculation for latli and latui needs to be switched because the value goes from north to south, 89 to -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] ) )  

Answer 3

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 4

To mirror the response from N1B4, you can also do it on one line with climate data operators (cdo):

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

Thus to loop over a set of file, I would do this in a BASH script, using cdo to process each file and then calling your python script:

#!/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 

I always try and do my file processing "offline" as I find it less prone to error. cdo is an alternative to ncks, I'm not saying it is better, I just find it easier to remember the commands. nco in general is more powerful and I resort to it when cdo can't perform the task I wish to carry out.

Answer 5

Note that this can be accomplished even quicker on the command line using NCO's ncks.

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

Answer 6

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