Working with Axes¶

Constructing variables in memory¶

It is also often useful to create PyGeode variables directly in memory; for instance, to evaluate analytical expressions on a given geophysical grid. A very useful set of building blocks for creating variables in this way are axes objects. If you already have a variable on the relevant grid (defined from a file for instance), it’s easy to simply use the existing axes, as has been done in many cases in these tutorials. However, most axes objects have simple constructors to create them directly.

# The simplest, default case
In [1]: print(pyg.NamedAxis(np.arange(15), 'myaxis'))
myaxis <NamedAxis 'myaxis'>:  0  to 14  (15 values)

# Some simple examples
In [2]: lon = pyg.Lon(180)

In [3]: print(lon)
lon <Lon>      :  0 E to 179 E (180 values)

In [4]: lat = pyg.Lat(92)

In [5]: print(lat)
lat <Lat>      :  EQ to 91 N (92 values)

In [6]: pres = pyg.Pres([1000, 900, 800, 700, 500, 300, 200, 100, 50, 30, 10])

# Gaussian latitudes (with appropriate weights)
In [7]: lat2 = pyg.gausslat(64)

In [8]: print(lat2)
lat <Lat>      :  88 S to 88 N (64 values)

Time axes are somewhat more complicated, as you need to specify the calendar, the reference date (startdate), offsets, and the native unit:

In [9]: time = pyg.ModelTime365(values=np.arange(3650), units='days', startdate=dict(year=2000, month=1))

but some simple helpers are available here as well:

# A time axis of fixed length
In [10]: print(pyg.modeltime365n('1 Jan 2000', 3650, units='days'))
time <ModelTime365>:  Jan 1, 2000 00:00:00 to Dec 31, 2009 00:00:00 (3650 values)

# A time axis spanning a range of dates
In [11]: print(pyg.modeltime365range('1 Jan 2000', '1 Jan 2001', step=6, units='hours'))
time <ModelTime365>:  Jan 1, 2000 00:00:00 to Dec 31, 2000 18:00:00 (1460 values)

Usually the easiest approach to creating variables in memory is to apply the relevant mathematical operations to the axes themselves:

In [12]: pyg.sin(2*np.pi*time / 365) * pyg.exp(-(lat2 / 20.)**2)
Out[12]: <Var '(sin(time)*exp(-lat**2))'>

In [13]: pyg.sin(2*np.pi*time / 365) * pyg.exp(-(lat2 / 20.)**2)
Out[13]: <Var '(sin(time)*exp(-lat**2))'>

However, if you have a numpy array that you want to turn in to a PyGeode variable, this can also be done.

In [14]: x = np.ones((64, 180), 'd')

In [15]: print(pyg.Var((lat2, lon), name = 'myvar', values=x))
<Var 'myvar'>:
  Shape:  (lat,lon)  (64,180)
  Axes:
    lat <Lat>      :  88 S to 88 N (64 values)
    lon <Lon>      :  0 E to 179 E (180 values)
  Attributes:
    {}
  Type:  Var (dtype="float64")

Time Axis Operations¶

PyGeode recognizes many types of geophysical dimensions, some of which we’ve already seen in these tutorials. By far the most complicated of these are those axes representing the time dimension. PyGeode has a reasonably sophisticated time axis, which is aware of four main types of Calandars - the standard Julian calendar (StandardTime), a 365-day calendar (ModelTime365), a 360-day calendar (ModelTime360), and a yearless calendar (Yearless) with no temporal organization at scales longer than days.

Creating Time Axes¶

Time axes will often be created automatically from datasets with CF-compliant metadata (of which PyGeode is aware); but they can also be created manually. There are two broad ways to create a time axis: 1) in terms of an offset from a reference date and 2) by directly specifying dates and times.

The first requires specifying a reference date (startdate), a list of values representing offsets from this date, and the units in which these offsets are specified (seconds, minutes, hours, or days); this is usually the more convenient way to construct an axis.

