Repository navigation
Support Two-Dimensional Coordinate Variables #605
Description
The CF Conventions supports the notion of a 2d coordinate variable in the case of irregularly spaced data. An example of this sort of dataset is below. The CF Convention is to add a "coordinates" attribute with a string describing the 2d coordinates.
dimensions:
xc = 128 ;
yc = 64 ;
lev = 18 ;
variables:
float T(lev,yc,xc) ;
T:long_name = "temperature" ;
T:units = "K" ;
T:coordinates = "lon lat" ;
float xc(xc) ;
xc:axis = "X" ;
xc:long_name = "x-coordinate in Cartesian system" ;
xc:units = "m" ;
float yc(yc) ;
yc:axis = "Y" ;
yc:long_name = "y-coordinate in Cartesian system" ;
yc:units = "m" ;
float lev(lev) ;
lev:long_name = "pressure level" ;
lev:units = "hPa" ;
float lon(yc,xc) ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
float lat(yc,xc) ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
I'd like to discuss how we could support this in xray. There motivating application for this is in plotting operations but it may also have application in other grouping and remapping operations (e.g. #324, #475, #486).
One option would just to honor the "coordinates" attr in plotting and use the specified coordinates as the x/y values.
ref: http://cfconventions.org/Data/cf-conventions/cf-conventions-1.6/build/cf-conventions.html#idp5559280
Activity
shoyer commentedon Oct 3, 2015
We do support multi-dimensional coordinates already, and even read the
coordinatesattributes. However, we don't do anything with them, beyond keeping them around: http://xray.readthedocs.org/en/stable/data-structures.html#coordinatesIt's definitely worth seeing if we can make
ds.T.plot.contourf('lat', 'lon')work, though.jhamman commentedon Oct 3, 2015
Okay, somehow I missed that we're already doing this to some extent. I just ran across the code in
conventions.pyand, yes, we don't do anything with them.@clarkfitzg - do you have any thoughts on how this could work?
shoyer commentedon Oct 3, 2015
We could probably modify our
contour/contourfwrapper to accept 2D arrays for X/Y. The matplotlib functions already handle all the necessary logic.clarkfitzg commentedon Oct 3, 2015
Are the
coordinatesnot dimensions of thoseDataArrays? Would it be better to make them dimensions when they get read in?shoyer commentedon Oct 3, 2015
These latitude and longitude coordinates can't be dimensions because they need to be 2D. A good example would be simulation data on some sort of projected grid (like in @jhamman's example above).
jhamman commentedon Oct 3, 2015
@shoyer - the idea of the coordinates convention is that we all ready know the x/y coordinate names. So we shouldn't need to supply them like this:
ds.T.plot.contourf('lat', 'lon'). I'm thinking it would be better to store the coord names coming from https://github.com/xray/xray/blob/master/xray/conventions.py#L882 in each DataArray. Maybe aDataArray.coord_keys = ('xc', 'yc')attribute is in order? Currently, we aren't really using the coordinates attribute in its intended way.In the plotting logic, we should just search to see if
coord_keys is not Nonebefore passing the coordinates to the plotting method.shoyer commentedon Oct 3, 2015
My hesitation with adding
coord_keysis two-fold:coordinatesattribute is meaningful according to CF conventions: "The value of the coordinates attribute is a blank separated list of the names of auxiliary coordinate variables. There is no restriction on the order in which the auxiliary coordinate variables appear in the coordinates attribute string."A better way to handle this might be to write a
xray-cfpackages (see #461 for another use case), which registers the.cfnamespace on Dataset/DataArray objects. Then we could haveds.cf.plot, which automatically plots according to CF conventions, and could even do things like automatically plot with cartopy if relevant. This means less domain specific details in xray, and could be a good model for other similar extension projects (e.g., xgcm).jhamman commentedon Oct 4, 2015
@shoyer - For the most part, I'm in agreement here with you. A few comments (reference
ncdumpof a irregular grid netCDF below).xray.niandnjare dimensions, not coordinates butxrayturns them into coordinates. I know why we do this in practice but, for example,frac(nj, ni)only usesnjandni,ycandxcare it's coordinates. So, you see, there is a bit of a disconnect there. Another way to think aboutxrayindexing is thatiselshould work relative to the dimension index, andselshould work on the coordinates. I know I'm out on a limb here...coordinatesattribute. Not restricting the order of the coordinates is a big deficiency here.coordinateattribute right now is that we use it to determine the coordinate names, and throw the attribute out.I seem to remember bringing this up in the past with our discussions of the calendar attribute and still don't like that we remove attributes that we decode.Especially this once since we don't do much with it.I have two related PRs from today just to show which way we could go with all this: #606, #608.
ncdump of domain file
(xray_dev)Joes-iMac$ ncdump -h domain.lnd.wr50a_ar9v4.100920.nc netcdf domain.lnd.wr50a_ar9v4.100920 { dimensions: n = 56375 ; ni = 275 ; nj = 205 ; nv = 4 ; variables: double xc(nj, ni) ; xc:long_name = "longitude of grid cell center" ; xc:units = "degrees_east" ; xc:bounds = "xv" ; double yc(nj, ni) ; yc:long_name = "latitude of grid cell center" ; yc:units = "degrees_north" ; yc:bounds = "yv" ; double xv(nj, ni, nv) ; xv:long_name = "longitude of grid cell verticies" ; xv:units = "degrees_east" ; double yv(nj, ni, nv) ; yv:long_name = "latitude of grid cell verticies" ; yv:units = "degrees_north" ; int mask(nj, ni) ; mask:long_name = "domain mask" ; mask:note = "unitless" ; mask:coordinates = "xc yc" ; mask:comment = "0 value indicates cell is not active" ; double area(nj, ni) ; area:long_name = "area of grid cell in radians squared" ; area:coordinates = "xc yc" ; area:units = "radian2" ; double frac(nj, ni) ; frac:long_name = "fraction of grid cell that is active" ; frac:coordinates = "xc yc" ; frac:note = "unitless" ; frac:filter1 = "error if frac> 1.0+eps or frac < 0.0-eps; eps = 0.1000000E-10" ; frac:filter2 = "limit frac to [fminval,fmaxval]; fminval= 0.1000000E-02 fmaxval= 1.000000" ; // global attributes: :title = "CCSM domain data:" ; :Conventions = "CF-1.0" ;xray representation of domain file
shoyer commentedon Oct 5, 2015
In the case of calendars, the information ends up on the
encodingattribute. I think you're right though that we currently throw away variable specificcoordinatesinformation -- that should probably also be preserved onencoding.There are two reasons why we do this sort of thing:
xcas a coordinate on the dataset, but it remained in thecoordinatesattribute onmask? Or if we deletedxcentirely?These sort of ambiguous scenarios are the part that worry me. In particular, I worry that it leads us down a path of maintaining arbitrary metadata through operations.
jhamman commentedon Oct 6, 2015
You are right, I had forgotten about this. My original statement has been amended accordingly.
I think we're square on my points 2 and 3. Does my first point make sense?
I'll open an issue regarding the coordinates encoding attribute.
jhamman commentedon Oct 9, 2015
closed via #608 and #610.