Calculus along a band axis, and calendars¶

Seven members, in two groups. What unites the first four is that they measure a variable against a band dimension's own coordinates, not against its step count — which is exactly what separates them from the members that only count steps. (cumulative is the odd one out: it is an accessor over cumsum / cumprod, which read no coordinates at all.)

Member	What it does
differentiate	the derivative along a band axis, by its coordinate spacing
integrate	the trapezoidal integral along a band axis, which it consumes
cumulative_integrate	the running integral, keeping the axis
polyfit	per-cell least-squares coefficients along a band axis
cumulative	the accessor form of the running total / product
convert_calendar	restamp a time axis onto another CF calendar
interp_calendar	interpolate the values onto another calendar's stamps

On an evenly spaced axis a wrong implementation that ignores the coordinates still looks right, so the examples that establish a member's contract use an uneven one. A few later examples keep an even axis deliberately, where the point is the shape of the answer rather than the spacing.

Setup¶

One notebook-relative data path. The sample files ship in the repo's examples/data/netcdf/ folder.

In [1]:

from pathlib import Path

import numpy as np

from pyramids.netcdf import ExtraDimensions, GeoReference, NetCDF

SAMPLES = Path('../../../examples/data/netcdf/samples')

# A COARDS cube: one variable, `rhum`, on a 12-step time axis and a 4-step pressure axis.
LEVELS = SAMPLES / 'coards__5v__1d4-4d1__y-desc.nc'

from pathlib import Path import numpy as np from pyramids.netcdf import ExtraDimensions, GeoReference, NetCDF SAMPLES = Path('../../../examples/data/netcdf/samples') # A COARDS cube: one variable, `rhum`, on a 12-step time axis and a 4-step pressure axis. LEVELS = SAMPLES / 'coards__5v__1d4-4d1__y-desc.nc'

The cube, and why the spacing matters¶

rhum sits on (time, level, lat, lon). Look at the level axis: the gaps are not equal, and the last one is twice the others. Any member that treats this axis as unit steps is simply wrong about it.

In [2]:

nc = NetCDF.read_file(str(LEVELS))
rhum = nc.get_variable('rhum')

levels = np.asarray(rhum.coords['level'])
rhum.dimension_names, levels.tolist(), np.diff(levels).tolist()

nc = NetCDF.read_file(str(LEVELS)) rhum = nc.get_variable('rhum') levels = np.asarray(rhum.coords['level']) rhum.dimension_names, levels.tolist(), np.diff(levels).tolist()

Out[2]:

(['time', 'level', 'lat', 'lon'],
 [1000.0, 925.0, 850.0, 700.0],
 [-75.0, -75.0, -150.0])

differentiate — a rate, not a step¶

diff subtracts neighbouring steps and shortens the axis. differentiate divides by the real distance between the coordinates, so it answers a rate and keeps the axis.

A straight line makes the difference obvious. Values 0, 10, 30 at positions 0, 1, 3 lie on one line of slope 10, so the derivative is 10 everywhere — while the raw differences are 10 and 20.

In [3]:

geo = GeoReference(geo=(0.0, 1.0, 0.0, 1.0, 0.0, -1.0), epsg=4326)
line = NetCDF.from_array(
    np.array([0.0, 10.0, 30.0]).reshape(3, 1, 1),
    geo_ref=geo,
    variable_name='v',
    dims=ExtraDimensions(name='level', values=[0.0, 1.0, 3.0]),
).get_variable('v')

(
    line.diff('level').read_array().ravel().tolist(),
    line.differentiate('level').read_array().ravel().tolist(),
)

geo = GeoReference(geo=(0.0, 1.0, 0.0, 1.0, 0.0, -1.0), epsg=4326) line = NetCDF.from_array( np.array([0.0, 10.0, 30.0]).reshape(3, 1, 1), geo_ref=geo, variable_name='v', dims=ExtraDimensions(name='level', values=[0.0, 1.0, 3.0]), ).get_variable('v') ( line.diff('level').read_array().ravel().tolist(), line.differentiate('level').read_array().ravel().tolist(), )

