Fitting a model, smoothing a series, and the cube as a dict¶

Four members, in two groups.

curvefit and rolling_exp both run along a band axis. curvefit is polyfit's general case — an arbitrary model rather than a polynomial — and rolling_exp is rolling's, with an infinite window instead of a fixed one.

to_dict and from_dict go the other way: out of the cube entirely, to plain dict, list, str and numbers, and back again.

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'

cube = NetCDF.read_file(str(LEVELS))
rhum = cube.get_variable('rhum')
rhum._band_dim_names, rhum._band_dim_sizes

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' cube = NetCDF.read_file(str(LEVELS)) rhum = cube.get_variable('rhum') rhum._band_dim_names, rhum._band_dim_sizes

Out[1]:

(('time', 'level'), (12, 4))

curvefit — an arbitrary model, per cell¶

polyfit fits a polynomial, which is linear in its coefficients, so a whole cube goes through one vectorised numpy.polyfit call. An arbitrary model is not linear in its parameters, so each cell runs its own scipy.optimize.curve_fit: about 195 µs per cell, with no vectorised form. xarray's curvefit goes per-cell for the same reason and measures slightly slower, so that cost is the method's rather than this route's.

Fit relative humidity against pressure with a straight line. p0 is positional — one entry per parameter — and it is required, because its length is what tells the member how many parameters the model has.

In [2]:

def line(x, slope, intercept):
    return slope * x + intercept


fit = rhum.curvefit('level', line, [0.0, 50.0])
fit._band_dim_names, fit._band_dim_sizes

def line(x, slope, intercept): return slope * x + intercept fit = rhum.curvefit('level', line, [0.0, 50.0]) fit._band_dim_names, fit._band_dim_sizes

Out[2]:

(('time', 'param'), (12, 2))

The fitted axis is replaced. level had four steps; param has one slot per parameter, stamped positionally 0 and 1. The stamps are integers rather than the parameter names slope and intercept because a text axis would be inert to the rest of the library — differentiate, integrate and interp all refuse a non-numeric one.

In [3]:

np.asarray(fit.coords['param']).tolist()

np.asarray(fit.coords['param']).tolist()

Out[3]:

[0, 1]

In [4]:

# `time` is untouched, so the result is (time, param, y, x): one fit per time step per cell.
coefficients = np.asarray(fit.read_array()).reshape(
    *fit._band_dim_sizes, *rhum.shape[-2:]
)
slope, intercept = coefficients[0, 0], coefficients[0, 1]
float(slope[18, 36]), float(intercept[18, 36])

# `time` is untouched, so the result is (time, param, y, x): one fit per time step per cell. coefficients = np.asarray(fit.read_array()).reshape( *fit._band_dim_sizes, *rhum.shape[-2:] ) slope, intercept = coefficients[0, 0], coefficients[0, 1] float(slope[18, 36]), float(intercept[18, 36])

Out[4]:

(0.08285714100514147, -1.1071324053087963)

full=True, and why it is worth asking for¶

curve_fit does not report a failed fit. On a degenerate series it returns p0 unchanged, does not raise, and reports ier=1 — a success code. Its covariance is no help and is in fact inverted: a degenerate cell comes back finite where a perfect fit comes back inf.

So an unfittable cell arrives as a smooth region of exactly p0, indistinguishable from a real result. On a map that is worse than an obvious gap.

full=True prepends one slot, stamped -1, holding each cell's residual sum of squares — the one quantity that separates the two. It goes first so the axis stays ascending, which keeps the coordinate-aware members usable on the result.

In [5]:

diagnosed = rhum.curvefit('level', line, [0.0, 50.0], full=True)
np.asarray(diagnosed.coords['param']).tolist()

diagnosed = rhum.curvefit('level', line, [0.0, 50.0], full=True) np.asarray(diagnosed.coords['param']).tolist()

Out[5]:

[-1, 0, 1]

In [6]:

planes = np.asarray(diagnosed.read_array()).reshape(
    *diagnosed._band_dim_sizes, *rhum.shape[-2:]
)
residual = planes[:, 0]

# Threshold the residual to build a quality mask, instead of trusting every cell equally.
float(np.nanmin(residual)), float(np.nanmedian(residual)), float(np.nanmax(residual))

planes = np.asarray(diagnosed.read_array()).reshape( *diagnosed._band_dim_sizes, *rhum.shape[-2:] ) residual = planes[:, 0] # Threshold the residual to build a quality mask, instead of trusting every cell equally. float(np.nanmin(residual)), float(np.nanmedian(residual)), float(np.nanmax(residual))

Out[6]:

(0.0, 125.32856582594769, 5968.285447482554)

In [7]:

# `sel` reaches the residual by its stamp, so the slot is addressable rather than positional.
only_residual = np.asarray(diagnosed.sel(param=-1).read_array())
only_residual.shape, float(only_residual.reshape(12, *rhum.shape[-2:])[0, 18, 36])

