Thanks much for the quick response. I updated both matplotlib andbasemap (now at 0.99.5) via svn and noticed the new netcdftime.py. First, from within site-packages/mpl_toolkits/basemap,
$ grep date2index *.py __init__.py::func:`date2index`: compute a time variable index corresponding to a date. __init__.py:def date2index(dates, nctime, calendar='proleptic_gregorian'): __init__.py: return netcdftime.date2index(dates, nctime, calendar=None) netcdftime.py:def date2index(dates, nctime, calendar=None, select='exact'): netcdftime.py: date2index(dates, nctime, calendar=None, select='exact') so there seems to be some disagreement between __init__.py and netcdftime.py concerning the presence of the "select" argument. When I call date2index with the "select" keyword arg I get In [24]: ix0 = date2index(date0,timedata,timedata.calendar,select='nearest') --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /home/amg/work/nhmm/<ipython console> in <module>() TypeError: date2index() got an unexpected keyword argument 'select' ----------------------- This detail aside, I am still having difficulty with date2index, butannoyingly, I seem to get different error messages with different datasets. I'll illustrate two here, starting with the one I initially posted about. (See note below regarding this data.)
In [3]: from mpl_toolkits.basemap import date2index,num2date,NetCDFFile as ncf
In [10]: from mpl_toolkits import basemap In [11]: print basemap.__version__ 0.99.5 In [24]: fname0 = 'http://esgcet.llnl.gov/dap/'In [25]: fname1 = 'ipcc4/20c3m/gfdl_cm2_1/pcmdi.ipcc4.gfdl_cm2_1.20c3m.run1.atm.mo.xml'
In [26]: fname = fname0+fname1 In [28]: datobj = ncf(fname) In [33]: datobj.variables['tas'].shape Out[33]: (1680, 90, 144) In [34]: timedata = datobj.variables['time'] In [35]: from datetime import datetime as dt In [36]: date0 = dt(1951,1,16,12,0,0) In [37]: print date0 1951-01-16 12:00:00 In [38]: nt0 = date2index(date0,timedata) --------------------------------------------------------------------------- ClientError Traceback (most recent call last) /home/amg/work/nhmm/<ipython console> in <module>()/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/__init__.pyc in date2index(dates, nctime, calendar)
3931 Returns an index or a sequence of indices. 3932 """ -> 3933 return netcdftime.date2index(dates, nctime, calendar=None) 3934 3935 def maskoceans(lonsin,latsin,datain,inlands=False):/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/netcdftime.pyc in date2index(dates, nctime, calendar, select)
1006 # If the times do not correspond, then it means that the times 1007 # are not increasing uniformly and we try the bisection method. -> 1008 if not _check_index(index, dates, nctime, calendar): 1009 1010 # Use the bisection method. Assumes the dates are ordered./home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/netcdftime.pyc in _check_index(indices, dates, nctime, calendar)
959 return False
960
--> 961 t = nctime[indices]
962 return numpy.all( num2date(t, nctime.units, calendar) == dates)
963
/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/netcdf.pyc
in __getitem__(self, index)
65
66 def __getitem__(self, index):
---> 67 datout = squeeze(self._var.__getitem__(index))
68 # automatically
69 # - remove singleton dimensions
/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/dap/dtypes.pyc
in __getitem__(self, key)
409 def __getitem__(self, key):
410 # Return data from the array.
--> 411 return self.data[key]
412
413 def __setitem__(self, key, item):
/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/dap/proxy.pyc
in __getitem__(self, index)
112
113 # Fetch data.
--> 114 resp, data = openurl(url, self.cache, self.username,
self.password)
115
116 # First lines are ASCII information that end with
'Data:\n'.
/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/dap/util/http.pyc in openurl(url, cache, username, password) 19 m = re.search('code = (?P<code>\d+);\s*message = "(?P<msg>.*)"', data, re.DOTALL | re.MULTILINE)
20 msg = 'Server error %(code)s: "%(msg)s"' % m.groupdict()
---> 21 raise ClientError(msg)
22
23 return resp, data
ClientError: 'Server error 0: "invalid literal for int(): [1113"'
-----------------------------
Note that this is a different error than previously reported. Also, the
correct time index is still 1080:
In [40]: taxvals = datobj.variables['time'][:] In [41]: num2date(taxvals[1080],timedata.units,timedata.calendar) Out[41]: 1951-01-16 12:00:00 -----------------------------This dataset, generated by one of the IPCC models, is password-protected, but could be a good target for decoding, since it is typical of a large class of climate models, that generate a lot of analytical activity. To get a password (they're free) one must register. Info is here: http://www-pcmdi.llnl.gov/ipcc/info_for_analysts.php. Follow "How to access..." then "Register to download output." Once you get a userid and password they can be inserted in the NetCDFFile call, voila. Note that there is a new iteration of IPCC coming down the pike; new model files to become widely available starting in 2010.
------------------------------The underlying data is available via ftp. I fetched it and extracted a small slab, which is available at http://iri.columbia.edu/~amg/test/gfdl_test.nc. The CDAT package can digest this file; first time step is plotted here: http://iri.columbia.edu/~amg/test/gfdl_test_time0.png. The dates can also be read by this package, and run from Aug 1950 to Aug 1951, inclusive (13 mos). So the file does not seem to be garbage.
In [16]: datobj = ncf('gfdl_test.nc')
In [17]: timedata = datobj.variables['time']
In [18]: date0 = dt(1951,1,16,12,0,0)
In [19]: nt0 = date2index(date0,timedata)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
/home/amg/work/nhmm/<ipython console> in <module>()
/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/__init__.pyc
in date2index(dates, nctime, calendar)
3931 Returns an index or a sequence of indices. 3932 """ -> 3933 return netcdftime.date2index(dates, nctime, calendar=None) 3934 3935 def maskoceans(lonsin,latsin,datain,inlands=False):/home/amg/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/netcdftime.pyc in date2index(dates, nctime, calendar, select)
1011 import bisect 1012-> 1013 index = numpy.array([bisect.bisect_left(nctime, n) for n in num], int)
1014 1015 nomatch = num2date(nctime[index], nctime.units) != dates TypeError: object of type 'netcdf_variable' has no len() Investigating the time axis, In [20]: taxvals = timedata[:] In [21]: taxvalsOut[21]: array([ 32712.5, 32743. , 32773.5, 32804. , 32834.5, 32865.5, 32895. , 32924.5, 32955. , 32985.5, 33016. , 33046.5, 33077.5])
In [22]: num2date(taxvals,timedata.units,timedata.calendar)Out[22]: array([1950-08-16 12:00:00, 1950-09-16 00:00:00, 1950-10-16 12:00:00,
1950-11-16 00:00:00, 1950-12-16 12:00:00, 1951-01-16 12:00:00,
1951-02-15 00:00:00, 1951-03-16 12:00:00, 1951-04-16 00:00:00,
1951-05-16 12:00:00, 1951-06-16 00:00:00, 1951-07-16 12:00:00,
1951-08-16 12:00:00], dtype=object)
Which agrees with what CDAT sees.
-------------------------
I think this is enough for now. I also had problems opening data files
whose time units were like "months since xxxx-xx-xx," since the "months"
unit does not seem to be supported. ("years since..." could also be
useful in some cases.) But maybe one or two things at a time is enough!
Thanks for any assistance/advice! Best, Arthur Jeff Whitaker wrote:
David Huard wrote: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 HuardArthur: I've just updated basemap svn with David's latest version of date2index, so another option is to update basemap from svn. Or, even simpler, just drop the attached netcdftime.py file in lib/mpl_toolkits/basemap (replacing the old one) and run python setup.py install.-JeffOn Mon, Sep 7, 2009 at 9:27 AM, Arthur M. Greene <a...@iri.columbia.edu <mailto:a...@iri.columbia.edu>> 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/<ipython console> in <module>()/usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/__init__.pycin 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.pycin date2index(dates, nctime, calendar) 986 987 # Perform check again. --> 988 _check_index(index, dates, nctime, calendar) 989990 # convert numpy scalars or single element arrays to pythonints./usr/local/cdat/trunk/lib/python2.5/site-packages/mpl_toolkits/basemap/netcdftime.pycin _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 thecorrect 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 <mailto:Matplotlib-users@lists.sourceforge.net> https://lists.sourceforge.net/lists/listinfo/matplotlib-users ------------------------------------------------------------------------------------------------------------------------------------------------------ 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
-- *^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~* Arthur M. Greene, Ph.D. The International Research Institute for Climate and Society (IRI) The Earth Institute, Columbia University, Lamont Campus Monell Building, 61 Route 9W, Palisades, NY 10964-8000 USA Tel: 845-680-4436 | Fax: 845-680-4865 a...@iri.columbia.edu | http://iri.columbia.edu *^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*^*~*
<<attachment: amg.vcf>>
------------------------------------------------------------------------------ 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