Out[3]:

([10.0, 20.0], [10.0, 10.0, 10.0])

On the real cube it keeps every dimension, so the result is still a cube you can sel into:

In [4]:

rate = rhum.differentiate('level')

rate.dimension_names, rate.shape == rhum.shape

rate = rhum.differentiate('level') rate.dimension_names, rate.shape == rhum.shape

Out[4]:

(['time', 'level', 'lat', 'lon'], True)

One thing worth knowing about gaps: a gap spoils the steps that read it, not its own. A central difference at step i reads i-1 and i+1 and never i, so the step holding the gap is the one step the gap does not spoil.

integrate — an area, and the axis is consumed¶

reduce(dim, how="sum") adds the step values and ignores the spacing. integrate weights each trapezoid by its own gap, so widen the axis and the integral grows while the sum does not move.

In [5]:

flat = NetCDF.from_array(
    np.full((3, 1, 1), 2.0),
    geo_ref=geo,
    variable_name='v',
    dims=ExtraDimensions(name='level', values=[0.0, 1.0, 10.0]),
).get_variable('v')

(
    float(flat.reduce('level', 'sum').read_array()[0, 0]),
    float(flat.integrate('level').read_array()[0, 0]),
)

flat = NetCDF.from_array( np.full((3, 1, 1), 2.0), geo_ref=geo, variable_name='v', dims=ExtraDimensions(name='level', values=[0.0, 1.0, 10.0]), ).get_variable('v') ( float(flat.reduce('level', 'sum').read_array()[0, 0]), float(flat.integrate('level').read_array()[0, 0]), )

Out[5]:

(6.0, 20.0)

The integrated axis is consumed, exactly as a full reduce collapses one — a (time, level) cube integrated over level comes back on time alone:

In [6]:

column = rhum.integrate('level')

rhum._band_dim_names, column._band_dim_names

column = rhum.integrate('level') rhum._band_dim_names, column._band_dim_names

Out[6]:

(('time', 'level'), ('time',))

The sign follows the coordinates rather than being normalised away. A pressure axis descends, so integrating down it is negative — that is what the coordinates say, and it is not second-guessed. Reverse the axis, or negate the result, when the magnitude is what you want.

In [7]:

round(float(column.read_array().ravel()[0]), 3)

round(float(column.read_array().ravel()[0]), 3)

Out[7]:

-14653.128

cumulative_integrate — the running integral¶

integrate's length-preserving twin: the integral up to each step. Its last step is what integrate answers for the whole axis.

It starts at zero, because no interval has been traversed yet — a real asymmetry with cumsum, whose first step holds the first value.

In [8]:

steady = NetCDF.from_array(
    np.full((4, 1, 1), 2.0),
    geo_ref=geo,
    variable_name='v',
    dims=ExtraDimensions(name='level', values=[0.0, 1.0, 2.0, 3.0]),
).get_variable('v')

(
    steady.cumsum('level').read_array().ravel().tolist(),
    steady.cumulative_integrate('level').read_array().ravel().tolist(),
    float(steady.integrate('level').read_array()[0, 0]),
)

steady = NetCDF.from_array( np.full((4, 1, 1), 2.0), geo_ref=geo, variable_name='v', dims=ExtraDimensions(name='level', values=[0.0, 1.0, 2.0, 3.0]), ).get_variable('v') ( steady.cumsum('level').read_array().ravel().tolist(), steady.cumulative_integrate('level').read_array().ravel().tolist(), float(steady.integrate('level').read_array()[0, 0]), )

Out[8]:

([2.0, 4.0, 6.0, 8.0], [0.0, 2.0, 4.0, 6.0], 6.0)

polyfit — a per-cell trend¶

The least-squares fit, vectorised over every cell. This is the first member that replaces a dimension with a differently-named one: fitting degree n over k steps answers n+1 coefficients, so level becomes degree of length deg + 1.

In [9]:

