Accessing/Plotting Curvilinear Grids

[rsignell]

The USGS COAWST (Coupled Ocean, Atmosphere, Wave and Sediment Transport) forecast model has a curvilinear grid in the horizontal dimension and a ocean_s_coordinate coordinate in the vertical dimension.
The data is CF-Compliant and loads up fine in NetCDF-Java based apps like Unidata's IDV and the NCTOOLBOX for Matlab.  

I'd love to get this going in Iris.   

My first try is at: http://nbviewer.ipython.org/5284711

It seems that Iris is correctly finding the 2D coordinate variables for lat and lon, which is cool.

I'm not sure how to get Iris to plot using these coordinates, however.

Also, it seems that constraints do not work for 2D coordinates.  

Is that expected behaviour for now?

Thanks,
Rich

[pp-mo]

This look like a very interesting application !
Just trying to give some feedback for now...

>>  Is that expected behaviour for now?

Basically, yes.
At present 2D coordinates can exist, but there isn't a whole lot you can do with them in terms of analysis or data manipulation.
In Iris, there is at present essentially no special support for handling auxiliary horizontal coordinates, so they are just 'associated values' with each gridpoint.

>>  I'm not sure how to get Iris to plot using these coordinates, however.

With data like this, although the data is structured as a 2D array, the exact 'meaning' of the array dimensions in terms of a geospatial projection is not available.
So, although you have a full set of (true-)lats and lons for the gridpoints, this doesn't let you precisely calculate the outlines of gridboxes to make a filled colour plot (i.e. like pcolor).
Because of this, directly plotting the data as filled patches over a map is not possible.
However, contouring can be done, using matplotlib contour(xx, yy, data), with xx,yy as full 2d arrays matching the data shape.

It is also possible to 'guess' gridboxes for a filled plot by calculating arrays of x and y values interpolated halfway between the given points, and extrapolated "half a step" beyond at the edges
This can then be used effectively with pcolormesh to draw a patch for each gridcell, but there may be 'awkward' cases for which the results are not good.

>>  Also, it seems that constraints do not work for 2D coordinates. 

Constraints basically provide a link between coordinate values and array indice.  For now at least, this is only implemented for one-dimensional coordinates
 -- i.e. if a coordinate is mapped to one array dimension, you can then convert a constraint on coordinate values into an indexing on the dimension in the obvious way.

The 'extract' operation also works by converting coordinate values to index ranges, so that doesn't extend to multidimensional coordinates either.
However, you can define a region of interest with a mask calculated from the lat+lon values.
You can also, with some effort, calculate the rectangular boundary of a required region and extract that.
Perhaps this approach would be useful? (I got this code to work with something like your data) ...

data = cube.data[0,0]
lons = cube.coord('longitude').points
lats = cube.coord('latitude').points

def minmax(v):
    return np.min(v), np.max(v)

# calculate where in the data
lons_range = -71.5, -65.0
lats_range = 39.5, 46.0
lonslo, lonshi = minmax(lons_range)
latslo, latshi = minmax(lats_range)
inregion = np.logical_and(np.logical_and(lons > lonslo, lons < lonshi),
                          np.logical_and(lats > latslo, lats < latshi))

# mask out all data that isn't in our chosen area -- or is NAN
data_selected = np.ma.array(data,
                            mask=np.logical_or(np.isnan(data),
                                               np.logical_not(inregion)))

# extract the rectangular subarray containing the whole valid region ...
region_inds = np.where(inregion)
imin, imax = minmax(region_inds[0])
jmin, jmax = minmax(region_inds[1])
subcube = cube[imin:imax+1, jmin:jmax+1]

[rsignell]

pp-mo,

Thanks - that works!  
I used your code and was able to extract the region with curvliinear coords contained by a bounding box and plot it using Cartopy:

http://nbviewer.ipython.org/5297644

There is a minor problem with placement of tick labels being clipped in IPython Notebook and the title overwrites the tick lables, but other than that, it looks great!

-Rich

P.S. The clipping of tick labels reminds me of this: https://github.com/matplotlib/matplotlib/issues/688

[pp-mo]

Hi Rich

That looks nice !
I do like the notebook approach.

I hope you're using pcolor with caution, though ?
I tried that yesterday and was reminded of this problem..

If you pass NxM arrays of X, Y and DATA to pcolormesh, it treats the X and Y as defining the corners of N-1xM-1 quadrilateral cells.
It then colours each cell according to the value aligned with its lower-left corner - and the whole last row + last column of DATA are ignored.
I.E.  points ([X[i],Y[j]), (X[i+1],Y[j]), (X[i],Y[j+1]), (X[i+1],Y[j+1])) define a region colored by DATA[i,j]

(As usual, all this is there in the documentation if you look really close!)

So each cell is offset from the actual point that defines its colour.  What it doesn't do, which is what you'd really like, is to paint cells surrounding each given point (like Iris guess_bounds might provide). 
You can achieve that by interpolating/extrapolating the points arrays, but it's only an approximation and may not work in tricky cases (e.g. near poles and across longitude seams)

Regards
Patrick

[Phil Elson]

To visualise this data you could also use iris.analysis.cartography.project (http://scitools.org.uk/iris/docs/latest/iris/iris/analysis/cartography.html?highlight=project#iris.analysis.cartography.project).

We do this for data on the ORCA grids from NEMO (http://www.nemo-ocean.eu/About-NEMO). The output is only useful for visualisation purposes, but might be of use to you here.

> P.S. The clipping of tick labels reminds me of this: https://github.com/matplotlib/matplotlib/issues/688

It does look like the same problem. I've recently fixed this in matplotlib (https://github.com/matplotlib/matplotlib/pull/1448).

HTH