[Matplotlib-users] Fetching a data slab with NetCDFFile

[Arthur M. G. <am...@ir...>]

Hi All,

The problem is not with fetching the data slice itself, but finding the 
correct indices to specify, particularly with the time dimension. The 
below examples refer to a remote dataset that I can open and slice using 
indices, as in

slice = remoteobj.variables['tas'][:120,20:40,30:50].

However, I have problems when trying to use the syntax in plotsst.py or 
pnganim.py (from the examples) to find time indices:

In [107]: from datetime import datetime as dt
In [108]: date0 = dt(1951,1,1,0)
In [110]: print date0
1951-01-01 00:00:00

In [125]: timedata = remoteobj.variables['time']
In [126]: nt0 = date2index(date0,timedata)
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)

/home/amg/work/nhmm/<ipython console> in <module>()

/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/__init__.pyc 

in date2index(dates, nctime, calendar)
    3924     Returns an index or a sequence of indices.
    3925     """
-> 3926     return netcdftime.date2index(dates, nctime, calendar=None)
    3927
    3928 def maskoceans(lonsin,latsin,datain,inlands=False):

/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/netcdftime.pyc 

in date2index(dates, nctime, calendar)
     986
     987         # Perform check again.
--> 988         _check_index(index, dates, nctime, calendar)
     989
     990     # convert numpy scalars or single element arrays to python
ints.

/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/netcdftime.pyc 

in _check_index(indices, dates, nctime, calendar)
     941     for n,i in enumerate(indices):
     942         t[n] = nctime[i]
--> 943     assert numpy.all( num2date(t, nctime.units, calendar) == dates)
     944
     945

AssertionError:

---------------------------------------------------------

It turns out that date0 corresponds best to index 1080:

In [139]: remoteobj.variables['time'][1080]
Out[139]: 32865.5

In [141]: num2date(32865.5,timedata.units,timedata.calendar)
Out[141]: 1951-01-16 12:00:00

This isn't the _exact_ date and time I had specified, but

In [142]: date0 = dt(1951,01,16,12,00,00)
In [143]: print date0
1951-01-16 12:00:00

In [144]: date2index(date0,timedata,timedata.calendar)

produces the same AssertionError. Where is the problem?

What I would _like_ to do is to issue a simple call using coordinates 
rather than the indices, of the form:

slice = variable[date0:date1,[plev],lat0:lat1,lon0:lon1],

or similar, preferably without writing a whole module just to find the 
correct indices. I need to fetch similar slices from a group of models, 
having time axes that may each be defined slightly differently -- 
different calendars, time point set at a different day of the month, 
etc. (It's monthly data and I'm specifying only monthly bounds, even 
though the calendar may be defined as "days since 1860...") I need to 
automate the process so I get back the correct slab regardless.

Suggestions appreciated!

Thx,

Arthur


*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*
Arthur M. Greene, Ph.D.
The International Research Institute for Climate and Society (IRI)
The Earth Institute, Columbia University, Lamont Campus

amg at iri dot columbia dot edu | http://iri.columbia.edu
*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*