# `sel` reaches the residual by its stamp, so the slot is addressable rather than positional. only_residual = np.asarray(diagnosed.sel(param=-1).read_array()) only_residual.shape, float(only_residual.reshape(12, *rhum.shape[-2:])[0, 18, 36])

Out[7]:

((12, 1, 37, 72), 223.9107042761254)

A caller mistake is kept separate from a data problem. An initial guess outside its bounds, bounds whose lower exceeds their upper, or a model taking fewer parameters than p0 supplies is refused up front with a message naming the argument — rather than running the fit, failing in every cell, and returning an all-NaN raster with a warning that blames the data.

Seeing the quality mask¶

The residual is a (y, x) field per time step, so it maps. Plotting goes through cleopatra, which every plot in pyramids delegates to.

In [8]:

%matplotlib inline

worst = float(np.nanpercentile(residual[0], 99))
mask_me = residual[0] > worst
print(f'cells above the 99th percentile of residual: {int(mask_me.sum())}')

%matplotlib inline worst = float(np.nanpercentile(residual[0], 99)) mask_me = residual[0] > worst print(f'cells above the 99th percentile of residual: {int(mask_me.sum())}')

cells above the 99th percentile of residual: 27

In [9]:

glyph = (
    diagnosed.sel(param=-1)
    .isel(time=0)
    .plot(title='curvefit residual, first time step')
)
type(glyph).__name__

glyph = ( diagnosed.sel(param=-1) .isel(time=0) .plot(title='curvefit residual, first time step') ) type(glyph).__name__

Out[9]:

'ArrayGlyph'

rolling_exp — an infinite window¶

rolling averages a fixed window with equal weights, which costs it the first window - 1 steps to no-data and makes it forget everything older than the window. rolling_exp weights every earlier step instead, decaying geometrically, so it answers from step 0 and never fully forgets. On a short or irregular satellite series that is the difference between discarding the earliest scenes and keeping them.

It needs no new dependency: pandas.ewm does the work, and a cube's cells are already columns, so one call smooths the whole cube in C at about 19 µs per cell. (xarray's own rolling_exp is built on numbagg and raises ImportError without it.)

It is also step-based, not coordinate-aware — it counts steps along the axis and never reads its coordinates, exactly as rolling does.

In [10]:

smoothed = rhum.rolling_exp('time', 0.3)
series_raw = np.asarray(rhum.read_array()).reshape(12, 4, *rhum.shape[-2:])[
    :, 0, 18, 36
]
series_smooth = np.asarray(smoothed.read_array()).reshape(12, 4, *rhum.shape[-2:])[
    :, 0, 18, 36
]

np.round(series_raw, 2).tolist()

smoothed = rhum.rolling_exp('time', 0.3) series_raw = np.asarray(rhum.read_array()).reshape(12, 4, *rhum.shape[-2:])[ :, 0, 18, 36 ] series_smooth = np.asarray(smoothed.read_array()).reshape(12, 4, *rhum.shape[-2:])[ :, 0, 18, 36 ] np.round(series_raw, 2).tolist()

Out[10]:

[79.25, 80.5, 80.25, 80.25, 82.25, 80.75, 81.0, 85.0, 86.5, 83.5, 78.75, 79.75]

In [11]:

np.round(series_smooth, 2).tolist()

np.round(series_smooth, 2).tolist()

Out[11]:

[79.25,
 79.99,
 80.11,
 80.16,
 80.92,
 80.86,
 80.91,
 82.21,
 83.55,
 83.53,
 82.07,
 81.36]

The smoothed series answers at step 0, where rolling(3) would still be no-data, and it approaches each change rather than snapping to it.

Four reductions are offered — mean, sum, std and var. xarray also exposes cov and corr, which need a second cube and therefore a different signature, so they are not here. std and var are no-data at step 0, one step having no spread.

In [12]:

spread = rhum.rolling_exp('time', 0.3, how='std')
np.round(
    np.asarray(spread.read_array()).reshape(12, 4, *rhum.shape[-2:])[:4, 0, 18, 36], 3
).tolist()

spread = rhum.rolling_exp('time', 0.3, how='std') np.round( np.asarray(spread.read_array()).reshape(12, 4, *rhum.shape[-2:])[:4, 0, 18, 36], 3 ).tolist()

Out[12]:

[nan, 0.884, 0.591, 0.443]

A gap does not blank the rest of the series — it contributes nothing and the decay carries through it — but the gap step itself stays a gap. That is the one place this departs from raw ewm, which fills it with the carried value: a finite, plausible number with nothing marking it unobserved. For a cloud-masked scene that would be fabricated data.

In [13]:

gappy = NetCDF.from_array(
    np.array([10.0, np.nan, 10.0, 20.0]).reshape(4, 1, 1),
    geo_ref=GeoReference(geo=(0.0, 1.0, 0.0, 1.0, 0.0, -1.0), epsg=4326),
    variable_name='t',
    dims=ExtraDimensions(name='time', values=[0.0, 1.0, 2.0, 3.0]),
).get_variable('t')

np.round(gappy.rolling_exp('time', 0.5).read_array().ravel(), 4).tolist()

