WRF netcdf file - subset smaller array out of dataset based on coordinate boundaries in python

Question

I have two netcdf files from WRF runs, one with hourly data and the other smaller file has the coordinates (XLAT and XLONG). I am trying to retrieve a subset of the data based on certain coordinates.

An example of one of the variables is temperature 'T2' which has the dimensions (1,1015,1359) being the (time, south_north, west_east), respectively.

The XLAT and XLONG have the same dimensions (1,1015,1359).

There was an identical question asked (please see netcdf4 extract for subset of lat lon), because my lat/long dimensions are a little different the script did not work for me and I haven't been able to figure out why. I tried to change the coordinates into 1D arrays, so that it would be similar to the previous question, but the script does not work and I get an indexing error.

If someone could please help me out that would be awesome! Thanks in advance :)

import numpy as np
from netCDF4 import Dataset  
import matplotlib.pyplot as plt

lons = b.variables['XLONG'][:]
lats = b.variables['XLAT'][:]

lons2d =lons.reshape((1015,1359))
lons1d = lons2d.reshape((1379385))

lats2d =lats.reshape((1015,1359))
lats1d = lats2d.reshape((1379385))

lat_bnds, lon_bnds = [49,53], [-125,-115]
lat_inds = np.where((lats1d > lat_bnds[0]) & (lats1d < lat_bnds[1]))
lon_inds = np.where((lons1d > lon_bnds[0]) & (lons1d < lon_bnds[1]))

T_subset = a.variables['T2'][:,lat_inds,lon_inds]

However I get the following error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-2-0f8890d3b1c5> in <module>()
 25 lon_inds = np.where((lons1d > lon_bnds[0]) & (lons1d < lon_bnds[1]))
 26 
---> 27 T_subset = a.variables['T2'][:,lat_inds,lon_inds]
 28 
 29 

netCDF4/_netCDF4.pyx in      netCDF4._netCDF4.Variable.__getitem__(netCDF4/_netCDF4.c:35672)()

/Users/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/netCDF4/utils.pyc in _StartCountStride(elem, shape, dimensions, grp, datashape, put)
197         # Raise error if multidimensional indexing is used.
198         if ea.ndim > 1:
--> 199             raise IndexError("Index cannot be multidimensional")
200         # set unlim to True if dimension is unlimited and put==True
201         # (called from __setitem__)
IndexError: Index cannot be multidimensional

Answer 1

I see an obvious problem with lat_inds as it has max shape 1015*1359, but You try to use it as an index for latitude, which has size 1015. So IMO You should first find similar values for lat_inds and lon_inds, points which satisfy both lon and lat limits, and then use this array for flattened data. Something like:

uni_ind=numpy.intersect1d(lat_inds,lon_inds)
T_subset=np.ravel(a.variables['T2'])[uni_ind]

Converting array back to 2D may contain some more issues, because I assume Your original data is not in cylindrical coordinates and thus the resulting subset may not be rectangular. This code is not tested, if You share the original data files, I could do that as well.

EDIT: For correct plotting it is easier to use masking, this example should be informative enough.

import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt

b = Dataset('wrfout_conus_constants.nc')
a = Dataset('wrf2d_d01_2010-01-11_000000')

## Data coords
xlong = b.variables['XLONG'][0]
xlat = b.variables['XLAT'][0]
## Data var
temp = a.variables['T2'][0]
## Data bounds
longmax, longmin = -115, -125
latmax, latmin = 53, 49
## Mask coordinates according to bounds
latmask=np.ma.masked_where(xlat<latmin,xlat).mask+np.ma.masked_where(xlat>latmax,xlat).mask
lonmask=np.ma.masked_where(xlong<longmin,xlong).mask+np.ma.masked_where(xlong>longmax,xlat).mask
totmask = lonmask + latmask
## Show mask compared to full domain
plt.pcolormesh(totmask)
## Apply mask to data
temp_masked = np.ma.masked_where(totmask,temp)
## plot masked data
fig=plt.figure()
plt.contourf(temp_masked)
## plot full domain
fig=plt.figure()
plt.contourf(temp)
plt.show()

Answer 2

I'm not sure why it's not working, but I think this does what you want and is cleaner:

import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt

# By indexing at 0 along first dimension, we eliminate the time
# dimension, which only had size 0 anyway.
lons = b.variables['XLONG'][0]
lats = b.variables['XLAT'][0]
temp = a.variables['T2'][0]

lat_bnds, lon_bnds = [49,53], [-125,-115]

# Just AND together all of them and make a big mask
subset = ((lats > lat_bnds[0]) & (lats < lat_bnds[1]) & 
          (lons > lon_bnds[0]) & (lons < lon_bnds[1]))

# Apply mask--should apply to trailing dimensions...I think
T_subset = temp[subset]