Skip to content

Migration guide#

This guide helps downstream packages that depend on pyramids migrate across breaking changes. It is organized by subpackage, and under each subpackage by the release that introduced the changes.

isinstance(x, NetCDF) keeps working throughout, so most code needs no change. Everything below either emits a DeprecationWarning (discoverable with the command in the next section) or is a hard behavior change (called out explicitly).

Finding what affects you (one command)#

Run your own test suite against the new pyramids with deprecation warnings turned into errors — this surfaces every deprecated call site:

python -W error::DeprecationWarning -m pytest

Or, in code under test:

import warnings
warnings.simplefilter("error", DeprecationWarning)

This catches the deprecations; the hard behavior changes do not warn — search for them manually.

base#

unreleased#

abfs:// now means Azure Data Lake Storage Gen2, not Blob. It mapped to GDAL's /vsiaz/ (Blob) and now maps to /vsiadls/ (Gen2), matching what abfs means everywhere else in the Azure and Hadoop ecosystems — it is the Azure Blob File System driver, which is Gen2.

  • If you used abfs:// against a Gen2 account, this is the fix: you were being routed to the Blob handler and now get the right one, with directory semantics.
  • If you used abfs:// against a flat Blob account, switch to az://, which is unchanged and still names Blob. Credentials are the same AZURE_STORAGE_* family either way.
  • abfss:// is accepted too, if you prefer to be explicit about TLS. There is no adls:// scheme: the registered Azure Data Lake name is adl:// (Gen1), which GDAL does not handle.
  • The canonical abfss://<filesystem>@<account>.dfs.core.windows.net/<path> form is understood: the filesystem becomes the container and the account is dropped from the path, because GDAL takes it from configuration. The rewrite is deterministic — the same URL always yields the same /vsi* path — and if the URL names an account that disagrees with the configured one, a UserWarning says which account will actually be read. Set AZURE_STORAGE_ACCOUNT (or AccountName= inside AZURE_STORAGE_CONNECTION_STRING) to match, or use the bare abfs://<filesystem>/<path> form.

/vsiadls/ is now recognised as a remote, network-backed path. Previously is_remote() and is_network_backed() both returned False for a Gen2 path, so it was classified as a local file — opening still worked, but anything reasoning about the path (credential handling, archive chaining) got the wrong answer.

Archive chaining now covers every network handler, including for raw /vsi* paths. A raster inside a zip or tar is rewritten to /vsizip/<handler>/... for /vsiadls/, /vsioss/ (Alibaba), /vsiswift/ (OpenStack), /vsihdfs/ and /vsiwebhdfs/, not only for S3, Google Cloud, Azure Blob and /vsicurl/. Four of those have no URL scheme at all, so the raw /vsi* spelling is the only way to name them — paths already in /vsi* form now go through the chaining too, where before they were returned untouched. If you were chaining by hand, the manual prefix is no longer needed, and a hand-built /vsizip//vsioss/... still works unchanged.

Credentials are unchanged. /vsiadls/ reads the same AZURE_STORAGE_ACCOUNT / AZURE_STORAGE_ACCESS_KEY / AZURE_STORAGE_SAS_TOKEN / AZURE_STORAGE_CONNECTION_STRING family as /vsiaz/, and honours AZURE_NO_SIGN_REQUEST identically (verified: with it set both handlers skip the instance-metadata credential lookup). CloudConfig gains an azure_no_sign_request=True flag so anonymous Azure reads have the same knob the S3 side already had.

0.47.0#

Dataset.epsg no longer reports EPSG:4326 for a raster that has no CRS. It previously substituted WGS 84 whenever the projection was empty, so an ungeoreferenced grid claimed a georeference it did not have — and that claim propagated into to_file, to_crs, bounds and alignment checks. It now returns None, matching GDAL, rasterio and geopandas, none of which substitute a default.

What changed, and what did not:

  • No CRS and no evidence (a GeoTIFF or ASCII grid with an empty projection) → epsg is None. Operations that genuinely need a CRS raise CRSError naming the fix instead of proceeding.
  • A CRS with no EPSG authority (geostationary, spherical-earth GRIB) → unchanged by this release. Note that epsg is not uniformly None here: NetCDF reports None for a geostationary grid (see #706), while a spherical-earth GRIB still resolves through epsg_from_wkt's 4326 default. Read crs when you need the authoritative answer for such a grid.
  • A CF NetCDF with degrees_east / degrees_north axes and no grid_mapping → still 4326. CF leaves the datum implicit for these and the whole ecosystem reads them as WGS 84, so this is a reading of the file's metadata rather than an assumption. Three things narrow it:
    • Only a coordinate variable's units are evidence — one that declares axis, is named in a coordinates attribute (a curvilinear grid's 2-D lat/lon), or carries a conventional coordinate name. A data variable in degrees_east (a wind direction, a solar angle) is not evidence.
    • A horizontal axis carrying a linear unit (m, km) or plain degrees/rad vetoes it, so a projected, rotated-pole or geostationary grid is not relabelled on the strength of auxiliary lat/lon arrays. A vertical axis does not veto — its units say nothing about the horizontal frame. Vertical axes are recognised from CF's axis: Z, falling back to conventional names (depth, lev, …) only for a file that declares no axes at all.
    • The grid's own extent must fit in the lon/lat range. A metre-scale grid is not lat/lon whatever its metadata says.

A CRS with no EPSG code now reports epsg is None too. Previously any projection the EPSG register does not name — an orthographic or geostationary projection, a rotated pole, a spherical-earth GRIB GEOGCS — fell back to 4326. That claimed WGS 84 for grids that are not WGS 84: an orthographic frame is not lat/lon at all, and a spherical datum differs from the WGS 84 ellipsoid by up to ~20 km.

The CRS itself is unaffected — ds.crs still returns the WKT, crs_spec() falls back to it, and reprojection, cropping and alignment all keep working. Only the code is absent, because there is not one. If your code reads ds.epsg on such a grid, read ds.crs (or crs_spec(ds.epsg, ds.crs)) instead.

If you relied on the old default, set the CRS explicitly — dataset.epsg = <code> in process, or gdal_edit.py -a_srs EPSG:<code> <file> on disk. To find affected code, look for .epsg used without a None check, and for the dataset.epsg or dataset.crs idiom, which is now crs_spec(dataset.epsg, dataset.crs).

0.46.1#

import pyramids no longer configures logging. Previously, importing the package ran Config(), which added a handler to the root logger, called root.setLevel(logging.INFO), pinned six third-party loggers to WARNING, and printed Logging is configured. to stderr. A library reconfiguring the root logger silently takes over logging for the whole host application, so that is gone.

Change Kind What you do
Import no longer installs a console handler behavior Configure logging yourself, or ask pyramids for it
Root logger untouched (no handler, no setLevel) behavior Nothing — this is the fix
Config / LoggerManager default level=INFO → None behavior Pass a level to get console output
Logging is configured. announce moved from INFO to DEBUG behavior Nothing
Third-party loggers only quieted on request behavior Pass quiet_third_party=True to restore

If you saw pyramids log lines and want them back, pick one:

# Option A — you own logging (recommended for applications)
import logging

logging.basicConfig(level=logging.INFO)
import pyramids  # pyramids records propagate to your handler
# Option B — let pyramids install its coloured console handler
from pyramids.base.config import Config

Config(level="INFO")

Option B sets propagate = False on the pyramids logger, so records are emitted once by pyramids' own handler rather than twice (once more through your root handler). Do not combine it with Option A expecting both to print.

Option B is reversible: a bare Config() removes the handlers pyramids installed and restores propagation. That matters for test suites — pytest's caplog attaches to the root logger, so while propagate is False a caplog assertion on a pyramids.* logger captures nothing and passes vacuously. If a dependency opts in on your behalf, call Config() to hand the namespace back.

If you relied on the third-party quieting, pass it explicitly:

Config(level="INFO", quiet_third_party=True)  # pins fiona/rasterio/shapely/matplotlib/urllib3/osgeo to WARNING

Exception types on three dtype helpers changed to a descriptive ValueError, from the incidental error that leaked out of an empty table lookup. Only affects code catching the old type:

Call Before After
numpy_to_gdal_dtype(<unmapped dtype>) IndexError ValueError
gdal_to_numpy_dtype(gdal.GDT_Unknown) AttributeError ValueError
gdal_to_ogr_dtype(<complex band>) TypeError ValueError

dataset#

unreleased#

