[bjlittle]

This PR provides support for the CF ocean sigma over z dimensionless vertical coordinate.

As per the CF Metadata Conventions the standard_name of the dimensonless vertical coordinate is ocean_sigma_z_coordinate. However, the CF community seem to have failed to register this as a valid standard_name 😢

As a follow-on task, I'll communicate this discrepancy and lobby for its definition and inclusion in the official CF Standard Name Table.

I'm also not 100% sure whether the derived coordinate standard_name of sea_surface_height_above_reference_ellipsoid is appropriate. Confirmation welcomed.

This PR has been replaced by #1105.

[bjlittle]

Yeah, yeah @esc24 ... I forgot to update the year (again).

[rsignell-usgs]

Does anyone know a CF-compliant model output that uses "ocean_sigma_over_z" available via OPeNDAP?

[bjlittle]

💤 Yawn! ... this PR is starting to feel quite lonely and sad ... 😢

[esc24]

Anyone care to take up this challenge?

[rsignell-usgs]

I'm ready to try this, but I still don't know of any sample sigma_over_z model output to try it on.

[rsignell-usgs]

Okay, I finally created what I believe to be a CF compliant sigma_over_z ocean model output, and tested out this pull request. If I do this:

cubes = iris.load('http://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/models/ncom/ncom_sigma_over_z.nc')

