Fetching a data slab with NetCDFFile

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
*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*

Arthur,

I wrote the date2index function and I think what you are seeing is a bug that I fixed a couple of months ago. By using the latest version of netcdf4-python, not only should this bug disappear, but you’ll also find that date2index now supports different selection methods: ‘exact’, ‘before’, ‘after’, ‘nearest’, that should help with your use case.

If this does not fix the problem you are seeing, I’d appreciate having a copy of the file and code to reproduce the problem and find a solution.

HTH,

David Huard

···

On Mon, Sep 7, 2009 at 9:27 AM, Arthur M. Greene <amg@…1985…> wrote:

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/ in ()

/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

^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~*


Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day

trial. Simplify your report design, integration and deployment - and focus on

what you do best, core application coding. Discover what’s new with

Crystal Reports now. http://p.sf.net/sfu/bobj-july


Matplotlib-users mailing list

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users