CF-packed data is now unpacked on read, and the keyword is unpack. Hard change, silent for the values — nothing raises. A raster that declares scale_factor / add_offset (GDAL's GetScale / GetOffset) now reads back in physical units as float64 instead of raw stored counts.

v = NetCDF.read_file(path).get_variable("VHM0")
v.read_array().max()
# before -> 1435.0   raw counts
# after  ->   14.35  metres, the honest answer

stats() and plot() follow the same rule; before, neither could produce a physical value at all. stats() does it without re-reading — CF packing is affine, so min, max and mean take raw * scale + offset while std takes |scale|.

A raster that is not packed is untouched. The identity is short-circuited: same values, same dtype, no copy, no float64 promotion. It counts in every spelling it arrives in — GDAL answers None for a band that was never packed, while Dataset.scale / .offset report that as 1.0 / 0.0, and a file may store the identity outright. If your rasters are unpacked — nearly all GeoTIFFs — nothing changes.

scaled= is gone; use unpack=. The two names were the same concept in two places (IO.read_array(scaled=) for rasters, NetCDF.read_array(unpack=) for CF variables), both defaulting off and both routing to the same primitive. There is one keyword now, defined on Dataset so NetCDF inherits it, and no alias:

ds.read_array(scaled=True)   # before
ds.read_array()              # after — unpacked by default
ds.read_array(unpack=False)  # after — the raw store, when you want it

apply() no longer truncates, and no longer corrupts a packed raster. Two separate faults. It built the output at the source band's type, so a float-valued function on an integer raster was written back as integers — 1435 * 0.01 stored 14, not 14.35. It also dropped the band's scale / offset while keeping the raw values, so apply(lambda a: a) — the identity — changed a packed variable's physical value from 2.5 to 100.0. The output now takes the function's own result type, and reads reaching apply are already physical, so the packing is genuinely spent rather than silently discarded. This holds for both arms: the out-of-core elementwise=True path sizes its per-tile buffer from the same probe as the destination band, where it used to round the result back into the source's type before the wider band ever saw it. A scalar-only callable (math.sqrt) is probed the way it is run, through np.vectorize, so it is not truncated either.

A function whose result GDAL has no type for (an object array) still writes at the source type, as before.

An operation that copies the store keeps the packing; one that computes new values spends it. This is the rule that makes the default safe, and it is worth knowing if you subclass or read the internals. A crop, a seam join across the antimeridian, the container-wide fan-out — these move stored bytes, so they read with unpack=False and carry scale_factor / add_offset onto the result: a cropped int16 CMEMS variable stays int16 and still declares its recipe, rather than quadrupling in size as float64. apply, from_array and the reductions produce values the packing has already been spent on, so their results declare none.

Two consequences you may see:

  • a packed raster survives crop, to_crs (both the plain and maintain_alignment=True forms), resample and align with its scale and offset intact. Several of these dropped the packing before: the cutline border-trim silently discarded it (a Dataset lost it outright; a NetCDF variable kept it only in a Python attribute), and the three ReprojectImage paths never carried it at all;
  • for a NetCDF variable, _scale / _offset are authoritative for the read when set, and the band's own GetScale / GetOffset are the fallback. Everything that applies the packing resolves it through that one rule — both read arms, stats, get_histogram, point / read_part / preview / read_overview_array / get_tile, set_variable and the streaming transforms — so none of them can answer in different units for the same band. On a NetCDF variable the public scale / offset properties answer from the same rule, so a sel() result -- whose band declares nothing and whose recipe lives only in Python -- reports the factor its values are actually unpacked with, and assigning to either updates that recipe rather than being outranked by it.

no_data_value stays a stored value, and that is the second half of the contract. CF puts _FillValue in the packed datatype, GDAL reports it that way, and it is what gets written back — so it is left alone. The consequence is that you cannot compare it against a default read:

values = ds.read_array()          # physical: the gap reads as -98.49
ds.no_data_value                  # stored:   (-9999.0,)
values == ds.no_data_value[0]     # never true on a packed raster

Inside the library every such comparison is now made in one consistent unit. Where a mask is wanted, it is built against the stored counts and the values unpacked afterwards; where the sentinel is handed on as a value — to cleopatra, to a feature table, to an ASCII grid's NODATA_value header — it is expressed in the physical units the array holds. The two pick out the same cells, because the packing is affine and a zero or non-finite scale_factor is refused. If you were doing that comparison yourself, use read_array(masked=True), which masks in stored units and hands back physical values, or read with unpack=False and mask before you scale.

A transform always works in physical units; the destination only decides how the result is stored. apply, map_blocks and stream_transform hand your function physical values, and a gap in the source comes back as the result's declared sentinel whatever the function does to it. stream_transform(out=...) into a destination that declares a recipe — including the in-place out=ds form — packs the result with that destination's recipe, (value - offset) / scale, rounded to its dtype, so source and destination need not share a packing. An in-place lambda t: t * 2 therefore doubles the physical values, not the stored counts.

merge_rasters inherits its no-data from the sources instead of defaulting to 0. A hard behavior change, and the reason is that 0 is real data in most rasters worth merging: sea-level land in an elevation or bathymetry model, the zero crossing of an anomaly or difference raster. The old default stamped 0 on the mosaic, so every such cell disappeared from masked reads, from stats() and from any plot — a Copernicus DEM mosaic of the Nile Delta lost 45.9% of the scene to it. It also discarded whatever the sources declared, in both directions: their marker was not inherited, and a cell that was no-data arrived in the mosaic as a real value.

no_data_value now defaults to inheriting:

  • Sources agree — the mosaic declares that value. Nothing else changes.
  • Sources disagree — the first one wins and a UserWarning names both, instead of silently picking.
  • No source declares one, and the sources leave a gap — a marker is chosen, because those uncovered pixels are still written and leaving them undeclared makes them read as measurements. The choice is the value such a pixel already holds where the output's data type can store it — NaN for a floating mosaic, and always for method="min", "max" and "sum", which write Float64. An integer band has no NaN, so there the marker is instead a value that band can store and that the mosaic's own cells do not use (65535 for a UInt16 scene holding small numbers, -9999 for a signed one), and the compositing step is filled with it so the gaps really hold what the output declares. Only a mosaic whose data uses every value its type could spare is written without a marker, and that warns.
  • No source declares one, the sources tile their area, and the band has no NaN — nothing is declared, because there is nothing to mark. Choosing a value would mean reading every source to prove it unused, which is the expensive half of the paragraph above and buys nothing when no pixel is uncovered. Whether they tile is answered from the footprints, without reading a pixel. This case is the integer z-order mosaic only: a floating one still declares NaN, and so do method="min", "max" and "sum" whatever the footprints, because NaN costs nothing to prove unused.

Previously all of this was 0: on an integer band 0 was also the only thing hiding the NaN that init puts in the VRT, which the band cannot store.

Related behaviour that changed with it:

  • The reduction methods' uncovered pixels changed value, not just their declaration — they used to be 0.0 and are now the inherited marker (e.g. -9999.0) or NaN.
  • method="min", "max" and "sum" over integer sources were wrong before and are now right. Each source was warped into its own data type before being folded, where the NaN marking its uncovered area was rounded to 0; those zeros then won every fmin and were added by every sum, so a mosaic of integer tiles that did not tile contiguously was corrupted whatever no_data_value said: min collapsed to that value everywhere, sum had it added into every cell, and max kept its covered data but lost the gaps. Sources are now warped into the Float64 the reduction writes.
  • The marker now reaches the pixels it marks. Uncovered pixels used to hold init (NaN by default) while the mosaic declared something else, so a mosaic inheriting -9999 declared -9999 over gaps holding NaN and read_array(masked=True) masked none of them. The gaps are now filled with whatever is declared.
  • An inherited value the output band cannot store is replaced, with a warning. Integer sources can declare NaN — Dataset.no_data_value reports it for an integer raster asked for one — and GDAL then refused the marker outright ("Nodata value was not set to output band"), leaving the mosaic unmarked.
  • Each source's own no-data now marks that source's holes. n is an override, and its default "nan" was passed to gdal.BuildVRT as one, replacing every source's declaration. Tiles that declared different values therefore had all but the winner's holes composited as real measurements — a tile declaring -32768 beside one declaring -9999 read its holes back as -32768. The default now means "no override". If you were relying on n's default to ignore NaN cells in sources that declare nothing, pass n= explicitly.
  • A numeric init is a preference, not an instruction. While the marker is being inherited, init is offered to the sentinel search as the preferred candidate and taken only if the mosaic's own cells do not hold it. init=0 over tiles containing real zeros would otherwise have become the declared marker and masked them — #1086 again, through a different argument. Being passed over warns.
  • no_data_value is validated. A value that names no number used to reach gdal.Translate, which answered "Nodata value was not set to output band" and wrote no marker at all; it now raises ValueError. The string "inherit" is accepted as the default, since that is how the default renders in help() and in the API docs.

Migrating:

  • If you relied on the 0 default, pass it explicitly: merge_rasters(src, dst, no_data_value=0) restores the old declaration, and is worth a second look — it masks every genuine 0 in your inputs. It does not restore the old pixels: the three changes in the third bullet below apply whatever you pass, so the gaps now hold 0 rather than init. Add init= alongside it to keep the old fill.
  • If you want no marker at all, pass no_data_value=None. Both write paths honour that, and neither stamps anything.
  • If you were passing no_data_value= already, the marker is unaffected — but the pixels can still change. Three of the bullets above apply whatever you passed: each source's own holes are now skipped rather than composited, the integer reductions were corrected, and — new — the pixels no source covers are now filled with your marker instead of with init, so that the marker you asked for actually masks them. Pass init= alongside it to keep the old fill. On a z-order merge these three are the only changes; a min/max/sum merge of integer sources changes more, because it was wrong.
  • The value is read from the sources' band 1, in the order you pass them, so ordering decides a disagreement.
  • DatasetCollection.merge forwards the same default, so it answers identically.

count_domain_cells refuses an out-of-range band with ValueError, not IndexError. Asking for a band the raster does not have used to reach self._ds.no_data_value[band] and surface IndexError: tuple index out of range, which names neither the band nor the dataset. It now raises the same message the rest of the package gives for a bad band, as convert_units and get_overview_dataset already did:

>>> dataset.count_domain_cells(band=5)
ValueError: band 5 is out of range for a 1-band dataset.
  • Change except IndexError to except ValueError if you wrapped the call. Nothing else about the count changes, and a valid band behaves exactly as before.
  • Only the too-large case changes. A negative band already raised ValueError with this same message, from the read layer further down, so count_domain_cells(band=-1) is unaffected. It is simply refused earlier now, before any data is read.
  • The new domain_area refuses the same bands the same way, so the pair stays consistent.

The web-service readers accept a bbox that crosses the antimeridian. Additive if you pass an ordinary box; a hard change if you relied on minx > maxx being rejected. Dataset.from_wcs, from_wms, from_wmts and from_ogc_coverages used to raise ValueError: bbox must have minx < maxx and miny < maxy for (170, -10, -170, 10) — while Dataset.crop had served exactly that box for a long time, by splitting at the 180 degree seam. The readers now split too: two requests, one per side, concatenated into one raster whose geotransform continues past the seam (170 .. 180 then 180 .. 190) instead of jumping back to -180.

  • A transposed longitude pair is no longer caught, and that is the price. A network reader has no grid to measure the bbox against — it is fetching the grid — so (6, 51, 5, 52), a west/east pair entered the wrong way round, now reads as a 359 degree wrap and issues two large requests rather than raising. If you were relying on that ValueError to catch swapped arguments, check the bbox yourself before the call. Latitude is unaffected: miny >= maxy is still an error, because there is no seam in latitude.
  • Two wrapping bboxes are refused rather than approximated, both on the from_wms / from_wmts side: one with a projected crs, where there is no 180 degree meridian and minx > maxx really is inverted; and one with a corner outside -180 .. 180, where "west of the seam" stops meaning anything.
  • from_wms's size is the width of the whole result, divided between the two halves at one shared resolution, with the seam snapped to the nearest pixel boundary. The two widths sum to the size you asked for, odd or even, so a wrapping request returns the image dimensions a non-wrapping one would. A wrapping request needs at least 2 pixels of width, since one pixel cannot straddle the seam.
  • from_wms(resolution=...) is now capped at 25,000 px per axis. Hard change, and it applies to ordinary non-wrapping reads too: from_wms(bbox=(5, 51, 6, 52), resolution=0.00001) used to return a 100,000 x 100,000 image and now raises ValueError. The cap exists because accepting a wrap made a transposed bbox silently span 359 degrees — at a fine resolution that becomes two GetMap requests of a few hundred thousand columns each, which no server will serve. size=(width, height) is not capped: a size you state outright cannot be amplified by a mistake in bbox, so it is still taken verbatim. If you were relying on a very fine resolution, pass size= instead.
  • from_wmts measures the seam offset in the layer's native CRS rather than assuming 360 — about 40,075,017 m for a Web Mercator layer — and only computes it when there is a seam, so a non-wrapping read is byte-identical to before.
  • from_ogc_coverages sizes the two halves together when no resolution is given: the combined span is capped once and the resulting pixel size applied to both. Sizing each half against the cap on its own would give them different pixel sizes and row counts, which cannot be concatenated at all.
  • A wrapping from_wcs(direct=True) read should pass resolution=. Direct mode has no descriptor to read, so the two halves are independent GetCoverage requests; with a resolution they are snapped onto one lattice, and without one nothing constrains the server to grid them alike. A mismatch is refused with antimeridian halves were produced at different resolutions rather than stitched into a misaligned raster.
  • A seam read through from_wcs or from_ogc_coverages returns a plain Dataset, even when called on a subclass, because those two merge by rebuilding through Dataset.from_array. from_wms and from_wmts stitch into a GDAL handle and wrap it with the class they were called on, so a subclass survives there. Only relevant if you call these classmethods on a Dataset subclass; nothing in pyramids does.

Two rasters can now be combined directly, and Dataset gained arithmetic operators. Additive — nothing that worked before behaves differently. ds.combine(other, func) runs a binary function over two aligned rasters and returns a Dataset on the left operand's grid, and -, +, *, / between two rasters are thin wrappers over it. Previously this meant reading both bands into numpy, computing, rebuilding a GeoReference from one of the inputs and wrapping the result with from_array — the step where georeferencing gets lost.

  • The operands must already share a grid; a mismatch raises AlignmentError rather than resampling. The new ds.same_grid(other) predicate answers the question in advance, and align() is the explicit fix.
  • A scalar operand is accepted on either side: ds * 2, 2 * ds, 20 - ds and ds >= 5 all work. A scalar takes the same route as a raster operand, so ds * 2 and ds * other agree on band count, dtype and sentinel. That is deliberately not the same as ds.apply(lambda v: v * 2), which transforms one band and keeps the source's sentinel — reach for apply when you want that contract. bool and complex scalars are still refused.
  • A cell that is no-data in either operand is no-data in the result. An integer result's sentinel is derived against the values func computed, so no in-range number is claimed as a gap by the arithmetic that produced it, and one that masked nothing declares no sentinel at all. A floating result always declares NaN: a cell func computed as NaN (0/0 in a normalised difference) has no value, so it is a gap and the result says so. Pass no_data_value= to choose one, or no_data_value=None for no masking.
  • sum(rasters) works: __radd__ absorbs the 0 that sum() seeds with, returning a copy so a one-element sum never aliases its input. Absorption is symmetric and is not limited to the integer seed — adding any real scalar equal to zero, or multiplying by any equal to one, short-circuits to a copy on either side, so ds + 0, 0.0 + ds and ds * 1 are copies too. ds - 0 is absorbed as well, on the right only — 0 - ds negates. ds / 1 is not: true division widens an integer band to float64, so it is no no-op. Absorption is what keeps the two spellings of one commutative expression from disagreeing: routed through combine, ds + 0.0 would widen an int16 band to float64 where 0.0 + ds is byte-identical, and ds + 0 would drop the band's declared sentinel. Every other scalar computes through combine, so 1 + ds returns a raster rather than raising. The copy is a real cost: sum() and math.prod() each materialise one extra full raster that functools.reduce(operator.add, rasters) does not, which matters near the memory limit.
  • combine is whole-array: both operands are read in full. For rasters near the memory limit use apply(elementwise=True) or read_array(chunks=).
  • The module-private _same_grid helper in pyramids.dataset.dataset moved to Spatial.same_grid, faced on Dataset. It was never public, but anyone importing it directly must switch to a.same_grid(b).
  • <, <=, > and >= between two rasters, or between a raster and a real scalar, return a Byte mask (1/0, and 255 wherever either operand was no-data). == and != are not overridden — they stay identity-based, so Dataset remains usable in assert, in sets and as a dict key; use a.combine(b, np.equal) for the mask.
  • math.prod(rasters) folds like sum(rasters); both absorb their identity scalar from either side, and compute for every other scalar.

bool(ds) now raises — replace if ds: with if ds is not None:. Hard change. A raster holds one value per cell, and a comparison between two rasters is itself a raster, so there is no honest single truth value. Reduce explicitly — bool(np.asarray((a >= b).read_array()).all()) — and use is not None for presence checks. numpy, pandas and xarray all refuse a truth value for the same reason.

The refusal is what keeps the new comparison operators safe. Before this release a >= b raised TypeError, so if a >= b: and sorted(rasters) were errors; making them return a raster without also refusing a truth value would have turned those errors into silently wrong answers, which is why the two ship together.

np.<ufunc>(ds, ...) now raises TypeError. Hard change. Dataset sets __array_ufunc__ = None, so numpy defers instead of treating a raster as an opaque object. Previously np.array([0.0]) + ds returned an object array holding a raster copy — a silent wrong answer rather than an error. Use the operators, or combine.

Dataset.dtype reports numpy's spelling, so a Byte raster reads uint8 rather than byte. Hard change, silent — the property still returns one string per band, but two of the catalog's names moved, two half-precision types became reachable, and one band type that used to return a value now raises. The names are numpy's wherever numpy has one, which is what STAC's raster:bands[].data_type wants; GDAL's complex-integer and complex-half types have no numpy name, so those keep the catalog's own spelling rather than all collapsing onto complex64.

Band type Before After
GDT_Byte ['byte'] ['uint8']
GDT_CFloat64 ['complex-float64'] ['complex128']
GDT_Float16 IndexError from the catalog lookup ['float16']
GDT_CFloat16 IndexError from the catalog lookup ['complex-float16']
GDT_Unknown the catalog's empty name — a float NaN, not a str ValueError naming the code
everything else (uint16, int32, float32, complex-int16, complex-float32, …) unchanged unchanged

byte and complex-float64 are the two that break working code: a ds.dtype[0] == "byte" branch stops matching, and anything that feeds dtype into STAC, a manifest or a filename emits a different string. Search for string comparisons against "byte" and "complex-float64" and update them.

Feeding the value back is better off than it was. numpy_to_gdal_dtype("byte") returned 14 — GDT_Int8, a different band type — so numpy_to_gdal_dtype(ds.dtype[0]) on a Byte raster silently produced the wrong code; "uint8" gives 1, which is right. "complex-float64" raised TypeError there and "complex128" gives 11. The two half-precision rows are an addition, not a break — they raised before. GDT_Unknown turns a value into an exception, but the value it returned was a NaN that violated the property's own list[str] contract, so nothing could have used it.

Dataset.copy() with no path returns a dataset in "write" mode, even when the source was read-only. Hard change, silent — the copy used to inherit the source's access mode:

ds = Dataset.read_file("raster.tif")     # read_only
clone = ds.copy()

# Before
clone.access          # 'read_only'   <- a write raised ReadOnlyError
# After
clone.access          # 'write'

An in-memory copy is a fresh MEM raster that nothing else references, so refusing writes on it protected nothing — the source file is untouched either way — while a copy() taken specifically to get a mutable working copy of a read-only file could not be written to, and the ReadOnlyError advised reopening the source with read_only=False, which is not what the caller wanted. Copying to a path is unchanged in spirit: a format that supports Create gives "write", and a write-by-copy-only format (.png, .jpg / .jpeg, .jp2 / .j2k, .asc) still gives "read_only", because CreateCopy hands back a read-only dataset for those. If you relied on the copy being read-only as a guard, guard it yourself — there is no parameter to ask for the old mode.

A no-data value the band cannot store is no longer replaced with the dtype's maximum. Hard change, silent on the read path and loud on the write path. A NaN sentinel cannot be stored in an integer band, and the unsigned types wider than a byte used to answer that by substituting their own maximum: a uint16 band asked for NaN came back declaring 65535. Every genuinely-65535 cell in it — the saturated end of a scaled reflectance product, the void value of a DEM — was then out of domain, and align handed the raster to gdal.ReprojectImage, which rewrote those real 65535s to 65534 to keep them distinguishable from the sentinel. uint8 was excluded from the substitution for exactly this reason (255 is white in 8-bit imagery); the wider types now follow the same rule, so there is one rule across all of them.

dtype and no-data Before After
uint16 + NaN (65535,) (nan,)
uint32 + NaN (4294967295,) (nan,)
uint64 + NaN (18446744073709551615,) (nan,)
uint8, signed, floating + NaN (nan,) (nan,) — unchanged
any dtype, no_data_value=None (None,) (None,) — unchanged

Two knock-on effects to look for:

  • change_no_data_value(None) on a uint16 / uint32 / uint64 raster now raises NoDataValueError, where it used to succeed and leave the dtype maximum behind. uint8 and the signed types already raised. Pass a sentinel the band can hold, or leave the band without one.
  • A raster that declares no storable sentinel now really has no no-data, so count_domain_cells, apply, fill and the statistics treat every cell as data. That is the honest answer for a band with no free value, and it is what a uint8 raster already did.
  • Dataset.create(..., no_data_value=np.nan) returns a different canvas, not just a different declaration. Creating a band fills it with the sentinel it just declared, so the substitution used to leave every cell out of domain; with nothing storable to declare there is nothing meaningful to fill with, and the band keeps the driver's zeros:
create(rows=3, columns=4, dtype=…, no_data_value=np.nan) Before After
uint16 declares 65535, holds 65535, 0 domain cells declares NaN, holds 0, 12 domain cells
uint32 declares 4294967295, holds that, 0 domain cells declares NaN, holds 0, 12 domain cells

If you used that call as an out-of-domain canvas and wrote only the cells you cared about, the untouched ones are now real zeros that every consumer counts as data. Create the raster with a sentinel the dtype can hold instead, and the old behaviour returns unchanged.

What has not changed is the repair applied when a caller asks for a sentinel that overflows the band: Dataset.from_array(uint8_array) still warns that the default -9999 is out of range and stores 255. There the caller asked for a sentinel and it did not fit, so picking a storable one honours the request; the removed substitution fired when the caller asked for no sentinel, and answered a question nobody posed.

crop derives the value it writes into masked cells, and declares it on the result. Hard change, and the fix for two bugs. Cropping has to put something in the cells the mask excludes, and it used to write the source's declared sentinel — which is why a multi-band uint8 or int16 crop raised ValueError: cannot convert float NaN to integer outright, and why a uint16 crop silently relied on the substitution above.

The cropped output now carries a fill chosen against the band's own data: the band's sentinel when it can store one, otherwise the first value that fits the dtype and occurs nowhere in the band (NaN, then -9999, then the dtype's extremes, then a scan of the remaining range for the narrow integer types). That value is set as the output's no_data_value, so the result says which cells are absent instead of leaving a bare 0 or a fabricated maximum indistinguishable from a measurement. A band that holds every candidate raises NoDataValueError naming the fix rather than reclassifying a real observation as a gap.

A band that declares no sentinel is derived for too, on this path. That is not inventing a property the data lacked — the crop is what makes those cells absent — and what it replaces was collide-or-crash. Measured on the previous release, a raster-mask crop of a band declaring nothing:

dtype single band multi-band
uint16 / uint32 TypeError wrote and declared the dtype maximum
uint8 / int16 TypeError TypeError
floating NaN, declared NaN, declared

Because the output now declares a sentinel the trim recognises, a raster-mask crop also removes rows and columns that the mask excluded entirely — so the result can be smaller than the mask's extent. That is not new behaviour, only newly reachable: a float32 band, whose NaN was always storable, was trimmed the same way before. A 4x4 raster masked along its first row and column comes back 3x3 for every dtype now, where the integer cases used to raise instead.

The one case that completed is the one this issue exists for: _check_no_data_value turned the band's None into 65535, wrote it into the excluded cells and declared it, so every genuinely-65535 cell in the band became out-of-domain on the output. Every other integer case raised TypeError: int() argument must be ... not 'NoneType', single- and multi-band alike. A floating band is unchanged — it gets NaN, which is what that path already wrote and declared.

# uint8 source asked for a NaN sentinel, data in 1..7
out = src.crop(mask)
out.no_data_value      # (255.0,) -- derived, and declared
# Before: ValueError: cannot convert float NaN to integer (multi-band)

# uint8 source declaring nothing at all
out = src.crop(mask)
out.no_data_value      # (255.0,) -- derived, and declared
# Before: TypeError: int() argument must be ... not 'NoneType'

A polygon (cutline) crop is the same story from the other side, with one deliberate difference: a band that declares no sentinel is left alone there. That path does not write the excluded cells itself — GDAL fills them — so stamping a -dstnodata would put a sentinel on the output that the source never had, which the netCDF fan-out would then carry onto a rebuilt variable. Only a band with a sentinel it cannot store gets one.

For that band, the border outside the cutline was GDAL's default 0 — a real observation on most bands, and declared nowhere. The trim that removes the fully-outside rows and columns works by reading the sentinel back, so it recognised nothing there and left the border in place: an integer raster came back larger than the same crop of a float one. It now gets the derived fill, is trimmed to the polygon, and declares what its border holds.

crop(polygon) of an 8x8 raster asked for a NaN sentinel Before After
uint8 / int16 6x6, border of undeclared 0 4x4, border declared
uint16 / uint32 4x4, border declared 65535 4x4, border declared — now derived
int64 6x6, border of undeclared 0 4x4, border declared
uint64 6x6, border of undeclared 0 6x6 — unchanged, see below
float32 (a storable NaN) 4x4 4x4 — unchanged

uint16 and uint32 were already trimmed and already declared, because the substitution being removed had given them 65535 / 4294967295 before the crop ever ran. What changes for them is only where the value comes from: derived against the band's data rather than fabricated, so it can no longer be a value the band holds.

uint64 is left alone deliberately. -dstnodata reaches GDAL as text it parses into a C double, and that dtype's maximum does not survive the trip — it arrives rounded to 2**63, which GDAL announces and then writes into the pixels. Declaring the value we asked for would describe cells holding something else, so no fill is offered and the border stays undeclared.

The declared sentinel of a well-formed raster is unchanged: a band that already declares a storable value keeps it, its data is never read to make the decision, and the warp options it produces are the ones it produced before.

A crop can now come back larger, for any raster whose cells come close to its sentinel. Hard change, silent, and unrelated to whether a fill was derived. The trim that removes wholly-excluded rows and columns asked is_no_data, whose operational rtol of 0.001 treats everything within a part per thousand of the sentinel as no-data; it now asks is_stored_no_data, the tolerance the storage forces. Cells that merely lie near the sentinel are data again, so fewer rows qualify for removal and the result only ever grows:

float32 declaring -9999, border holding -9995 Before After
crop(mask=polygon) 2x2 6x6
crop(bbox=…) 2x2 6x6

The old shape was reached by deleting real observations — elevation, depth and scaled-reflectance products sit in exactly that window — so the new one is the correct answer. But code that indexes, reprojects or asserts against a crop's shape will see it change. It applies to both crop routes and to bands with an ordinary, storable, declared sentinel.

create_from_array is now from_array, and takes a GeoReference. Hard change, no deprecation alias — the old name and the old flat keywords are gone. The same rename applies to UgridDataset.create_from_arrays -> from_arrays (still plural: it takes three arrays, not one).

The four flat georeferencing keywords are replaced by one geo_ref argument, and driver_type is gone entirely — the driver is now derived from the path extension:

# Before
Dataset.create_from_array(arr, geo=GEO, epsg=4326)
Dataset.create_from_array(arr, top_left_corner=(0, 10), cell_size=0.05)
Dataset.create_from_array(arr, driver_type="GTiff", path="out.tif")
Dataset.create(cell_size=0.05, rows=r, columns=c, dtype="float32", bands=1,
               top_left_corner=(0, 10), epsg=4326)
Dataset.create_empty(rows, cols, geo=GEO, epsg=4326, driver_type="MEM")
UgridDataset.create_from_arrays(node_x, node_y, faces)

# After
from pyramids.dataset import Dataset, GeoReference

Dataset.from_array(arr, geo_ref=GeoReference(geo=GEO, epsg=4326))
Dataset.from_array(arr, geo_ref=GeoReference(top_left_corner=(0, 10), cell_size=0.05))
Dataset.from_array(arr, geo_ref=GeoReference(geo=GEO), path="out.tif")  # driver from the extension
Dataset.create(rows=r, columns=c, dtype="float32", bands=1,
               geo_ref=GeoReference(top_left_corner=(0, 10), cell_size=0.05, epsg=4326))
Dataset.create_empty(rows, cols, geo_ref=GeoReference(geo=GEO, epsg=4326))
UgridDataset.from_arrays(node_x, node_y, faces)

The module-level pyramids.netcdf.engines.variables.create_from_array is renamed from_array too. It carries no underscore but is an engine internal — reach it through NetCDF.from_array, its public facade.

GeoReference is importable from pyramids.dataset, pyramids.netcdf and pyramids.netcdf.array_options — all three resolve to the same class.

Why there is no alias. The signature and the name were both changing; keeping the old spelling working would have meant migrating twice. Both failure modes are loud — the old name raises AttributeError, a flat keyword raises TypeError — and because NetCDF.from_array now declares its real parameters instead of *args, **kwargs, mypy catches both statically, which the old facade prevented.

driver_type removed from create_from_array and create_empty. path alone decides memory-vs-disk, and the extension names the format, so the parameter could only agree with path or contradict it. path=None gives an in-memory raster; path="out.tif" a GTiff; path="out.nc" a netCDF — the last of which was rejected on both routes before: the default driver_type="MEM" by the .tif suffix check, and an explicit driver_type="netcdf" by the unconditional COMPRESS=LZW the netCDF driver refuses. Removing the parameter and gating that option is what makes it work. Creation options still require a path, which is the invariant the old driver_type="GTiff" guard was really enforcing.

Formats that GDAL can only write by copy (.png, .jp2) now raise FileFormatNotSupportedError naming the extension and the driver, instead of failing inside GDAL. An extension the driver catalog does not know still raises DriverNotExistError.

merge_rasters and stack_bands now refuse write-by-copy-only destinations. Both have two internal write paths and only one of them can produce such a format, so accepting it would make the destination's legality depend on method (or on whether the inputs share a dtype) — the same defect as letting an unrelated argument pick the format. Both take the stricter answer, which refuses .png, .jpg / .jpeg, .jp2 / .j2k, .asc and .vrt. The one this costs is .asc, which merge_rasters could write before through its z-order (method="first" / "last") path. Write a GTiff and convert it.

Dataset.to_file now accepts more extensions, and one of them can lose data. This method is untouched by the rename, but it reads the same driver catalog, so correcting the catalog's extension rows widened what it accepts. .tiff, .png, .jpg, .jpeg, .img, .jp2 and .j2k previously raised DriverNotExistError and now write a file. The sharp case is dtype. PNG and JPEG store 8-bit data, so a float32 raster written to .png is converted to Byte — values clipped, fractional parts gone. GDAL reports that only as a RuntimeWarning, so pyramids now raises it to a DtypeNarrowingWarning naming both dtypes and the driver. JPEG2000 is stricter still: .jp2 / .j2k accept only their own dtypes and raise rather than converting, so a float32 raster there warns and then fails with FailedToSaveError — write those from a uint8 / int16 source.

import warnings
from pyramids.errors import DtypeNarrowingWarning

warnings.simplefilter("error", DtypeNarrowingWarning)   # make the conversion fail loudly

Note the deliberate asymmetry with the constructors: to_file("x.png") writes (PNG supports CreateCopy) while from_array(path="x.png") raises FileFormatNotSupportedError (PNG does not support Create). Whether a format is reachable depends on how the call writes, not only on the extension.

You can no longer name a driver whose extension the catalog does not list. This is the one capability the driver_type removal costs. On main, driver_type="EHdr", path="out.bil" worked; now .bil raises DriverNotExistError and there is no escape hatch. If you need a format the catalog does not carry, write a GTiff (or build in memory) and convert, or open an issue asking for the extension to be added. Sibling spellings of a catalogued format now resolve alike: .tiff as well as .tif, .nc4 as well as .nc, .jpg as well as .jpeg, and .j2k as well as .jp2. Resolving is not the same as being writable by these constructors: all four JPEG/JPEG2000 spellings resolve and are then refused with FileFormatNotSupportedError, because they are write-by-copy-only formats. (.jpeg did not resolve on main at all — the catalog carried it with a leading dot, so the lookup never matched.)

Dataset.create no longer requires a CRS. epsg was a required positional; omitting it raised TypeError. It is now a field of geo_ref defaulting to 4326, so a call that says nothing about the CRS silently stamps WGS 84 instead of refusing. from_array and create_empty always behaved this way — only create changes, and it changes from "state the CRS" to "we assume WGS 84". Pass epsg explicitly, or GeoReference(..., epsg=None) for a deliberately CRS-less raster.

Dataset.create's positional order changed, and its CRS handling widened. cell_size is gone from slot 0, so rows and columns each shift up one — but geo_ref is required and keyword-only, so a ported positional call fails loudly with TypeError rather than silently misreading its arguments. Separately, create used to take epsg: int straight through sr_from_epsg; it now shares the helper the other constructors use, so GeoReference(epsg=None) yields a raster with no CRS, and non-EPSG CRS strings are accepted.

create_empty is the only constructor where geo_ref is optional. A header-only allocation often does not care where it sits, so omitting geo_ref — or passing one that carries no transform at all, such as GeoReference(epsg=3857) — keeps the identity transform. A partially specified reference is not covered by that convenience: a top_left_corner without a cell_size (or the reverse) raises, exactly as it does in from_array and create, rather than silently discarding the half you supplied.

RasterLike and RasterBase changed shape. create_from_array is from_array on the pyramids.base.protocols.RasterLike protocol and on the RasterBase ABC; bands_values is gone, and variable_name moved down to NetCDF, where it actually applies. If you subclass RasterBase, rename your override and adopt the new signature. RasterLike is @runtime_checkable, so this one bites at runtime too: an isinstance(obj, RasterLike) guard on a stale duck type carrying create_from_array flips from True to False and silently takes the else-branch, rather than raising anywhere you can see it.

pyramids.netcdf.array_options.GeoTransform is gone, with no re-export. It was a local alias for the plain 6-tuple, and pyramids.dataset.GeoTransform is a richer NamedTuple of the same shape but a different type, so re-exporting that under the old name would have silently changed what the name meant. If you imported the alias for annotation, use pyramids.base.georeference.GeoTransformTuple (the structural 6-tuple) or the pyramids.dataset.GeoTransform NamedTuple — whichever you actually meant.

NetCDF.from_array returns a Container. It always did; the annotation said NetCDF. Container adds no public API over NetCDF, so nothing changes at runtime — the type is just honest now.

create_overviews now refuses a plain VRT whose description is not a path. It raises OverviewTargetError instead of returning normally. That is a new exception, importable as from pyramids.errors import OverviewTargetError, which subclasses ValueError, so an existing except ValueError around the call keeps catching it — catch the new type to tell "this dataset can never work" apart from "these arguments were wrong". recreate_overviews raises it for the same shape, where it previously reported a misleading ReadOnlyError advising a reopen that a handle with no path cannot perform — or, when the VRT exposed no levels at all, returned normally with a UserWarning saying to call create_overviews() first. That advice is now refused too, so a pipeline that called recreate_overviews() defensively across a mix of handles and ignored the warning goes from a silent no-op to a raise. A plain VRT owns no pixel storage, so GDAL can only write its overviews to an external .ovr sidecar named after the dataset description — and when that description is not a path, the sidecar landed as a file called literally .ovr in the process's working directory, attached to nothing, while the levels the handle was already exposing were dropped. The call reported success and silently did neither thing it promised.

Three descriptions are refused: an empty one, a blank one, and any that begins with < once stripped — in practice an inline VRT XML document passed to Dataset.read_file(...), which previously surfaced as a raw GDAL RuntimeError naming a file whose name was the whole document. No real path starts with <.

These calls hand back such a handle, and so start raising:

  • Dataset.get_overview_dataset(...) — only its lazily described form, taken when the parent can be reopened by name. When the parent cannot be (a from_array raster, for instance), the method materialises a MEM level instead, and that still builds.
  • anything wrapping a gdal.Translate(..., format="VRT") or gdal.BuildVRT("", ...) result kept in memory.
  • NetCDF.get_variable(...) when the classic view comes back in index space — an irregularly spaced coordinate defeats the geotransform guess, and the view is then wrapped in a plain pathless VRT. A view over a regular grid is not VRT-wrapped and is unaffected.
  • Dataset.wrap_longitude() on a file-backed source, whose lazy roll is a plain pathless VRT. On an in-memory source the roll is a MEM raster and is unaffected. This one previously produced the damage above verbatim — a stray .ovr and no levels — so the refusal is a fix, not a regression.
  • Dataset.to_zarr(..., overview_factors=[...]) on any handle the refusal covers, since it builds the pyramid levels through create_overviews. The target is checked pre-flight, so the call leaves no store behind, and before the compute argument, so a call that is also compute=False reports the refusal rather than the ValueError. Without overview_factors the refusal does not apply; a description-less VRT then fails where it always did, in to_zarr's base-array write.

Write the view out first and build the overviews on the saved raster:

view.to_file("level.tif")
view.close()  # on Windows an open handle keeps the parent file locked
saved = Dataset.read_file("level.tif", read_only=False)
saved.create_overviews(overview_levels=[2, 4])

Unaffected by the create_overviews refusal: warped VRTs keep their overviews in RAM and need no sidecar. These produce one — Dataset.warped_view(...), Dataset.to_crs(...) in its warping form, Dataset.crop(mask, touch=False) with a vector mask, and the lazy georeference / orthorectify forms. (to_crs(..., maintain_alignment=True) and crop with a raster mask are different paths: both return a MEM dataset, exempt as a non-VRT rather than as a warped one.) Tests cover warped_view and to_crs; the others are exempt by the same root-element check. A VRT with a real path (including under /vsimem/) names its sidecar after that path, and MEM rasters are not VRTs at all. Warped VRTs are not exempt from the recreate_overviews refusal — see the next entry.

recreate_overviews now raises OverviewTargetError instead of ReadOnlyError when a VRT computes the levels. GDAL reports every unwritable overview target with the same CPLE_NoWriteAccess it uses for a genuinely read-only dataset, so the old message told callers to reopen with read_only=False whenever the write was refused. For a level a VRT computes rather than stores, that advice cannot help and often cannot even be followed. Two shapes are affected:

  • a warped VRT — its levels land on VRTWarpedRasterBands, which are never writable. A warped view taken from a writable parent fails identically, and a pathless view has no path to reopen. create_overviews() still builds a warped view's levels; only in-place regeneration is refused.
  • a plain VRT that inherits its levels from the source it wraps — the common gdal.Translate(..., format="VRT") or mosaic over an already-overviewed raster. GDAL serves each such level from an implicit read-only VRT it builds for the level, not from the handle you opened, so the read-only dataset in its message is one read_only=False cannot reach — holding the source open writable does not help either. Give the VRT its own sidecar with create_overviews(), or regenerate on the source raster instead.

The two are told apart by who owns the level: a level owned by a VRT is computed, one owned by a real raster (the dataset itself for an internal overview, the .ovr GTiff for an external one) is stored.

ReadOnlyError is now raised only when a reopen is still worth trying — a stored level on a handle that is itself open read-only. It is not a promise that reopening will succeed: a VRT serving an explicit <Overview> owns a real, on-disk-writable .ovr, yet GDAL opens VRT sources read-only and refuses however the parent was opened. What has changed is that a dataset already open for writing never gets told to reopen — that shape now raises OverviewTargetError too, since the access mode demonstrably is not the blocker.

file_i counts only real files, refuses a negative index, and now works for a tar. Hard change, silent for zip and loud for tar. file_i is public on Dataset.read_file and NetCDF.read_file, and it selects a member of a compressed archive. Three behaviours moved. Measured on a zip and a tar each holding a sub/ directory entry, sub/a.asc and b.asc:

Call Before After
zip, file_i=0 .../probe.zip/sub/ — the directory entry .../probe.zip/sub/a.asc
zip, file_i=1 .../probe.zip/sub/a.asc .../probe.zip/b.asc
zip, file_i=2 .../probe.zip/b.asc FileFormatNotSupportedError naming the count and members
zip, file_i=-1 .../probe.zip/b.asc — Python's from-the-end FileFormatNotSupportedError
zip, file_i=-5 IndexError: list index out of range FileFormatNotSupportedError
tar, any file_i /vsitar/<archive> — the index was ignored the member at that index
  • Directory entries no longer consume an index. A zip written with explicit directory entries shifted every member by one and made file_i=0 name a directory GDAL cannot open. If you compensated by adding one, remove the adjustment.
  • A negative index is refused rather than counted from the end. file_i=-1 never meant "the last member" — it happened to work through Python's list indexing while -5 raised a bare IndexError from the same line. Both now raise FileFormatNotSupportedError, which names how many members the archive holds and lists them. To read the last member, pass its index: len(members) - 1.
  • A tar now honours file_i at all. Every index resolved to the bare /vsitar/<archive> path before, and GDAL picked a member itself. A tar read that relied on that pick can now land on a different file — pass the index you want, or open the member path directly.

to_zarr writes the band's sentinel as the array's fill_value. Hard change, silent, and it changes bytes on disk. The store's array metadata carried fill_value: 0.0 regardless of the raster's no-data, which was recorded only in a pyramids-specific no_data_value attribute. It now carries the real sentinel, in both the Zarr fill_value slot and a CF _FillValue attribute. Written from a float32 raster whose no-data is -9999.0:

# Before                                    # After
data/zarr.json  fill_value = 0.0            data/zarr.json  fill_value = -9999.0
attrs: no_data_value = [-9999.0]            attrs: _FillValue = -9999.0, no_data_value = [-9999.0]

Dataset.read_file(store).no_data_value      Dataset.read_file(store).no_data_value
  -> (0.0,)                                   -> (-9999.0,)

The cell values are untouched — only the declared sentinel moved — so pyramids reading its own store back is the fix: the raster's no-data survives the round trip instead of being replaced by 0.0. The cost is borne by every other reader. Any consumer that honours the Zarr fill_value or the CF _FillValue — zarr-python and GDAL do, and so does anything layered on them — reads -9999.0 from a store written by the new version where it read 0.0 from one written by the old, and masks accordingly. Stores already on disk are unaffected until they are rewritten. If a downstream pipeline depended on 0.0 being the fill value, rewrite the store or set the reader's fill explicitly.

No-data masking in the analysis engine now takes its tolerance from the band's dtype. Hard change, silent. count_domain_cells, apply, fill, footprint, plot_histogram and to_image all ask one question — does this cell hold the band's no-data sentinel — and each asked it with a relative tolerance of its own: rtol=1e-3 in count_domain_cells and apply, 1e-6 in fill, 1e-5 in the two renderers. So one band could be counted one way and drawn another, and every one of them discarded ordinary cells that merely lay near the sentinel. A relative window scales with the sentinel: rtol=1e-3 around -9999 masks everything down to -10009, and rtol=1e-5 around an int32 sentinel of 2e9 masks everything within 20 000 of it.

The tolerance is now the band's own — none at all for an integer or boolean band, whose sentinel is stored exactly, and single precision's eps (about 1.2e-7 relative) for a floating band of any width, which is the slack a sentinel picks up passing through float32 storage or a driver's decimal text. numpy.isclose's default atol=1e-8 is dropped with it, so a sentinel of 0 no longer swallows genuinely small cells.

Two consequences of "the band's own" worth stating, because the obvious reading of each is the wrong one:

  • The float tolerance does not narrow further for a wide band, and does not widen for a narrow one. It is single precision's eps for float16, float32 and float64 alike. Taking each dtype's own eps would make a float16 band the loosest of the three — one ULP of half precision is 9.8e-4, a window of ±9.8 around a -10000 sentinel, which swallows both of that sentinel's neighbouring representable values.
  • An integer sentinel is compared as an exact integer, not through a float. That matters only at the 64-bit limits, where it is the whole answer: float(2**63 - 1) rounds up past int64's maximum, so a band whose sentinel is its own dtype maximum would be judged unable to hold it and mask nothing at all — with fill then overwriting the very cells it was told to leave alone. A 64-bit maximum is a sentinel a caller reasonably declares, and the one a crop derives for a band that has none, so this is a path a 64-bit raster reaches rather than a corner of one.

For a caller, a cell within ~0.1% of the sentinel that used to be no-data to count_domain_cells / apply is now data; nothing that was masked before is unmasked less accurately, the change only ever keeps cells. On the test suite's own rasters nothing moves at all — all 473 fixture bands carrying a concrete sentinel (of 519 scanned) mask identically under the old rule and the new one — so this surfaces only on data that holds values close to its own sentinel. Ask pyramids.base._domain.is_no_data(arr, nodata, rtol=...) directly if you want the old, looser window.

RasterBase.get_x_lon_dimension_array / get_y_lat_dimension_array are removed. Hard change. They were thin wrappers over GeoTransform.x_axis / GeoTransform.y_axis; use those directly instead — e.g. ds.transform.x_axis(ds.columns) and ds.transform.y_axis(ds.rows), or GeoTransform(*geo).y_axis(rows) when you hold a bare geotransform. Two behaviour differences to know when you switch:

  • The axis honours the signed step. x_axis / y_axis read the signed geotransform[1] / geotransform[5], so a positive geotransform[5] (south-up) ascends where the old get_y_lat_dimension_array forced a descending axis (it negated its argument and documented it as a positive pixel height). A caller that passed abs(geotransform[5]) must now pass the signed value. This is what fixes Dataset.y / Dataset.lat reporting coordinates outside a south-up raster's extent, and Dataset.bbox is likewise normalised now (always [min_x, min_y, max_x, max_y], no longer inverted for a south-up or east-left grid). The common north-up case is unchanged.
  • Values can differ in the last ULP. The axis is origin + (index + 0.5) * step, not a repeated element-wise addition, so a printed coordinate can move (10.35 where the walk produced 10.350000000000001). Compare axis arrays with numpy.allclose, not == or a repr.

0.48.0#

recreate_overviews now rebuilds a band's levels in one cascading pass. Each deeper level is decimated from the level above rather than from the full-resolution band, matching what create_overviews (BuildOverviews) has always done — previously recreate_overviews rebuilt every level directly from the source. Nothing warns; the call still succeeds.

Level 0 is unaffected. A level >= 1 keeps its previous values only where the resampling survives being applied twice:

  • nearest — unchanged everywhere; GDAL does not cascade it.
  • average, rms — unchanged on a floating-point band with no no-data value. An integer band picks up per-level rounding of up to 1 DN, and a no-data gap changes the result regardless of dtype.
  • cubic, cubicspline, bilinear, lanczos, gauss, mode — changed on ordinary data. mode can move by whole classes, because the mode of modes is not the mode.

There is no way to get the old per-level values back: no API rebuilds a deep level directly from the source any more. If your pipeline depends on them, rebuild the levels you need yourself from Dataset.read_array() and write them out.

0.47.0#

Operations that need a CRS now refuse instead of assuming one. Each raises CRSError naming the operation and how to fix it, where previously the missing CRS was silently filled in with WGS 84:

  • Dataset.to_crs(...) — reprojection has no source frame to warp from.
  • Dataset.align(...) — both the receiver and the reference must carry a CRS; two rasters that both report epsg is None are no longer treated as matching.
  • Dataset.crop(bbox=...) / NetCDF.crop(bbox=...) without an explicit epsg=.

A cube with no CRS no longer writes a fabricated one. Previously both writers recorded 4326 and a WGS 84 WKT, claiming a projection the data never had. Now:

  • to_zarr records epsg: 0 with an empty crs_wkt. Readers should treat that pair as "no CRS"; geobox_crs() does.
  • to_netcdf omits both attributes entirely, since a NetCDF has no geobox slot to put a null in.

On read, a NetCDF's crs_wkt / epsg root attributes are adopted only when written beside the GeoTransform the same writer emits, so a stray attribute in a third-party file no longer defines the CRS.

read_part(bbox=...) and point(...) changed their default CRS, for every raster. bbox_crs / point_crs used to default to 4326, so an unqualified bbox or point was interpreted as lon/lat and reprojected into the raster's CRS. They now default to None, meaning "already in the raster's own coordinates", and nothing is transformed. On a raster that is not in EPSG:4326 this changes which pixel an unqualified call reads.

  • To keep the old behaviour, pass the CRS explicitly: ds.read_part(bbox, bbox_crs=4326), ds.point(x, y, point_crs=4326).
  • To read in the raster's own coordinates — usually what you want — leave it out.

Neither method refuses a raster with no CRS, so a windowed read of an ungeoreferenced raster keeps working. Passing an explicit bbox_crs / point_crs against a raster that has none raises CRSError, rather than being ignored: there is no frame to transform into.

feature#

unreleased#

FeatureCollection.from_wfs and from_ogc_features accept an antimeridian bbox. Same change as the raster readers, for the same reason, and with the same caveat about a transposed longitude pair no longer being caught. The vector side has a sharper motivation, though: a wrapping rect handed to OGR whole does not raise. It is silently normalised into its complement, so a filter written as (170, -10, -170, 10) — the 20 degrees around Fiji — quietly returned the other 340 degrees of the planet. The bbox is now split at the seam, each half requested on its own, and the two results merged with duplicates removed.

  • A wrapping bbox costs two requests. A non-wrapping one still costs exactly one; the split is a no-op for it.
  • Duplicates are removed by geometry and attributes, never by FID. A feature that straddles the seam comes back from both halves, and servers assign FIDs per request, so the same feature can arrive under two different ids — keying on the FID would drop real features rather than duplicate ones. Matching on the geometry's WKB plus the attribute values is what actually de-duplicates it. The accepted cost: two rows the source genuinely holds twice, identical in geometry and in every attribute, are indistinguishable and collapse to one. This only applies to a wrapping bbox — a single-request read is passed through untouched.
  • FeatureCollection.fishnet still refuses a wrapping bbox. Its bounds are in the target crs, which may be projected and so has no seam; there is no coherent column numbering across a wrap; and a cell containing ±180 would have to be a MultiPolygon.

pyramids.feature.bbox no longer carries Transformer and crs_from_user_input. Both were incidental re-exports — names the module imported to implement transform, never advertised in its docs or __all__. The reprojection itself moved down to pyramids.base._bbox so pyramids.base could stop importing pyramids.feature (and, with it, geopandas), and its imports went with it. Nothing that the module documents changed: from pyramids.feature.bbox import Bbox, transform still works and returns the same objects as before.

If you imported either name through pyramids.feature.bbox, import it from its own home instead:

from pyproj import Transformer
from pyramids.base.crs import crs_from_user_input

cli#

unreleased#

pyramids calc refuses inputs that do not share the first input's grid. Hard change — it used to let the bound names broadcast against each other and wrote the answer on the first input's grid, a georeferenced result to a question the inputs never agreed on. Warp them onto a common grid first (pyramids warp). An input carrying no CRS tag at all is compared on its pixel grid alone, so an untagged mask or QA layer sitting on the template's cells is still accepted.

pyramids calc preserves a rotated or anisotropic geotransform. Bug fix, silent before. The output was rebuilt from the template's top-left corner and cell size, which collapses to a north-up square-pixel grid: an input at (100, 2, 0.5, 200, 0.25, -3) was written out as (100, 2, 0, 200, 0, -2), with its skew zeroed and its y-resolution changed. The whole geotransform is now copied. The output's no-data value is unchanged — still -9999 whatever the inputs declare, falling back per dtype exactly as before when -9999 is out of range (255 for uint8, 65535 for uint16).

0.47.0#

pyramids calc refuses a first input with no CRS. The result cannot be georeferenced and pyramids will not stamp a default; set a CRS on the input first. pyramids georeference is unaffected — its GCPs and --gcp-crs replace the georeference wholesale.

netcdf#

unreleased#

A container now duck-types as a sequence, because it gained __iter__ and __len__. NetCDF had neither dunder before, so every isinstance(x, Iterable) test answered False for one. That flipped:

expression before now
isinstance(nc, collections.abc.Iterable / Sized / Container) False True
isinstance(nc, collections.abc.Mapping / Sequence) False False, unchanged
  • The case to check is a helper of your own that accepts "anything iterable". It used to reject a container outright and now quietly receives a list of name strings, so a mistake fails later and less clearly than it used to. Gate on isinstance(x, NetCDF) first if that matters.
  • Array coercion now raises, where a container used to answer uselessly. np.asarray(nc) returned a 0-d object array wrapping the container; it now raises TypeError pointing at nc["name"].read_array(). That is a hard change, but the old value was a box around the object rather than any of its data, so nothing could have been computing with it.

The refusal exists for the variable case, where leaving the coercion alone would have been worse than useless: iterating a variable yields nothing, so np.mean(nc["t2m"]) would have answered nan for a cube of real values. origin/main raised there, and so does this — only the exception type differs. - np.array([nc], dtype=object) still boxes the dataset and is still (1,): that request reads nothing and fabricates nothing, so it is the one __array__ honours. - bool(nc) is unaffected: Dataset.__bool__ still refuses, and takes precedence over __len__. Ask len(nc) == 0 whether a container is empty.

to_xarray() decodes the CF time axis, so the coordinate is datetime64[ns] rather than float64. Hard change, silent for the read side — nothing raises and nothing warns when the axis decodes. The bridge exists to hand you to xarray, and resample, .dt and groupby("time.<component>") all raised on the numeric index it used to produce.

xds = nc.to_xarray()
xds.coords["time"].values[:2]   # was [0., 6.]  ->  now ['2024-01-01T00:00', '2024-01-01T06:00']
xds.coords["time"].attrs["units"]      # was 'hours since 2024-01-01'  ->  now KeyError
xds.coords["time"].encoding["units"]   # 'hours since 2024-01-01'
  • units and calendar move from attrs to encoding. The values are no longer expressed in them, so carrying them as attributes would describe the coordinate wrongly and invite a double decode on write. That is where xarray itself keeps them. Code reading xds.time.attrs["units"] must read .encoding["units"].
  • to_xarray(decode_times=False) restores the old output exactly — the stored offsets, with units and calendar still in attrs. It is the spelling that now matches get_dimension_values("time").
  • A CF bounds array follows the coordinate that names it. CF says a bounds variable declares no units of its own and inherits the parent's, so a decoded time no longer sits beside a numeric time_bnds.
  • Only a standard calendar inside datetime64[ns] decodes. A 360_day or noleap axis would decode to cftime objects, which GDAL has no band type for, and an instant outside the type's range would wrap rather than raise. Those axes keep their offsets — and now say so, with a new TimeDecodingWarning naming the dimension and the reason. Five of the repo's own fixtures take that path, so a project running -W error::UserWarning will see it. Silence it with warnings.filterwarnings("ignore", category=TimeDecodingWarning) from pyramids.errors.

from_xarray() writes the time axis in the CF units the array's encoding names, not seconds since 1970-01-01. Hard change, silent — the instants are the same, the file is byte-for-byte different. A datetime64 coordinate carrying units (and optionally calendar) in its encoding is encoded back into them, so a to_xarray() → from_xarray() round trip returns the axis it started with instead of rebasing it on the epoch and stamping a proleptic_gregorian calendar the source never declared. An array with no such encoding still takes the epoch, unchanged.

get_variable returns a LabeledArray, not a raw gdal.MDArray, for a variable with no raster plane. Two shapes are affected: a 1-D array (a profile axis, a bounds array, a hybrid-sigma coefficient), and a string or compound array. GDAL cannot expose either as a raster, and get_variable used to hand the MDArray straight back.

bounds = nc.get_variable("time_bounds")   # was gdal.MDArray, now LabeledArray
bounds.values                             # np.ndarray  (was: bounds.ReadAsArray())
bounds.dims                               # ('number_of_time_bounds',)
bounds.shape                              # (2,)
  • Replace the GDAL calls with attributes. ReadAsArray() → .values, GetDimensions() → .dims / .shape, GetDataType() → .values.dtype. For a string variable Read() was the reader (ReadAsArray cannot serve one), and its replacement is .values — or .values.tolist() for the plain list Read() returned.
  • The labels come with it. GetUnit() → .unit, GetNoDataValueAsDouble() → .no_data_value, GetScale() / GetOffset() → .scale / .offset, and the attribute API → .attributes. Without these a -9999.0 in values would be indistinguishable from real data, and a packed variable's unit would label the wrong numbers. .name is the name you asked for, group path included.
  • Values are as stored. No CF time decoding, and no scale_factor / add_offset applied — matching read_array's own unpack=False default. Unpack with .values * .scale + .offset, or ask read_array(name, unpack=True).
  • A string variable is object, holding str and None for a missing entry. A <U array cannot represent a missing entry, so letting NumPy choose made the dtype depend on whether every record happened to be written.
  • A compound variable decodes to a structured ndarray — index fields with values["code"]; it is not an np.recarray, so values.code does not work. GDAL's netCDF driver reduces a fixed-length character-array member of a compound to its first byte (u1); that is a driver limitation this does not correct.
  • Windowed and strided reads are gone. The MDArray allowed Read(array_start_idx=..., count=...) and GetView; the wrapper is materialised whole, so slice .values. There is no windowed route for these variables now: read_array(variable=...) also reads the whole array, and refuses window, bbox and chunks for them.
  • nc.variables[name].values is read-only for these variables. The mapping caches one object per name and hands it to every caller, so a write would have been visible through variables but not get_variable or read_array. Take .copy() to modify it; get_variable returns a fresh, writeable wrapper each call.
  • Two things now raise that did not. A compound record with a string field cannot be read through GDAL's Python bindings by any means, so get_variable refuses it by name with a ValueError, where it used to return an unreadable handle. And the raster-only operations — crop_variable, reproject_variable, resample_variable, plot(variable=), open_mfdataset(variable=) — refuse a non-raster variable with a ValueError naming it, where they used to fail with AttributeError: 'MDArray' object has no attribute ....
  • LabeledArray is the class LabeledDataset["var"] already returned, exported from pyramids.netcdf. The two agree on type but not always on values: LabeledDataset decodes a CF time axis and applies its store's current selection. The private _read_md_array still returns the MDArray; only the public accessor wraps.

A time axis outside datetime64[ns]'s range now decodes to cftime objects instead of wrapping. Soft change, warned — a UserWarning names the axis and its units. Only arrays that were previously wrong change: decode_cf_time cast to datetime64[ns] under a guard that cannot fire, because a date beyond the type's 1677-09-21 to 2262-04-11 range does not raise on the cast, it wraps.

import numpy as np
from pyramids.netcdf.utils import decode_cf_time

decode_cf_time(np.array([400_000]), "days since 1970-01-01", "standard")
# before -> np.datetime64('1896-01-21T00:50:52.580896768')   # the date is year 3065
# after  -> cftime.real_datetime(3065, 3, 1)                 # + UserWarning

A missing timestep is now NaT rather than the epoch. cftime marks the offsets it cannot decode — a NaN or an inf — and that mark used to be dropped, leaving the fill value, which is the origin, so a missing timestep read back as a real date. It is now NaT in a datetime64 result and None in an object one. Values change wherever an axis has missing offsets, and they change from wrong to right; the integer decode path has always behaved this way, so only the cftime path moves.

A "months since …" axis on a non-360_day calendar now raises a clearer error. It always raised — a calendar month has no fixed length, so cftime allows the unit only on 360_day — but the ValueError came from inside cftime, naming neither the axis nor the store. It now names the axis, units and calendar and chains the original. Nothing that previously succeeded now fails.

Both bounds are affected, and no large offset is needed to reach one: a store written against a days since 0001-01-01 epoch is out of range at offset zero. Its in-range dates are unaffected — a 20th-century date on that epoch still decodes to datetime64[ns] exactly as before.

Which object you get back is cftime's choice, and it decides what still works downstream. It follows the origin, not the dates: when proleptic Gregorian rules already cover the origin — a proleptic_gregorian calendar, or a mixed-calendar origin after the 1582 reform — you get cftime.real_datetime, a datetime subclass, so pandas gives datetime64[us] and to_dataframe / to_parquet / to_csv all keep working and now carry the correct date. That covers the common far-future case above. A pre-1582 origin on a mixed calendar yields a true cftime datetime instead, for every value on the axis — including its post-1582 ones.

LabeledDataset.to_parquet raises on that second case, where it previously wrote a file full of wrapped dates. Parquet has no type for a cftime datetime, so the write now fails with a FailedToSaveError naming the offending columns rather than an ArrowInvalid from inside pyarrow. This was already the behaviour for a non-standard calendar (360_day, noleap), which has always produced cftime objects — the change is that the error explains itself. Use to_csv, which writes these stores unchanged, or select an in-range window first.

If you need datetime64 regardless, the values were never trustworthy in this range; convert deliberately from the returned objects, or read the axis with a calendar-aware library. Everything inside the range is unchanged.

str(nc) is a different shape, and the summary now depends on whether you hold a container or a variable. Hard change, silent — nothing raises and nothing warns. Anything scraping the old text (log parsing, notebook snapshots, doctests) will stop matching.

__str__ was defined once on NetCDF and inherited by Container and Variable alike, so a container printed raster fields it has no raster for: for a 12x5x5 cube at 0.25 degrees it reported Dimension: 512 * 512 at Cell size: 1.0 — GDAL's in-memory placeholder for a multidimensional store. The summary is now chosen by type.

Before, both types printed the same block — Cell size, Dimension, EPSG, projection, Variables, Metadata, File — every line indented by 12 spaces.

After, a container describes the store — this is the real output for tests/data/netcdf/cf__5v__1d4-3d1__geog__y-desc.nc:

<Container cf__5v__1d4-3d1__geog__y-desc.nc>
  dimensions : valid_time=12, latitude=5, longitude=5
  variables  :
    expver      (12)  unknown
    latitude    (5)  float64  degrees_north
    longitude   (5)  float64  degrees_east
    t2m         (12, 5, 5)  float32  K
    valid_time  (12)  int64  seconds since 1970-01-01
  groups     : none
  CRS        : EPSG:4326
  attributes : 6 global

The listing is every array in the store, so the coordinate and auxiliary variables appear alongside the data variables — a wider set than variable_names returns. On a grouped store each row is labelled by its full group/name path, the way variable_names spells its names. At most ten rows are shown, followed by a ... N more line; the dimensions and groups lines truncate the same way once they get long.

and a variable describes the raster it is:

<Variable t2m - cf__5v__1d4-3d1__geog__y-desc.nc>
  grid    : 5 x 5 @ 0.25, EPSG:4326
  bands   : 12 along valid_time
  units   : K
  dtype   : float32
  no-data : nan

Three fields are gone from the text on purpose:

  • the raw projection WKT, which ran to thousands of characters on one line — the CRS now prints as EPSG:4326, or the CRS name when there is no authority code;
  • the whole meta_data object, replaced by a count of global attributes;
  • the leading 12 spaces on every line, which came from the f-string being indented inside the method.

repr() is unchanged, and still shows the placeholder. Only str() / print(nc) route to the new summary. repr() remains gdal.Info on the underlying raster, so a bare nc in a notebook cell, a debugger's variable pane and pytest's assertion output all still report Size is 512, 512 for a container. Use print(nc) where you want the summary.

If you were parsing str(), read the properties instead: nc.cell_size, nc.rows, nc.columns, nc.epsg, nc.crs, nc.variable_names, nc.meta_data.

NetCDF.variable_names answers a different question, so the set of names changes. Hard change, silent — nothing raises and nothing warns, and both the membership and the order can move. The property used to hand back the CF classification's own list; it now filters the store's declared list, which means it walks sub-groups and it drops the arrays CF says are not data. nc.variables has the same keys, so it moves with it.

Two things change, and they pull in opposite directions:

  • Group-qualified names appear. A grouped netCDF-4 store's sub-group arrays are enumerated now, named by their path with / as the separator. On the repository's own none__35v__1d35__groups-nc4.nc fixture the list goes from 1 name to 29, 28 of which carry a /:
# Before
nc.variable_names     # ['UTC_time']

# After
nc.variable_names     # ['UTC_time',
                      #  'mozaic_flight_2012030403540535_ascent/air_press',
                      #  'mozaic_flight_2012030403540535_ascent/CO', ...]   <- 29 in all
  • CF's non-data arrays leave. A bounds array, a 2-D curvilinear coordinate field, an ancillary variable is still in the store and still readable, but it is no longer enumerated as a data variable. Measured over the 26 netCDF fixtures in the repository: 18 lists are byte-identical, 8 changed, and none changed order alone. What left, by fixture: lat_bnds / lon_bnds / time_bnds (two fixtures), lon_rho / lat_rho / Cs_r / h (a curvilinear ROMS grid, 6 names → 2), xc / yc, DQF / band_id / band_wavelength, expver, and UGRID's face_node_connectivity.

Nothing became unreachable. Every name that left is still accepted by get_variable, and asking variables for one says so rather than raising a bare KeyError:

nc.variable_names                 # ['tos']            <- was ['lon_bnds', 'lat_bnds', 'time_bnds', 'tos']
nc.get_variable("lat_bnds")       # still works, shape (170, 2)
nc.variables["lat_bnds"]          # KeyError: "'lat_bnds' is not a data variable, so it is not a key of
                                  #  `variables`; read it with `get_variable('lat_bnds')`."

A group-qualified name works the same way, and a group can also be opened on its own:

nc.get_variable("mozaic_flight_2012030403540535_ascent/CO")        # shape (74,)
nc.get_group("mozaic_flight_2012030403540535_ascent").variable_names
# ['air_press', 'CO', 'O3', 'altitude']

The order within the list is the store's declaration order rather than the classification's. That matters beyond iteration: _fan_out_eager templates a container-wide result from the first spatial variable — taking its geotransform, CRS, no-data and extra dimensions — so a container-wide to_crs, resample or crop can template from a different variable than it did before, and the order propagates into to_netcdf and to_xarray output.

What to do: if you iterated variable_names to convert, export or plot a whole store, re-measure what you get. A grouped store now hands you every group's arrays under /-qualified names — filter on "/" not in name if you want only the root's. A CF store hands you fewer names — name the bounds or ancillary arrays you need explicitly and read them with get_variable. If you depended on the sorted order, sort explicitly: sorted(nc.variable_names).

NetCDF.to_xarray() on a grouped store exports every group's arrays, under flattened names. Hard change, silent. It followed variable_names, so where that returned one name the export had one data variable. On the same grouped fixture:

# Before
list(nc.to_xarray().data_vars)    # ['UTC_time']

# After
list(nc.to_xarray().data_vars)    # ['UTC_time',
                                  #  'mozaic_flight_2012030403540535_ascent_air_press',
                                  #  'mozaic_flight_2012030403540535_ascent_CO', ...]   <- 9 in all

Three rules decide those keys, and they are worth knowing before you index the result:

  • The / becomes _. An xr.Dataset is one flat namespace and netCDF forbids / in a variable name, so a path is flattened with the group kept as a prefix — flight_a/CO and flight_b/CO stay apart as flight_a_CO and flight_b_CO. Exporting the group path verbatim produced a Dataset that could not be written back at all.
  • The store's own name wins a collision. Where a flattened sub-group name lands on a name the store really holds — a root array literally called flight_a_CO, or a dimension coordinate of that name — the root array keeps the plain key and the sub-group array takes _2, _3, …, whichever order the two happen to be enumerated in. A UserWarning lists each rename as store name -> export name; the name before the arrow is the one get_variable takes, and the renamed variable also carries it as a pyramids_store_name attribute, since the suffixed key alone cannot be turned back into it.
  • A variable whose dimensions clash is skipped, with a warning naming it. An xr.Dataset has one size per dimension name, which a grouped store need not honour: the fixture above holds 29 arrays across eight flight groups that all declare a recNum of their own, so 9 are exported and the other 20 are named in a UserWarning telling you to read them with get_variable() instead.

If you keyed the exported Dataset, or round-tripped it through to_netcdf, re-derive the keys rather than hard-coding them — and catch the warnings if you need to know what was renamed or skipped.

A reprojected, resampled or cropped container reports its own grid, not the source's. Hard change, silent. to_crs, resample and crop fan out over a container and build a new one, and the result's x / y axis arrays are what its georeference is read from — pyramids prefers a coordinate pair over the stored transform. Two things about those axes were wrong, in two different arms, and both are fixed.

The in-memory result's axes are float64. They took the data variable's dtype, so reprojecting an int16 raster into EPSG:3857 wrote metre coordinates into an int16 axis, where every one of them saturated. coards__4v__1d3-3d1__y-desc.nc → to_crs(3857):

# Before                                          # After
x  dtype int16,  [-32768, -32768, -32768, …]      x  dtype float64, [-17788309.8, -17464393.7, …]
y  dtype int16,  [ 32767,  32767,  32767, …]      y  dtype float64, [ 13331117.2,  13007201.1, …]
geotransform (-32768.0, 0.0, 0, 32767.0, 0, -0.0) geotransform (-17950267.9, 323916.1, 0,
             <- a zero pixel size                              13493075.2, 0, -323916.1)

The written file no longer carries the source's lat / lon. to_crs, resample and crop with a path= argument stream the result to disk, and that writer copied the source's spatial coordinate arrays into the output beside the result's own x / y. Reopening such a file read the lon/lat pair in preference to the stored transform and reported the source's degrees under the result's CRS. Reprojecting cf__7v__1d3-2d3-3d1__y-asc.nc to EPSG:3857 and reading the file back:

# Before                                          # After
epsg          3857                                epsg          3857
geotransform  (0.0, 2.0, 0, 90.0, 0, -1.0)        geotransform  (-20037477.8, 1054539.6, 0,
              <- degrees, under a metre CRS                      242528680.9, 0, -1054817.1)
arrays        […, 'lat', 'lat_bnds', 'lon',       arrays        […, 'lat_bnds', 'lon_bnds',
               'lon_bnds', 'x', 'y']                             'x', 'y']

A spatial axis describes the grid it was written for, and every operation that reaches here has changed that grid, so the result derives its own x / y and no source axis is carried. Bounds arrays are still carried, verbatim and on a bare dimension — netCDF allows a dimension with no coordinate variable — so a lat_bnds in the output still holds the source's rows and no longer has a lat beside it to be read against. That is deliberate: stale metadata next to an array is the lesser of the two costs, against a file that is silently georeferenced wrong.

If you read a written file's lat / lon arrays, read x / y instead, or take the coordinates from the container's geotransform. If you consumed the carried lat_bnds / lon_bnds, note that they describe the source grid and always did — they are unchanged, only the axis that used to sit beside them is gone. And if you compared a reprojected container's geotransform against a stored value, re-derive it: on a raster whose data dtype could not hold the new coordinates, the old number was not merely imprecise but degenerate.

0.37.0#

Introduces the Container / Variable type split plus a wave of API consolidation. Each change below shows what changed and the before → after so you can update at a glance.

At a glance#

breaking* = only breaks exact-type checks; isinstance(x, NetCDF) still holds.

Change Kind What you do
read_file / get_variable return Container / Variable breaking* use isinstance(x, NetCDF)
subset() returns a NetCDF, not a Dataset breaking* nothing unless you used type(x) is Dataset
CFInfo is frozen breaking dataclasses.replace(cf, ...)
LabeledDataset.read_file(engine=...) validates engine breaking pass a valid engine name
NetCDF(gdal_dataset) direct construction deprecated use read_file / get_variable
get_variable_names() deprecated variable_names property
NetCDF.plot(colour=...) breaking loose colour kwargs + axes=CoordinateSpec(...)
MetaData / DimMetaData deprecated ClassicDimMetadata / ClassicDimensionInfo
kerchunk backend="kerchunk" deprecated backend="legacy"
_LabeledArray / _apply_unpack renamed LabeledArray / apply_unpack

Breaking changes (update required)#

1. Opening a store returns Container; extracting a variable returns Variable. NetCDF is now a base class with two concrete subclasses. Both still pass isinstance(x, NetCDF), so only exact-type checks break.

# Before — everything was a NetCDF
nc  = NetCDF.read_file("cube.nc")     # NetCDF
var = nc.get_variable("t")            # NetCDF
type(nc) is NetCDF                    # True

# After
nc  = NetCDF.read_file("cube.nc")     # Container
var = nc.get_variable("t")            # Variable
isinstance(nc, NetCDF)                # True   <- use this
type(nc) is NetCDF                    # False  <- this no longer holds

2. subset() returns a NetCDF (a Variable), not a plain Dataset.

# Before
ds = nc.subset("t", time=0)           # Dataset

# After
var = nc.subset("t", time=0)          # Variable (a NetCDF); read_array()/crop()/sel() all still work

3. CFInfo is immutable (frozen). In-place assignment now raises dataclasses.FrozenInstanceError.

# Before
cf.<field> = new_value                # mutated in place

# After
import dataclasses
cf = dataclasses.replace(cf, <field>=new_value)

4. LabeledDataset.read_file(engine=...) validates the engine. An unrecognised value now raises ValueError instead of being silently ignored.

# Before
ds = LabeledDataset.read_file(store, engine="zar")    # typo silently ignored

# After — valid: "zarr", "netcdf", "netcdf4", "hdf5", "h5netcdf", or None
ds = LabeledDataset.read_file(store, engine="zarr")

Deprecations (old still works, warns — update when convenient)#

5. Constructing the base NetCDF(...) directly is deprecated. Use the typed entry points.

# Before
nc = NetCDF(gdal_dataset)

# After
nc  = NetCDF.read_file("cube.nc")     # -> Container
var = nc.get_variable("t")            # -> Variable

6. get_variable_names() -> the variable_names property.

# Before
names = nc.get_variable_names()
# After
names = nc.variable_names

7. NetCDF.plot mirrors Dataset.plot — colour= is gone; colour is loose kwargs and coordinates group into axes=CoordinateSpec(...). The colour= parameter and both the ColorOpts dataclass and its British alias ColourOpts have been removed from pyramids.netcdf. Pass the colour knobs (cmap, vmin, vmax, robust, center, extend, levels, norm) as loose keyword arguments, and group the curvilinear-coordinate parameters (coords, x_dim, y_dim) into the new frozen CoordinateSpec dataclass passed as axes=.

# Before
from pyramids.netcdf import ColourOpts
nc.plot("t2m", colour=ColourOpts(cmap="viridis", robust=True), x_dim="rlon", y_dim="rlat")
# After
from pyramids.netcdf import CoordinateSpec
nc.plot("t2m", cmap="viridis", robust=True, axes=CoordinateSpec(x_dim="rlon", y_dim="rlat"))

For colour-bar styling and data-style presets, pass the cleopatra bags exactly as Dataset.plot does — colorbar=ColorBar(...) (or colorbar=False to hide it) and data_style=DataStyle(...).

8. Classic dimension models renamed.

# Before
from pyramids.netcdf.dimensions import MetaData, DimMetaData
# After
from pyramids.netcdf.dimensions import ClassicDimMetadata, ClassicDimensionInfo

9. kerchunk backend="kerchunk" -> backend="legacy".

# Before
nc.to_kerchunk("refs.json", backend="kerchunk")
# After
nc.to_kerchunk("refs.json", backend="legacy")   # or omit backend= for the native default

Renames where the old name still works (no rush)#

10. Promoted internals (underscore aliases kept).

# Before
from pyramids.netcdf.labeled import _LabeledArray
from pyramids.netcdf._lazy import _apply_unpack
# After
from pyramids.netcdf import LabeledArray
from pyramids.netcdf._lazy import apply_unpack

Also: the private modules _kerchunk.py -> _kerchunk_facade.py and _kerchunk_native.py -> _kerchunk_builder.py were renamed. They are internal — use NetCDF.to_kerchunk / NetCDF.combine_kerchunk rather than importing them.

New, opt-in capabilities (not breaking)#

  • Typed dispatch: branch on isinstance(x, Container) vs isinstance(x, Variable) instead of inspecting is_subset / band_count. Import from pyramids.netcdf or pyramids.netcdf.variable.
  • Cloud read tuning: CloudConfig(vsicurl_tuning=True, curl_cache_size=...) enables the fast single-file /vsicurl/ read preset.

Naming note#

The new public types are named Container and Variable. They read cleanly when namespace-qualified (pyramids.netcdf.Container). If you from pyramids.netcdf import Variable in a module that also uses a labeled-array library, the name collides conceptually with that library's own Variable type — prefer the namespace-qualified form, or alias on import:

from pyramids.netcdf import Variable as NcVariable

Migration checklist#

  1. Pin the new pyramids version in your dependencies.
  2. Run your suite under -W error::DeprecationWarning and fix everything it flags (items 5-10 above).
  3. Fix the hard-behavior-change items (1-4) — they do not warn: search for type(x) is NetCDF, type() / Dataset checks on subset() results, CFInfo mutation, and unrecognised LabeledDataset engines.
  4. Done — isinstance(x, NetCDF) keeps working, so most code needs no change.