In [16]: import pygeode as pyg, numpy as np

# A time axis representing the (leap) year 2008 at 6 h increments
In [17]: print(pyg.StandardTime(values=np.arange(0, 24*366, 6), units = 'hours', startdate=dict(year=2008, month=1, day=1)))
time <StandardTime>:  Jan 1, 2008 00:00:00 to Dec 31, 2008 18:00:00 (1464 values)

# 2000 to 2010 in 360-day years
In [18]: print(pyg.ModelTime360(values=np.arange(10*360), units = 'days', startdate=dict(year=2000, month=1, day=1)))
time <ModelTime360>:  Jan 1, 2000 00:00:00 to Dec 30, 2009 00:00:00 (3600 values)

The second involves providing an explicit list of the dates and times of each element of the axis. PyGeode recognizes the auxiliary arrays year, month, day, hour, minute and second, though not all need to be explicitly provided:

# A time axis representing the months of 1890
In [19]: print(pyg.ModelTime365(year=12*[1890], month=range(1, 13), units='days'))
time <ModelTime365>:  Jan , 1890 :: to Dec , 1890 :: (12 values)

Combining Time Axes¶

Internally the axis values (used for mapping two axes, selecting subranges, and plotting) are represented as offsets from the reference startdate. PyGeode will map between the two if a valid map exists (Broadcasting), though note that which reference frame the result will be returned in can depend on the order of operations:

# A time axis representing the year 2005 at 6 h increments
In [20]: t1 = pyg.modeltime365range('1 Jan 2005', '1 Jan 2006', step=6, units = 'hours')

# Another time axis representing the same dates but with a different reference date and different units
In [21]: t2 = pyg.modeltime365range('1 Jan 2005', '1 Jan 2006', step=0.25, units = 'days', ref = '1 Jan 2004')

# Add these in two different ways
In [22]: ts12 = t1 + t2

# Since either axis can be mapped to each other, these will differ
In [23]: ts21 = t2 + t1

# Both of these work and refer to the same time period...
In [24]: print(ts12.time)
time <ModelTime365>:  Jan 1, 2005 00:00:00 to Dec 31, 2005 18:00:00 (1460 values)

In [25]: print(ts21.time)
time <ModelTime365>:  Jan 1, 2005 00:00:00 to Dec 31, 2005 18:00:00 (1460 values)

# ...but are defined in different reference frames, with different units.
In [26]: print(ts12.time.units, ts21.time.units)
hours days

# Accessing absolute times will refer to the same value in each case
In [27]: print(ts12(time='15 Jan 2005').time, ts21(time='15 Jan 2005').time)
time <ModelTime365>:  Jan 15, 2005 00:00:00 time <ModelTime365>:  Jan 15, 2005 00:00:00

# But relative times will differ
In [28]: print(ts12(time=400).time, ts21(time=400).time)
time <ModelTime365>:  Jan 17, 2005 18:00:00 time <ModelTime365>:  Feb 5, 2005 00:00:00

Accessing any results by an absolute date will therefore yield the same result, but accessing by value will not. Moreover since PyGeode used the values array for plotting, this can affect details of plots. At present this ambiguity only exists if the mapping can be done both ways; if one is a strict subset of the other the result will always be in the reference frame of the smaller axis.

Climatological statistics¶

A variety of time averaging operations are possible. In addition to computing climatologies, one can compute monthly means, daily means, diurnal means, and trends. These should generally work as expected; however, PyGeode also implements special mapping rules for time axes so that anomalies from these time averages can be easily defined. In general these are based on the

oOne common operation is to compute climatologies. As described in the first of these tutorials this is simply performed as follows

In [29]: from pygeode.tutorial import t2

In [30]: Tc = pyg.climatology(t2.Temp)    # Compute the climatology

Time axis utilities¶

There are a number of other useful utility functions for working with time axes in the timeutils module. For instance, it is sometimes convenient to convert data on a a standardtime axis

lagged variables