gappy = NetCDF.from_array( np.array([10.0, np.nan, 10.0, 20.0]).reshape(4, 1, 1), geo_ref=GeoReference(geo=(0.0, 1.0, 0.0, 1.0, 0.0, -1.0), epsg=4326), variable_name='t', dims=ExtraDimensions(name='time', values=[0.0, 1.0, 2.0, 3.0]), ).get_variable('t') np.round(gappy.rolling_exp('time', 0.5).read_array().ravel(), 4).tolist()

Out[13]:

[10.0, nan, 10.0, 16.1538]

to_dict — the whole structure as plain Python¶

Every other export goes to a file or to a foreign object. This one goes to dict, list, str and numbers, which is what inspection, diffing two cubes, JSON transport and asserting on structure in a test all want.

The schema is a superset of xarray's: dims, coords, data_vars and attrs keep their xarray meanings, and everything GDAL needs rides under one extra pyramids key. xarray's schema has no slot for a CRS or a geotransform, so matching it exactly would make the round trip silently drop the georeferencing.

In [14]:

payload = cube.to_dict()
sorted(payload)

payload = cube.to_dict() sorted(payload)

Out[14]:

['attrs', 'coords', 'data_vars', 'dims', 'pyramids']

In [15]:

payload['dims']

payload['dims']

Out[15]:

{'time': 12, 'level': 4, 'lat': 37, 'lon': 72}

In [16]:

payload['pyramids']

payload['pyramids']

Out[16]:

{'schema': 1,
 'epsg': 4326,
 'geotransform': [-2.5, 5.0, 0.0, 92.5, 0.0, -5.0],
 'spatial_dims': ['lat', 'lon'],
 'no_data_value': {'rhum': 32766.0}}

A time axis keeps its CF attributes, which is the metadata whose loss would change what the numbers mean while every value still matched.

In [17]:

payload['coords']['time']['attrs']

payload['coords']['time']['attrs']

Out[17]:

{'units': 'hours since 1-1-1 00:00:0.0', 'calendar': 'standard'}

data=False describes the cube without reading any values — it makes no read_array call at all. The shape and dtype come from the declared metadata, and each band dimension's coordinate vector is still materialised, that being the structure the mode exists to report. Coordinates lose their data key here too, so the two modes stay consistent.

In [18]:

structure = cube.to_dict(data=False)
entry = structure['data_vars']['rhum']
entry['dims'], entry['shape'], entry['dtype'], 'data' in entry

structure = cube.to_dict(data=False) entry = structure['data_vars']['rhum'] entry['dims'], entry['shape'], entry['dtype'], 'data' in entry

Out[18]:

(['time', 'level', 'lat', 'lon'], [12, 4, 37, 72], 'float64', False)

from_dict — back again¶

The contract is the round trip, not the parse: values, coordinates, CRS, geotransform and no-data all come back.

In [19]:

rebuilt = NetCDF.from_dict(payload)
back = rebuilt.get_variable('rhum')

(
    back._band_dim_names,
    back._band_dim_sizes,
    back.epsg == rhum.epsg,
    tuple(back.geotransform) == tuple(rhum.geotransform),
)

rebuilt = NetCDF.from_dict(payload) back = rebuilt.get_variable('rhum') ( back._band_dim_names, back._band_dim_sizes, back.epsg == rhum.epsg, tuple(back.geotransform) == tuple(rhum.geotransform), )

Out[19]:

(('time', 'level'), (12, 4), True, True)

In [20]:

# Asserted rather than printed: with outputs stripped before commit, `--nbval-lax` compares
# nothing, so a cell that merely *displays* `False` would ship unnoticed. This one raises.
assert np.allclose(
    np.asarray(back.read_array(), dtype='float64'),
    np.asarray(rhum.read_array(), dtype='float64'),
    equal_nan=True,
)
print('every value round-tripped')

# Asserted rather than printed: with outputs stripped before commit, `--nbval-lax` compares # nothing, so a cell that merely *displays* `False` would ship unnoticed. This one raises. assert np.allclose( np.asarray(back.read_array(), dtype='float64'), np.asarray(rhum.read_array(), dtype='float64'), equal_nan=True, ) print('every value round-tripped')

every value round-tripped

Three things do not round-trip, in each case because from_array is the only constructor and takes neither: per-variable attributes, arbitrary global attributes (globals arrive as a fixed set of CF fields), and non-gridded variables, which have no raster plane and are dropped with a warning naming them.

And one caveat on JSON. Python's json emits NaN and Infinity and reads them back, so a payload round-trips between Python processes. Those tokens are not RFC-8259, so a strict parser rejects them: a gappy raster's payload is Python-JSON, not wire-JSON.

In [21]:

import json

text = json.dumps(payload)
restored = NetCDF.from_dict(json.loads(text))
len(text), restored.variable_names

import json text = json.dumps(payload) restored = NetCDF.from_dict(json.loads(text)) len(text), restored.variable_names

Out[21]:

(2438702, ['rhum'])