Broadcasting a band axis, and coordinate roles¶

Four members, plus one rule in the shared engine they rest on. The rule is the interesting part: a single-band operand now applies to every band of the other, so a time series times a land mask finally has a spelling.

Member	What it does
combine (the rule)	a length-one band axis stretches across the other operand
broadcast_like	give this cube another cube's band layout, repeating cells
broadcast_equals	compare two cubes after broadcasting both
set_coords	mark existing variables as CF auxiliary coordinates
reset_coords	demote auxiliary coordinates back to data variables

Two things this is not. It is not xarray alignment: joining axes by a coordinate index needs an index model, which pyramids does not have, so a shared axis at two real lengths is refused rather than outer-joined. And the spatial axes never broadcast - same_grid stays an exact check, so a one-cell raster is never stretched over a scene. Only the band (non-spatial) axes take part.

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.dataset import Dataset
from pyramids.dataset import GeoReference as DatasetGeoReference
from pyramids.netcdf import ExtraDimensions, GeoReference, NetCDF

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

# A CF store of 17 monthly ERA5-style variables on a shared 12-step `time` axis.
ERA = SAMPLES / 'cf__20v__1d3-3d17__y-desc.nc'

from pathlib import Path import numpy as np from pyramids.dataset import Dataset from pyramids.dataset import GeoReference as DatasetGeoReference from pyramids.netcdf import ExtraDimensions, GeoReference, NetCDF SAMPLES = Path('../../../examples/data/netcdf/samples') # A CF store of 17 monthly ERA5-style variables on a shared 12-step `time` axis. ERA = SAMPLES / 'cf__20v__1d3-3d17__y-desc.nc'

The cube we start from¶

p2t is one of the 17 variables: twelve monthly steps on a 73x144 grid. _band_dim_names is the non-spatial part of that - the axes broadcasting acts on.

In [2]:

nc = NetCDF.read_file(str(ERA))
t2m = nc.get_variable('p2t')

t2m.dimension_names, t2m.shape, t2m._band_dim_names

nc = NetCDF.read_file(str(ERA)) t2m = nc.get_variable('p2t') t2m.dimension_names, t2m.shape, t2m._band_dim_names

Out[2]:

(['time', 'latitude', 'longitude'], (12, 73, 144), ('time',))

A land mask on the same grid, with no band dimensions at all - a single plane:

In [3]:

geo = GeoReference(geo=nc.geotransform, epsg=nc.epsg)
mask_cells = (np.arange(73 * 144).reshape(73, 144) % 3 != 0).astype(float)
mask = NetCDF.from_array(mask_cells, geo_ref=geo, variable_name='mask').get_variable(
    'mask'
)

mask.dimension_names, mask.shape, mask._band_dim_names

geo = GeoReference(geo=nc.geotransform, epsg=nc.epsg) mask_cells = (np.arange(73 * 144).reshape(73, 144) % 3 != 0).astype(float) mask = NetCDF.from_array(mask_cells, geo_ref=geo, variable_name='mask').get_variable( 'mask' ) mask.dimension_names, mask.shape, mask._band_dim_names

Out[3]:

(['y', 'x'], (1, 73, 144), ())

The rule: a single band applies to every band¶

Before this, cube * mask raised - combine refused every band-count mismatch. Now the one-band operand stretches across all twelve steps. Nothing is tiled: the mask is read once and stretched with np.broadcast_to, which is a view.

In [4]:

masked = t2m * mask

masked.band_count, masked._band_dim_names, masked._band_dim_sizes

masked = t2m * mask masked.band_count, masked._band_dim_names, masked._band_dim_sizes

Out[4]:

(12, ('time',), (12,))

The result's band count, band names and band dimensions come from the multi-band operand, whichever side it is on, so the two orders agree:

In [5]:

reversed_order = mask * t2m

(
    reversed_order.band_count == masked.band_count,
    reversed_order._band_dim_names == masked._band_dim_names,
)

reversed_order = mask * t2m ( reversed_order.band_count == masked.band_count, reversed_order._band_dim_names == masked._band_dim_names, )

Out[5]:

(True, True)

Every binary operator and comparison routes through the same engine (Analysis._combine), so they all gained this at once - arithmetic, comparisons, and Dataset as well as NetCDF:

In [6]:

dataset_geo = DatasetGeoReference(top_left_corner=(0.0, 3.0), cell_size=1.0, epsg=4326)
stack = Dataset.from_array(np.arange(27.0).reshape(3, 3, 3), geo_ref=dataset_geo)
plane = Dataset.from_array(np.array([[1.0, 0.0, 1.0]] * 3), geo_ref=dataset_geo)

(stack * plane).band_count, (stack > plane).band_count

dataset_geo = DatasetGeoReference(top_left_corner=(0.0, 3.0), cell_size=1.0, epsg=4326) stack = Dataset.from_array(np.arange(27.0).reshape(3, 3, 3), geo_ref=dataset_geo) plane = Dataset.from_array(np.array([[1.0, 0.0, 1.0]] * 3), geo_ref=dataset_geo) (stack * plane).band_count, (stack > plane).band_count

Out[6]:

(3, 3)

Only 1-vs-n broadcasts. A 2-vs-3 pair is still refused with the same message as before, and so is a grid mismatch - those are genuine errors, not shapes to stretch.

where gains the rule through the same gate, but with an asymmetric twist: it answers in the receiver's shape, so a one-band condition covers a whole stack while a condition with more bands than the raster is refused.

In [7]:

stack.where(np.ones((3, 3), bool), 0.0).band_count

stack.where(np.ones((3, 3), bool), 0.0).band_count

Out[7]:

3

broadcast_like - materialise the layout¶

The operator route above never builds the stretched cells. When you need the broadcast result as a cube - to write it to a file, concat it, or hand it to to_xarray - broadcast_like materialises it:

In [8]:

lifted = mask.broadcast_like(t2m)

lifted.dimension_names, lifted.shape, lifted._band_dim_sizes

lifted = mask.broadcast_like(t2m) lifted.dimension_names, lifted.shape, lifted._band_dim_sizes

Out[8]:

(['time', 'y', 'x'], (12, 73, 144), (12,))

The added axis takes the donor's size and its coordinate values, and the cells simply repeat along it - the first step and the last are identical:

In [9]:

cells = np.asarray(lifted.read_array(squeeze=True))

np.array_equal(cells[0], cells[11])

cells = np.asarray(lifted.read_array(squeeze=True)) np.array_equal(cells[0], cells[11])

Out[9]:

True

An axis both cubes carry behaves by its length: length one here against the donor's n is stretched and takes the donor's stamps; equal lengths are left alone with their own stamps; an axis this cube has and the donor lacks is kept as it is.

In [10]:

one_step = NetCDF.from_array(
    np.full((1, 73, 144), 7.0),
    geo_ref=geo,
    variable_name='single',
    dims=ExtraDimensions(name='time', values=[0.0]),
).get_variable('single')

stretched = one_step.broadcast_like(t2m)
stretched._band_dim_sizes, np.asarray(stretched.coords['time']).tolist()[:4]

one_step = NetCDF.from_array( np.full((1, 73, 144), 7.0), geo_ref=geo, variable_name='single', dims=ExtraDimensions(name='time', values=[0.0]), ).get_variable('single') stretched = one_step.broadcast_like(t2m) stretched._band_dim_sizes, np.asarray(stretched.coords['time']).tolist()[:4]

Out[10]:

((12,), [898476, 898482, 898500, 898506])

A shared axis at two real lengths - neither of them one, so there is nothing to stretch - is refused rather than joined. That is the alignment pyramids does not do, and the error says so.

Note the cost, and one gotcha. The repeats are real cells: a mask broadcast over twelve steps holds twelve times the data. And on a container, broadcast_like goes through the same rebuild as transpose and pad, which does not carry variable attributes - including the CF coordinates attribute the next section writes. Promote after broadcasting, not before.

broadcast_equals - compare after broadcasting¶

equals is exact, so the mask and the twelve-step cube built from it are not equal. Broadcast both first and they are - which is precisely what broadcast_equals answers:

In [11]:

(
    mask.equals(lifted),
    mask.broadcast_equals(lifted),
)

( mask.equals(lifted), mask.broadcast_equals(lifted), )

Out[11]:

(False, True)