I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-151-c1613712bdeb> in <module>()
----> 1 cubes = iris.load('http://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/models/ncom/ncom_sigma_over_z.nc')

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/__init__.pyc in load(uris, constraints, callback)
    196 
    197     """
--> 198     return _load_collection(uris, constraints, callback).merged().cubes()
    199 
    200 

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/__init__.pyc in _load_collection(uris, constraints, callback)
    166     try:
    167         cubes = _generate_cubes(uris, callback)
--> 168         result = iris.cube._CubeFilterCollection.from_cubes(cubes, constraints)
    169     except EOFError as e:
    170         raise iris.exceptions.TranslationError(

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/cube.pyc in from_cubes(cubes, constraints)
    133         pairs = [_CubeFilter(constraint) for constraint in constraints]
    134         collection = _CubeFilterCollection(pairs)
--> 135         for cube in cubes:
    136             collection.add_cube(cube)
    137         return collection

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/__init__.pyc in _generate_cubes(uris, callback)
    157         elif scheme in ['http', 'https']:
    158             urls = [':'.join(x) for x in groups]
--> 159             for cube in iris.io.load_http(urls, callback):
    160                 yield cube
    161         else:

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/io/__init__.pyc in load_http(urls, callback)
    198     # Call each iris format handler with the appropriate filenames
    199     for handling_format_spec, fnames in handler_map.iteritems():
--> 200         for cube in handling_format_spec.handler(fnames, callback):
    201             yield cube
    202 

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/fileformats/netcdf.pyc in load_cubes(filenames, callback)
    465                 # Process any associated formula terms and attach
    466                 # the corresponding AuxCoordFactory.
--> 467                 _load_aux_factory(engine, cf, filename, cube)
    468 
    469                 # Perform any user registered callback function.

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/fileformats/netcdf.pyc in _load_aux_factory(engine, cf, filename, cube)
    418         terms_to_var_names = engine.requires['formula_terms']
    419         sigma = coord_from_var_name(terms_to_var_names['sigma'])
--> 420         eta = coord_from_var_name(terms_to_var_names['eta'])
    421         depth = coord_from_var_name(terms_to_var_names['depth'])
    422         depth_c = coord_from_var_name(terms_to_var_names['depth_c'])

/home/local/python27_epd/lib/python2.7/site-packages/Iris-1.5.0_dev-py2.7.egg/iris/fileformats/netcdf.pyc in coord_from_var_name(name)
    401                 return coord
    402         raise ValueError('Unable to find coordinate for variable '
--> 403                          '{!r}'.format(name))
    404 
    405     formula_type = engine.requires.get('formula_type')

ValueError: Unable to find coordinate for variable u'zeta'

@bjlittle: Do you see what is wrong?

Thanks,
Rich

[rsignell-usgs]

@bjlittle or @bblay , any chance you could take a quick look at this?

Thanks,
Rich

[rsignell-usgs]

@bjlittle or @bblay , should I post this problem to the google group?

I'm thinking that this problem:

ValueError: Unable to find coordinate for variable u'zeta'

seems possibly unrelated to this pull request.

[esc24]

@bjlittle has been away at a conference, but he is now back. I'm sure he'll respond shortly.

[bjlittle]

Hi @rsignell-usgs! Okay I'll have a look ...

[bjlittle]

@rsignell-usgs!

Whoop, here it is ...

>>> import iris
>>> cube = iris.load_cube('ncom_sigma_over_z_new.nc')
>>> print cube
sea_water_potential_temperature                  (time: 1; *ANONYMOUS*: 40; *ANONYMOUS*: 1280; *ANONYMOUS*: 2048)
     Dimension coordinates:
          time                                             x               -                -                  -
     Auxiliary coordinates:
          sea_surface_height                               x               -                x                  x
          siglay                                           -               x                -                  -
          zlay                                             -               x                -                  -
          angle                                            -               -                x                  x
          depth                                            -               -                x                  x
          dx                                               -               -                x                  x
          dy                                               -               -                x                  x
          latitude                                         -               -                x                  x
          longitude                                        -               -                x                  x
     Derived coordinates:
          sea_surface_height_above_reference_ellipsoid     x               x                x                  x
     Scalar coordinates:
          depth_c: 137.869 m
          nsiglay: 19
     Attributes:
          Conventions: CF-1.6
>>> 

I had to make some minor changes to your NetCDF file in order to get Iris to play ball. Here's the CDL ncdump of my tweaked version of your file ...

netcdf ncom_sigma_over_z_new {
dimensions:
    level = 40 ;
    time = 1 ;
    lon = 2048 ;
    lat = 1280 ;
variables:
    float time(time) ;
        time:units = "days since 2013-01-01 00:00" ;
        time:standard_name = "time" ;
    float zlay(level) ;
        zlay:units = "m" ;
    float siglay(level) ;
        siglay:standard_name = "ocean_sigma_z_coordinate" ;
        siglay:positive = "up" ;
        siglay:formula_terms = "sigma: siglay eta: zeta depth: depth depth_c: depth_c nsigma: nsiglay zlev: zlay" ;
    float depth_c ;
        depth_c:units = "m" ;
    int nsiglay ;
    float lon(lat, lon) ;
        lon:units = "degrees_east" ;
        lon:standard_name = "longitude" ;
    float lat(lat, lon) ;
        lat:units = "degrees_north" ;
        lat:standard_name = "latitude" ;
    float depth(lat, lon) ;
        depth:units = "m" ;
        depth:standard_name = "depth" ;
        depth:coordinates = "lat lon" ;
    float angle(lat, lon) ;
    float dx(lat, lon) ;
    float dy(lat, lon) ;
    float temp(time, level, lat, lon) ;
        temp:units = "degC" ;
        temp:standard_name = "sea_water_potential_temperature" ;
        temp:coordinates = "siglay zeta depth depth_c nsiglay zlay lat lon angle dx dy" ;
    float zeta(time, lat, lon) ;
        zeta:units = "m" ;
        zeta:standard_name = "sea_surface_height" ;

// global attributes:
        :Conventions = "CF-1.6" ;
}

Please note that I made the following changes:

• zeta - removed the coordinates attribute, it's not necessary as this variable is a CF auxiliary coordinate and not a CF data variable.
• depth_c - requires the appropriate units attribute.
• zlay - requires the appropriate units attribute.
• temp - this is the CF data variable and you need to explicitly state all its associated auxiliary coordinates in the coordinates attribute. This was the root of your problem 😉

And cube.coord('sea_surface_height_above_reference_ellipsoid') does indeed return a 4-d coordinate, hoo-rah! ... can you easily determine/verify whether the derived coordinate points are sensible ?

[bjlittle]

Hmmm it's not crystal clear to me from the CF standard whether all the formula_terms variables should be explicitly stated in the coordinates attribute of the associated data variable. I've just assumed that's the case. Clearly @rsignell-usgs naturally assumed otherwise.

So is it sufficient that only the variable name of the ocean_sigma_z_coordinate is stated in the coordinates attribute i.e. siglay in this case, and all the formula_terms variables are implied by association ?

I prefer our simple, yet explicit approach ... should Iris gracefully cater for both scenarios ?

Thoughts/guidance/clarity anyone @esc24, @rhattersley, @marqh, @rsignell-usgs ?

[esc24]

So is it sufficient that only the variable name of the ocean_sigma_z_coordinate is stated in the coordinates attribute i.e. siglay in this case, and all the formula_terms variables are implied by association ?

I can't offer an authoritative opinion with respect to CF, but my gut says the formula_terms variables should be implied by association and should not need to be listed twice. Having said that, if they are there twice it should still work.

[bjlittle]

Although I like to keep life simple, I tend to agree with you @esc24 ... implied through association would be better. Particularly from the user experience point of view ...

And this isn't the first time I've seen Iris choke with regards to this issue with formula_terms and coordinates.

[rsignell-usgs]

@bjlittle, I've been using Dimensionless Vertical Coordinate data (http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.html#dimensionless-vertical-coordinate) for five years or so, but always using the NetCDF-Java library to compute the vertical coordinates (e.g. z positions at a specified time step). I know that the NetCDF-Java library does not require all the formula_terms variables to appear in the coordinates attribute: only the vertical coordinate variable that contains the formula_terms attribute.

I'm also curious about the need to remove coordinates from the zeta variable, because although it's used in the computation of the vertical coordinate (independent variable), it also is a dependent variable (people often want to access and plot the water levels).

[rsignell-usgs]

Okay @bjlittle, I checked the vertical coordinate z values (points) and the values in the range determined by sigma are all zero, which is a problem. Check out what I did in this IPython Notebook:
http://nbviewer.ipython.org/urls/raw.github.com/rsignell-usgs/ipython-notebooks/master/ncom_coords.ipynb

[bjlittle]

Hi @rsignell-usgs!

Okay that's great. Thanks for taking the time to analyse the data. I'll have a dig to figure out what's going wrong in my PR.

Also, you shouldn't need to remove the coordinates attribute from the zeta variable. It's just that it's usage didn't make sense to me in the context of temp being the data variable; that's why I removed it. Adding it back in shouldn't make a difference here (well, at least I hope not!).

So let me resolve this zeroing issue of the data. Once done, I'd love to get this PR merged now that it's got some momentum again. I'll nudge another core developer.

With regards to the formula_terms issue, I think that this should be addressed in a follow-up PR. I'll make that happen 😉

[bjlittle]

@rsignell-usgs ... okay, I think I've ironed out the wrinkle. Phew! Completely my mistake. I misunderstood (i.e. didn't engage brain) the min(depth_c, depth(j, i)) part of the equation. Believe it or not I was calculating a global minimum between depth_c and depth! Hence a big, fat array or nothingness.

Your IPython notebook really helped to clarify the situation. So I'm completely indebted to your thoroughness.

Does this look any better @rsignell-usgs ?

I just need to re-check the associated unit tests given this change. Once I've done that I'll push the update to the PR so that you can give it further examination.

Cheers 😉

[bjlittle]

@rsignell-usgs fancy trying the latest cut of this PR with your data set ?

I'm keen to get this merged ...

[rsignell-usgs]

Ah, sorry, I didn't realize there was something for me to try. Thanks for the jingle.

I just tested and it looks good!
http://nbviewer.ipython.org/urls/raw.github.com/rsignell-usgs/ipython-notebooks/master/ncom_coords.ipynb

[bjlittle]

No worries!

Thanks @rsignell-usgs for taking the time to kick the tyres on this one. Very much appreciated.

Okay, for me that's good enough evidence that this PR is fit for purpose.

I'm going to re-base this PR against Iris upstream master, since this PR is now pretty stale. Then I need to grab the attention of another core developer to review and merge .... any takers @esc24 or @rhattersley ?

Once done, I'll look at getting Iris to automatically associate the formula_terms variables as (auxiliary) coordinates of a cube in a separate follow-on PR.

Then you never know ... I might get time to implement some of the easier ocean dimensionless vertical coordinates ...

Ocean sigma coordinate seems like the next natural target. What's at the top of your hit-list @rsignell-usgs ?

[rhattersley]

any takers @esc24 or @rhattersley ?

NB. We're both going to be at FOSS4G next week (and getting ready this week) so are likely to be out of the loop for a little while. Perhaps @bblay or @pelson might have a chance in the meantime.

[bjlittle]

Textbook body swerve ... lovely manoeuvre @rhattersley 😉

[rsignell-usgs]

The believe that the most-used vertical coordinates in the ocean are:
ocean_s_coordinate_g1 (used by most ROMS models > 2010)
ocean_s_coordinate (used by ROMS models < 2010)
ocean_s_coordinate_g2 ( used by some ROMS models > 2010)
ocean_sigma_coordinate (used by POM and ECOM models)

[bjlittle]

Re-based and ready for review.

Any takers @bblay or @pelson ?

[bjlittle]

@rsignell-usgs are ocean_s_coordinate_g1 and ocean_s_coordinate_g2 defined and accepted for CF v1.7 ?

[rsignell-usgs]

yes,
The big man himself (Jonathan Gregory) says so.
See last comment at:
https://cf-pcmdi.llnl.gov/trac/ticket/93

[rsignell-usgs]

And I just received further confirmation from Jonathan this morning. See the CF_trac_summary link below.

"There are a lot of tickets which have been agreed and should therefore be in
the next version of the CF standard. They are all the ones shown in green in
http://www.met.reading.ac.uk/~david/cf_trac_summary.html, including ticket 93.

Making a new version of CF depends on Jeff Painter or a colleague at PCMDI
having time to do it, I presume. I agree we need a new one, because there are
a lot of concluded tickets, and it's two years since CF 1.6 was published. We
do not close the tickets until the new version of CF is published which
includes them, because the process isn't finished until then. In CF 1.7, the
history will be updated to list all the tickets which were enacted in it.

However, even without a new version of CF, it can be assumed that ticket 93
(and all the other concluded ones) will be part of CF 1.7. Software developers
can act on that assumption."

[bjlittle]

Excellent, thanks @rsignell-usgs !

The Google Groups thread https://groups.google.com/forum/#!topic/scitools-iris/D4qjnKa7rB0 provides sufficient test data coverage for these ocean models. I'm assuming that you'd be happy to beta-test as and when we manage to generate the associated PRs? It's proved invaluable with this PR, so thanks again.

[pelson]

I don't think this needs a constant. Define it inline.

[pelson]

@bjlittle - I support this PR as a valuable feature, but its current implementation is worryingly long. I'm keen to reduce the number of lines that this PR adds to Iris - there are a lot of redundant checks in both the construct/update methods and in the tests themselves which I think are ripe for factorisation and/or deletion. Remember: "it is easier to ask forgiveness than permission" 😉

[rsignell-usgs]

@bjlittle , what does this mean? You come back later with a different pull request that is "leaner and meaner"?

[bjlittle]

@rsignell-usgs don't worry ... just taking it off to the side so that I can refactor.

[rsignell-usgs]

@bjlittle, have there been any developments here?

[bjlittle]

Thanks for the nudge @rsignell-usgs ... I've been shamelessly distracted. I'll blow the dust off this PR ...

[esc24]

@bjlittle - take a look at the changes I recently made to the hybrid pressure factory. I pulled out a lot of the common verification code into a static _check_dependencies() method.

[rsignell-usgs]

Where do we stand now?
We need this capability to create a tool that does underwater virtual vehicle transects through ocean model data.

[bjlittle]

@rsignell-usgs We're standing no more ... #1105 has already been merged! #1105 is the refactored PR based on lessons learnt from this PR.

Yup, ocean sigma over z is available from upstream/master right now ... so now go build your cool tool!

Any problems, then give me a jingle ...

[rsignell-usgs]

Awesome!