trend = NetCDF.from_array(
    np.array([1.0, 3.0, 5.0, 7.0]).reshape(4, 1, 1),
    geo_ref=geo,
    variable_name='v',
    dims=ExtraDimensions(name='time', values=[0.0, 1.0, 2.0, 3.0]),
).get_variable('v')

fit = trend.polyfit('time', 1)
(
    fit._band_dim_names,
    fit._band_dim_sizes,
    [round(c, 6) for c in fit.read_array().ravel().tolist()],
)

trend = NetCDF.from_array( np.array([1.0, 3.0, 5.0, 7.0]).reshape(4, 1, 1), geo_ref=geo, variable_name='v', dims=ExtraDimensions(name='time', values=[0.0, 1.0, 2.0, 3.0]), ).get_variable('v') fit = trend.polyfit('time', 1) ( fit._band_dim_names, fit._band_dim_sizes, [round(c, 6) for c in fit.read_array().ravel().tolist()], )

Out[9]:

(('degree',), (2,), [2.0, 1.0])

[2.0, 1.0] is 2x + 1 — highest power first, which is numpy.polyfit's own order, and the same order xarray's polyfit stamps its degree axis with. The two agree, so a result moves between them without reindexing, and the degree stamps say which power each coefficient belongs to:

In [10]:

np.asarray(trend.polyfit('time', 2).coords['degree']).tolist()

np.asarray(trend.polyfit('time', 2).coords['degree']).tolist()

Out[10]:

[2, 1, 0]

Two things to know. The fit uses the coordinates as sample positions, so an uneven axis fits correctly. And numpy.polyfit has no gap concept, so a cell whose series holds a gap answers all-NaN coefficients rather than a fit over the steps that remain — run dropna or interpolate_na first when a partial series should still fit.

cumulative — the accessor spelling¶

xarray spells the running total cumulative(dim).sum(). pyramids already answered it as cumsum, so this is the accessor over the shipped members: it forwards, and a test pins the equivalence so the two spellings cannot drift.

In [11]:

(
    steady.cumulative('level').sum().read_array().ravel().tolist()
    == steady.cumsum('level').read_array().ravel().tolist(),
    repr(steady.cumulative('level')),
)

( steady.cumulative('level').sum().read_array().ravel().tolist() == steady.cumsum('level').read_array().ravel().tolist(), repr(steady.cumulative('level')), )

Out[11]:

(True, "CumulativeAccessor(dim='level')")

Calendars¶

Climate-model output rarely uses a real calendar: 360_day, noleap and all_leap are bookkeeping calendars, and pyramids already reads them — cftime is a core dependency and the store's calendar attribute is parsed on the way in. These two members are the other half.

They divide the work, and the division is the whole point:

• convert_calendar keeps the values and moves the stamps. A date the target calendar has no day for is dropped, and its step with it.
• interp_calendar keeps the steps and moves the values, interpolating onto the target's own stamps. Nothing is lost.

In [12]:

def model_cube(offsets, calendar, values):
    """A small cube on a model calendar."""
    planes = np.stack([np.full((1, 1), float(v)) for v in values])
    return NetCDF.from_array(
        planes,
        geo_ref=geo,
        variable_name='v',
        dims=ExtraDimensions(
            name='time',
            values=[float(o) for o in offsets],
            attrs={'time': {'units': 'days since 2001-01-01', 'calendar': calendar}},
        ),
    ).get_variable('v')


on_360 = model_cube([0.0, 60.0], '360_day', [1.0, 2.0])
on_360._resolved_band_dim_time_attrs()['time']

def model_cube(offsets, calendar, values): """A small cube on a model calendar.""" planes = np.stack([np.full((1, 1), float(v)) for v in values]) return NetCDF.from_array( planes, geo_ref=geo, variable_name='v', dims=ExtraDimensions( name='time', values=[float(o) for o in offsets], attrs={'time': {'units': 'days since 2001-01-01', 'calendar': calendar}}, ), ).get_variable('v') on_360 = model_cube([0.0, 60.0], '360_day', [1.0, 2.0]) on_360._resolved_band_dim_time_attrs()['time']