It broadcasts both sides and reorders the right onto the left's axes, so the answer is symmetric, and a cube is always reflexively equal to itself. Unequal cells stay unequal - broadcasting changes the shapes, never the values:

In [12]:

(
    mask.broadcast_equals(mask),
    lifted.broadcast_equals(mask) == mask.broadcast_equals(lifted),
    mask.broadcast_equals(t2m),
)

( mask.broadcast_equals(mask), lifted.broadcast_equals(mask) == mask.broadcast_equals(lifted), mask.broadcast_equals(t2m), )

Out[12]:

(True, True, False)

set_coords - mark variables as coordinates¶

CF already expresses the coordinate/data duality: a variable is a coordinate of another because that other names it in its coordinates attribute. pyramids already read that - classify_variables assigns the auxiliary-coordinate role, and variable_names / data_vars filter themselves by it - so set_coords is the write side, and needed no new partition.

The store starts with 17 data variables:

In [13]:

len(nc.variable_names), nc.meta_data.cf.classifications['lsp']

len(nc.variable_names), nc.meta_data.cf.classifications['lsp']

Out[13]:

(17, 'data')

In [14]:

promoted = nc.set_coords('lsp')

len(promoted.variable_names), promoted.meta_data.cf.classifications['lsp']

promoted = nc.set_coords('lsp') len(promoted.variable_names), promoted.meta_data.cf.classifications['lsp']

Out[14]:

(16, 'auxiliary_coordinate')

lsp left the data variables, exactly as it does in xarray - so data_vars and everything built on it (merge, concat, apply, reduce, to_dataframe, to_xarray) follows without being told. The variable itself was never moved or copied: only its role changed, and it still reads:

In [15]:

(
    'lsp' in promoted.data_vars,
    promoted.get_variable('lsp').shape,
)

( 'lsp' in promoted.data_vars, promoted.get_variable('lsp').shape, )

Out[15]:

(False, (12, 73, 144))

Like rename_dims / assign_coords / drop_dims, this is non-mutating - the receiver is untouched and a new container comes back:

In [16]:

len(nc.variable_names)

len(nc.variable_names)

Out[16]:

17

The reference is written onto every gridded data variable whose non-spatial axes cover the promoted one's - that is what CF's coordinates attribute means ("these variables label my cells"). Promoting a dimension is refused rather than silently ignored: a dimension's coordinate is its same-named array by CF convention, and the netCDF driver cannot reassign it.

reset_coords - demote them back¶

The exact inverse. It removes the references, deleting the attribute outright rather than leaving an empty string behind:

In [17]:

restored = promoted.reset_coords('lsp')

len(restored.variable_names), restored.meta_data.cf.classifications['lsp']

restored = promoted.reset_coords('lsp') len(restored.variable_names), restored.meta_data.cf.classifications['lsp']

Out[17]:

(17, 'data')

With no names at all it sweeps the whole store, so a round trip through several promotions comes back where it started:

In [18]:

two_promoted = nc.set_coords(['lsp', 'cp'])

len(two_promoted.variable_names), len(two_promoted.reset_coords().variable_names)

two_promoted = nc.set_coords(['lsp', 'cp']) len(two_promoted.variable_names), len(two_promoted.reset_coords().variable_names)

Out[18]:

(15, 17)

Recap¶

• A single-band operand broadcasts across the other in every operator and comparison, on Dataset and NetCDF alike. Stretched as a view, never tiled. Only 1-vs-n; the spatial axes never broadcast.
• broadcast_like(other) materialises that layout as a cube - adding the donor's axes with its coordinates, stretching a length-one axis, leaving an equal-length one alone.
• broadcast_equals(other) broadcasts both sides, then compares; symmetric, and stricter than it looks - only the shapes move.
• set_coords(names) / reset_coords(names=None) write and remove the CF coordinates references. Container operations, non-mutating, and exact inverses.

What stays refused, and why: alignment. Joining a shared axis at two real lengths, join= semantics, and assign_coords of a brand-new coordinate all need an index model. Broadcasting a length-one axis is NumPy's rule and needs no index - which is why it shipped and aligning did not.

For the other structural members see the "Editing dimensions & variables" notebook; for the wider tour of the class, the NetCDF anatomy notebook.