I am attempting to get salinity data from one of the ocean models in our region and have run into a roadblock. This is ROMS model output in NetCDF form, but the grid is curvilinear (lon and lat arrays are 2D).
I want to get a subset of the data for the bounding box for the area I am interested in, however I am unsure how to go about doing this. Given a lat/lon/depth for a bounding box of interest, what are the grid cells and vice versa.
I didn't know if using PyDap would make it any easier to get the data of interest, or if there was another toolset I should look at. I saw Rob Hetland's Octant, but it looks like it works on NetCDF files and I want to avoid pulling the files down if possible since they are fairly substantial.
Thanks
Dan
Best Answer
If you want to subset some data from a NetCDF file using a lon/lat bounding box and that NetCDF file is not aligned with east/north, one strategy is to use a point-in-polygon routine and then find the min/max i,j indices of those points to define a subset to extract.
There are several different python packages to deal specifically with netCDF data from ROMS (a few examples are: https://github.com/kshedstrom/pyroms, and https://github.com/hetland/octant) but as far as I know they don't include subsetting by longitude/latitude. So here's an example:
which produces this plot here:
Update: I fixed the namespace issues reported by Luke below (I knew I would get burned by my
ipython notebook --pylab
eventually), and it occurs to me that I should probably post an example of subsampling velocity in ROMS, because it's trickier. Here I'm using theshrink
androt2d
functions from Rob Hetland's Octant package to average the velocities to grid cell centers and rotate for plotting. Basemap isn't necessary here, but that's the example I have lying around, so here it is:Produces the following plot: