Best way to interoperably return the closest time slice

[rsignell]

Iris gang,

I'm trying to demonstrate how the CF conventions lead to model interoperability, so I'd like to show a list of datasets where the users specifies a variable, and is able to extract the closest time slice without model-specific code. 

In the 1st dataset below, the time coordinate variable is "ocean_time", and doesn't have a long_name or standard_name = 'time', thus fails on cube.coord(time).  Yet cube.coord(axis='T') works.

In the 2nd dataset below, the time coordinate variable is "time",so cube.coord(time) works, but cube.coord(axis='T') doesn't work (it complains about finding two time coordinates)

So in the end, this was the workaround I came up with.  There has to be a better way, right?

import datetime as dt
import iris

urls=['http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2009_da/his',
      'http://geoport.whoi.edu/thredds/dodsC/coawst_4/use/fmrc/coawst_4_use_best.ncd']

mytime=dt.datetime.utcnow()

def time_near(cube,start):
    #    coord_names = [coord.name() for coord in cube.coords()]
    #    timevar = cube.coord(coord_names[0]))
    timevar=cube.coord('time')
    try:
        itime = timevar.nearest_neighbour_index(timevar.units.date2num(start))
    except:
        itime = -1
    return timevar.points[itime]


var='potential temperature' # long_name

for url in urls:
    cube = iris.load(url,iris.Constraint(name=var.strip()))[0]
    try:
        cube.coord(axis='T').rename('time')
    except:
        pass
    slice = cube.extract(iris.Constraint(time=time_near(cube,mytime)))

[Richard Hattersley]

> There has to be a better way, right?

There doesn't seem to be anything else that Iris should be doing automatically here. The two datasets just have different structures which can't be aligned without human intervention/interpretation.

That said, it does seem slightly restrictive that the only way to identify a coordinate in a constraint is by its name. As in this case it leads you to rename the ocean_time coordinate. I can imagine it might be useful to have something like:

def find_time_coord(cube):
    return [coord for coord in cube.coords() if coord.name().startswith('time')][0]

def time_near(time_coord, start):
    try:
        itime = time_coord.nearest_neighbour_index(time_coord.units.date2num(start))
    except:
        itime = -1
    return time_coord.points[itime]

for url in urls:
    cube = iris.load_cube(url, var)
    time_coord = find_time_coord(cube)
    iris.CoordConstraint(time_coord, time_near(time_coord, mytime))

But I wouldn't claim it's a dramatic improvement. And arguably what you're currently doing is better because by using the "rename()" method it's improving the quality of the metadata available for any subsequent operation on that cube.

NB. With your current rename('time') technique, you could still use "startswith" to identify the relevant coordinate if you so wish.

As an aside, I'd make a couple of other minor tweaks.
 - If you just want a single cube then I'd use iris.load_cube(), then you don't have to use a trailing [0] and it'll complain if you get multiple matches.
 - For a simple match on phenomenon name, you can replace your Constraint(name=...) with just the name. E.g. iris.load_cube(url, 'potential temperature')
 - It's a good habit to always avoid "bare" except clauses. It doesn't matter very much here, but under other circumstances it can cause problems, including making your Python code ignore your attempts to kill it with Ctrl-C.

[rsignell]

Richard,
Thanks for these tips.  Would it also be possible to have `cube.coord(axis='T')` return just the time coordinate variable?

Currently this bombs because it finds both the `time` coordinate and the `forecast_reference_time`, and is unhappy:

import iris
url='http://geoport.whoi.edu/thredds/dodsC/coawst_4/use/fmrc/coawst_4_use_best.ncd'

var='potential temperature' # long_name

cube = iris.load_cube(url,var)
cube.coord(axis='T')

CoordinateNotFoundError: u'Expected to find exactly 1 coordinate, but found 2. They were: time, forecast_reference_time.'

[Richard Hattersley]

On Wednesday, 4 December 2013 13:57:09 UTC, rsignell wrote:

Richard,
Thanks for these tips.  Would it also be possible to have `cube.coord(axis='T')` return just the time coordinate variable?

Up to now our interpretation of the CF conventions has been that any coordinate with units of "<something> since <reference time>" can be considered a time coordinate, as the conventions state, "A time coordinate is identifiable from its units string alone." And the reference to an optional standard name just says, "an appropriate value" rather than restricting it to "time".

Hence, "axis='T'" is selecting all coordinates that are "time-like", as defined above.

But this wouldn't be the first time that the conventions have been a touch ambiguous. ;-)

So I'm considering a post to the CF mailing list for clarification on intent and common usage. It seems just about possible the original intent was for all "time-like" quantities that are not "time" (i.e. the quantity with the standard name of "time") to require a standard name. And then any remaining coordinate with time-like units must be "time". (What happens if there are more than one of those I have no idea!)

[rsignell]

Richard,
Interesting.  The dataset in question is being automatically generated by the Forecast Model Run Collection capability of the THREDDS Data Server
<http://geoport.whoi.edu/thredds/dodsC/coawst_4/use/fmrc/coawst_4_use_best.ncd.html>
and hence many datasets will have these characteristics:
time:standard_name=time
time_run:standard_name=forecast_reference_time

Just to me sure I understand, are you proposing cube.coord(axis='T') could return only the variable with standard_name=time, and if that is not found, return the time coordinate that doesn't have a standard_name?

That would work for me!

-Rich

[rsignell]

I"m getting bit by this problem again.   All forecast model run collections from the Unidata THREDDS Data Server will fail here (we have 30 of these datasets here in the US Integrated Ocean Observing System).   Any further thoughts about resolution of this issue?

url='http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/his_Best/ESPRESSO_Real-Time_v2_History_Best_Available_best.ncd'
cube = iris.load_cube(url,'sea_water_potential_temperature')
timevar = cube.coord(axis='T')

results in:

---------------------------------------------------------------------------
CoordinateNotFoundError Traceback (most recent call last)
<ipython-input-53-fbe960bfde47> in <module>()
 1 url='http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/his_Best/ESPRESSO_Real-Time_v2_History_Best_Available_best.ncd'
 2 cube = iris.load_cube(url,'sea_water_potential_temperature')
----> 3 timevar = cube.coord(axis='T')

/home/local/python27_epd/lib/python2.7/site-packages/iris/cube.pyc in coord(self, name_or_coord, standard_name, long_name, var_name, attributes, axis, contains_dimension, dimensions, coord, coord_system, dim_coords, name)
 1252 'They were: %s.' % (len(coords), ', '.join(coord.name() for
 1253 coord in coords))
-> 1254 raise iris.exceptions.CoordinateNotFoundError(msg)
 1255 elif len(coords) == 0:
 1256 bad_name = name or standard_name or long_name or \

CoordinateNotFoundError: u'Expected to find exactly 1 coordinate, but found 2. They were: time, forecast_reference_time.'