Out[12]:

('days since 2001-01-01', '360_day')

convert_calendar restamps. Day 60 of a 360-day year is 1 March (three 30-day months); on a 365-day calendar 1 March is day 59, so the offset moves by one while the values stay put:

In [13]:

moved = on_360.convert_calendar('noleap')

(
    moved._resolved_band_dim_time_attrs()['time'][1],
    np.asarray(moved.coords['time']).tolist(),
    moved.read_array().ravel().tolist(),
)

moved = on_360.convert_calendar('noleap') ( moved._resolved_band_dim_time_attrs()['time'][1], np.asarray(moved.coords['time']).tolist(), moved.read_array().ravel().tolist(), )

Out[13]:

('noleap', [0.0, 59.0], [1.0, 2.0])

When a date has no counterpart the step is dropped. Under all_leap, day 59 is 29 February — a day noleap does not have — so converting loses that step and the axis gets shorter.

In [14]:

leap = model_cube([58.0, 59.0, 60.0], 'all_leap', [1.0, 2.0, 3.0])

import warnings

with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    dropped = leap.convert_calendar('noleap')

leap._band_dim_sizes, dropped._band_dim_sizes, dropped.read_array().ravel().tolist()

leap = model_cube([58.0, 59.0, 60.0], 'all_leap', [1.0, 2.0, 3.0]) import warnings with warnings.catch_warnings(): warnings.simplefilter('ignore') dropped = leap.convert_calendar('noleap') leap._band_dim_sizes, dropped._band_dim_sizes, dropped.read_array().ravel().tolist()

Out[14]:

((3,), (2,), [1.0, 3.0])

align_on='year' drops nothing instead, by keeping the position in the year rather than the calendar date — which is what you want when the dates are model bookkeeping rather than real days:

In [15]:

kept = leap.convert_calendar('noleap', align_on='year')

kept._band_dim_sizes, kept.read_array().ravel().tolist()

kept = leap.convert_calendar('noleap', align_on='year') kept._band_dim_sizes, kept.read_array().ravel().tolist()

Out[15]:

((3,), [1.0, 2.0, 3.0])

interp_calendar is the lossless route: it interpolates the values onto another cube's stamps over a calendar-independent decimal-year scale, so every target step is answered and the two cubes line up for concat or an operator.

In [16]:

source = model_cube([0.0, 180.0, 359.0], '360_day', [0.0, 10.0, 20.0])
onto = model_cube([90.0, 270.0], 'noleap', [0.0, 0.0])

aligned = source.interp_calendar(onto)
(
    aligned._band_dim_sizes,
    aligned._resolved_band_dim_time_attrs()['time'][1],
    [round(v, 3) for v in aligned.read_array().ravel().tolist()],
)

source = model_cube([0.0, 180.0, 359.0], '360_day', [0.0, 10.0, 20.0]) onto = model_cube([90.0, 270.0], 'noleap', [0.0, 0.0]) aligned = source.interp_calendar(onto) ( aligned._band_dim_sizes, aligned._resolved_band_dim_time_attrs()['time'][1], [round(v, 3) for v in aligned.read_array().ravel().tolist()], )

Out[16]:

((2,), 'noleap', [4.932, 14.821])

Recap¶

• differentiate(dim) is a rate — diff divided by the real coordinate spacing, axis kept. A gap spoils its neighbours, not its own step.
• integrate(dim) is an area — the trapezoid rule, consuming the axis, signed as the coordinates say.
• cumulative_integrate(dim) is the running integral, axis kept, first step zero.
• polyfit(dim, deg) fits per cell and replaces the axis with degree, highest power first.
• cumulative(dim).sum() / .prod() is the accessor over the shipped cumsum / cumprod.
• convert_calendar(cal) keeps the values and may lose a step; interp_calendar(other) keeps every step and moves the values. Reach for the second when losing a step is unacceptable.

All five numerical members read the dimension's coordinates, so they refuse an axis that carries none, a spatial axis, and an axis too short for the operation — rather than quietly treating step positions as measurements.