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:
Or, in code under test:
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 toaz://, which is unchanged and still names Blob. Credentials are the sameAZURE_STORAGE_*family either way. abfss://is accepted too, if you prefer to be explicit about TLS. There is noadls://scheme: the registered Azure Data Lake name isadl://(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, aUserWarningsays which account will actually be read. SetAZURE_STORAGE_ACCOUNT(orAccountName=insideAZURE_STORAGE_CONNECTION_STRING) to match, or use the bareabfs://<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) →
epsgisNone. Operations that genuinely need a CRS raiseCRSErrornaming the fix instead of proceeding. - A CRS with no EPSG authority (geostationary, spherical-earth GRIB) → unchanged by this release. Note that
epsgis not uniformlyNonehere:NetCDFreportsNonefor a geostationary grid (see #706), while a spherical-earth GRIB still resolves throughepsg_from_wkt's4326default. Readcrswhen you need the authoritative answer for such a grid. - A CF NetCDF with
degrees_east/degrees_northaxes and nogrid_mapping→ still4326. 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 acoordinatesattribute (a curvilinear grid's 2-D lat/lon), or carries a conventional coordinate name. A data variable indegrees_east(a wind direction, a solar angle) is not evidence. - A horizontal axis carrying a linear unit (
m,km) or plaindegrees/radvetoes 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'saxis: 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.
- Only a coordinate variable's units are evidence — one that declares
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 andmaintain_alignment=Trueforms),resampleandalignwith itsscaleandoffsetintact. Several of these dropped the packing before: the cutline border-trim silently discarded it (aDatasetlost it outright; aNetCDFvariable kept it only in a Python attribute), and the threeReprojectImagepaths never carried it at all; - for a NetCDF variable,
_scale/_offsetare authoritative for the read when set, and the band's ownGetScale/GetOffsetare 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_variableand the streaming transforms — so none of them can answer in different units for the same band. On aNetCDFvariable the publicscale/offsetproperties answer from the same rule, so asel()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
UserWarningnames 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 —
NaNfor a floating mosaic, and always formethod="min","max"and"sum", which writeFloat64. An integer band has noNaN, so there the marker is instead a value that band can store and that the mosaic's own cells do not use (65535for aUInt16scene holding small numbers,-9999for 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 declaresNaN, and so domethod="min","max"and"sum"whatever the footprints, becauseNaNcosts 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.0and are now the inherited marker (e.g.-9999.0) orNaN. 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 theNaNmarking its uncovered area was rounded to0; those zeros then won everyfminand were added by everysum, so a mosaic of integer tiles that did not tile contiguously was corrupted whateverno_data_valuesaid:mincollapsed to that value everywhere,sumhad it added into every cell, andmaxkept its covered data but lost the gaps. Sources are now warped into theFloat64the reduction writes.- The marker now reaches the pixels it marks. Uncovered pixels used to hold
init(NaNby default) while the mosaic declared something else, so a mosaic inheriting-9999declared-9999over gaps holdingNaNandread_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_valuereports 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.
nis an override, and its default"nan"was passed togdal.BuildVRTas 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-32768beside one declaring-9999read its holes back as-32768. The default now means "no override". If you were relying onn's default to ignoreNaNcells in sources that declare nothing, passn=explicitly. - A numeric
initis a preference, not an instruction. While the marker is being inherited,initis offered to the sentinel search as the preferred candidate and taken only if the mosaic's own cells do not hold it.init=0over 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_valueis validated. A value that names no number used to reachgdal.Translate, which answered "Nodata value was not set to output band" and wrote no marker at all; it now raisesValueError. The string"inherit"is accepted as the default, since that is how the default renders inhelp()and in the API docs.
Migrating:
- If you relied on the
0default, pass it explicitly:merge_rasters(src, dst, no_data_value=0)restores the old declaration, and is worth a second look — it masks every genuine0in 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 hold0rather thaninit. Addinit=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 withinit, so that the marker you asked for actually masks them. Passinit=alongside it to keep the old fill. On a z-order merge these three are the only changes; amin/max/summerge 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.mergeforwards 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:
- Change
except IndexErrortoexcept ValueErrorif 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
ValueErrorwith this same message, from the read layer further down, socount_domain_cells(band=-1)is unaffected. It is simply refused earlier now, before any data is read. - The new
domain_arearefuses 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 thatValueErrorto catch swapped arguments, check the bbox yourself before the call. Latitude is unaffected:miny >= maxyis still an error, because there is no seam in latitude. - Two wrapping bboxes are refused rather than approximated, both on the
from_wms/from_wmtsside: one with a projectedcrs, where there is no 180 degree meridian andminx > maxxreally is inverted; and one with a corner outside-180 .. 180, where "west of the seam" stops meaning anything. from_wms'ssizeis 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 thesizeyou 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 raisesValueError. The cap exists because accepting a wrap made a transposed bbox silently span 359 degrees — at a fine resolution that becomes twoGetMaprequests 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 inbbox, so it is still taken verbatim. If you were relying on a very fineresolution, passsize=instead.from_wmtsmeasures 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_coveragessizes the two halves together when noresolutionis 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 passresolution=. Direct mode has no descriptor to read, so the two halves are independentGetCoveragerequests; with aresolutionthey are snapped onto one lattice, and without one nothing constrains the server to grid them alike. A mismatch is refused withantimeridian halves were produced at different resolutionsrather than stitched into a misaligned raster. - A seam read through
from_wcsorfrom_ogc_coveragesreturns a plainDataset, even when called on a subclass, because those two merge by rebuilding throughDataset.from_array.from_wmsandfrom_wmtsstitch 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 aDatasetsubclass; nothing inpyramidsdoes.
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
AlignmentErrorrather than resampling. The newds.same_grid(other)predicate answers the question in advance, andalign()is the explicit fix. - A scalar operand is accepted on either side:
ds * 2,2 * ds,20 - dsandds >= 5all work. A scalar takes the same route as a raster operand, sods * 2andds * otheragree on band count, dtype and sentinel. That is deliberately not the same asds.apply(lambda v: v * 2), which transforms one band and keeps the source's sentinel — reach forapplywhen you want that contract.booland 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
funccomputed, 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 declaresNaN: a cellfunccomputed asNaN(0/0in a normalised difference) has no value, so it is a gap and the result says so. Passno_data_value=to choose one, orno_data_value=Nonefor no masking. sum(rasters)works:__radd__absorbs the0thatsum()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, sods + 0,0.0 + dsandds * 1are copies too.ds - 0is absorbed as well, on the right only —0 - dsnegates.ds / 1is not: true division widens an integer band tofloat64, so it is no no-op. Absorption is what keeps the two spellings of one commutative expression from disagreeing: routed throughcombine,ds + 0.0would widen anint16band tofloat64where0.0 + dsis byte-identical, andds + 0would drop the band's declared sentinel. Every other scalar computes throughcombine, so1 + dsreturns a raster rather than raising. The copy is a real cost:sum()andmath.prod()each materialise one extra full raster thatfunctools.reduce(operator.add, rasters)does not, which matters near the memory limit.combineis whole-array: both operands are read in full. For rasters near the memory limit useapply(elementwise=True)orread_array(chunks=).- The module-private
_same_gridhelper inpyramids.dataset.datasetmoved toSpatial.same_grid, faced onDataset. It was never public, but anyone importing it directly must switch toa.same_grid(b). <,<=,>and>=between two rasters, or between a raster and a real scalar, return a Byte mask (1/0, and255wherever either operand was no-data).==and!=are not overridden — they stay identity-based, soDatasetremains usable inassert, in sets and as a dict key; usea.combine(b, np.equal)for the mask.math.prod(rasters)folds likesum(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 auint16/uint32/uint64raster now raisesNoDataValueError, where it used to succeed and leave the dtype maximum behind.uint8and 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,filland the statistics treat every cell as data. That is the honest answer for a band with no free value, and it is what auint8raster 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 (afrom_arrayraster, for instance), the method materialises aMEMlevel instead, and that still builds.- anything wrapping a
gdal.Translate(..., format="VRT")orgdal.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 aMEMraster and is unaffected. This one previously produced the damage above verbatim — a stray.ovrand 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 throughcreate_overviews. The target is checked pre-flight, so the call leaves no store behind, and before thecomputeargument, so a call that is alsocompute=Falsereports the refusal rather than theValueError. Withoutoverview_factorsthe refusal does not apply; a description-less VRT then fails where it always did, into_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 oneread_only=Falsecannot reach — holding the source open writable does not help either. Give the VRT its own sidecar withcreate_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=0name 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=-1never meant "the last member" — it happened to work through Python's list indexing while-5raised a bareIndexErrorfrom the same line. Both now raiseFileFormatNotSupportedError, 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_iat 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
epsforfloat16,float32andfloat64alike. Taking each dtype's ownepswould make afloat16band the loosest of the three — one ULP of half precision is9.8e-4, a window of±9.8around a-10000sentinel, 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 pastint64's maximum, so a band whose sentinel is its own dtype maximum would be judged unable to hold it and mask nothing at all — withfillthen 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_axisread the signedgeotransform[1]/geotransform[5], so a positivegeotransform[5](south-up) ascends where the oldget_y_lat_dimension_arrayforced a descending axis (it negated its argument and documented it as a positive pixel height). A caller that passedabs(geotransform[5])must now pass the signed value. This is what fixesDataset.y/Dataset.latreporting coordinates outside a south-up raster's extent, andDataset.bboxis 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.35where the walk produced10.350000000000001). Compare axis arrays withnumpy.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.modecan 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 reportepsg is Noneare no longer treated as matching.Dataset.crop(bbox=...)/NetCDF.crop(bbox=...)without an explicitepsg=.
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_zarrrecordsepsg: 0with an emptycrs_wkt. Readers should treat that pair as "no CRS";geobox_crs()does.to_netcdfomits 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.fishnetstill refuses a wrapping bbox. Its bounds are in the targetcrs, 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 aMultiPolygon.
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:
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 raisesTypeErrorpointing atnc["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'
unitsandcalendarmove fromattrstoencoding. 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 readingxds.time.attrs["units"]must read.encoding["units"].to_xarray(decode_times=False)restores the old output exactly — the stored offsets, withunitsandcalendarstill inattrs. It is the spelling that now matchesget_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
timeno longer sits beside a numerictime_bnds. - Only a standard calendar inside
datetime64[ns]decodes. A360_dayornoleapaxis would decode tocftimeobjects, 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 newTimeDecodingWarningnaming the dimension and the reason. Five of the repo's own fixtures take that path, so a project running-W error::UserWarningwill see it. Silence it withwarnings.filterwarnings("ignore", category=TimeDecodingWarning)frompyramids.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 variableRead()was the reader (ReadAsArraycannot serve one), and its replacement is.values— or.values.tolist()for the plain listRead()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.0invalueswould be indistinguishable from real data, and a packed variable's unit would label the wrong numbers..nameis the name you asked for, group path included. - Values are as stored. No CF time decoding, and no
scale_factor/add_offsetapplied — matchingread_array's ownunpack=Falsedefault. Unpack with.values * .scale + .offset, or askread_array(name, unpack=True). - A string variable is
object, holdingstrandNonefor a missing entry. A<Uarray 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 withvalues["code"]; it is not annp.recarray, sovalues.codedoes 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
MDArrayallowedRead(array_start_idx=..., count=...)andGetView; 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 refuseswindow,bboxandchunksfor them. nc.variables[name].valuesis 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 throughvariablesbut notget_variableorread_array. Take.copy()to modify it;get_variablereturns 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_variablerefuses it by name with aValueError, 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 aValueErrornaming it, where they used to fail withAttributeError: 'MDArray' object has no attribute .... LabeledArrayis the classLabeledDataset["var"]already returned, exported frompyramids.netcdf. The two agree on type but not always on values:LabeledDatasetdecodes a CF time axis and applies its store's current selection. The private_read_md_arraystill returns theMDArray; 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_dataobject, 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 ownnone__35v__1d35__groups-nc4.ncfixture 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'sface_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_. Anxr.Datasetis one flat namespace and netCDF forbids/in a variable name, so a path is flattened with the group kept as a prefix —flight_a/COandflight_b/COstay apart asflight_a_COandflight_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. AUserWarninglists each rename asstore name -> export name; the name before the arrow is the oneget_variabletakes, and the renamed variable also carries it as apyramids_store_nameattribute, 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.Datasethas 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 arecNumof their own, so 9 are exported and the other 20 are named in aUserWarningtelling you to read them withget_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.
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)vsisinstance(x, Variable)instead of inspectingis_subset/band_count. Import frompyramids.netcdforpyramids.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:
Migration checklist#
- Pin the new
pyramidsversion in your dependencies. - Run your suite under
-W error::DeprecationWarningand fix everything it flags (items 5-10 above). - Fix the hard-behavior-change items (1-4) — they do not warn: search for
type(x) is NetCDF,type()/Datasetchecks onsubset()results,CFInfomutation, and unrecognisedLabeledDatasetengines. - Done —
isinstance(x, NetCDF)keeps working, so most code needs no change.