Reading, inspecting & validating

The read side of the COG surface: inspect structure without touching pixels, validate that a file is a valid COG, and read only the pixels you need via overview-decimated partial reads.

• Inspect — cog_info reads only headers/metadata (compression, predictor, blocksize, dtype, CRS/bounds/resolution, the overview pyramid, per-band tags, colour table) and returns a frozen COGInfo. Cheap even for a large remote COG.
• Validate — validate returns a ValidationReport (usable as a bool); Dataset.is_cog is a fast metadata-only probe and Dataset.validate_cog is the authoritative check.
• Partial reads — read_part / preview / point / read_tile request a smaller output size so GDAL serves from the nearest overview, fetching only the relevant byte ranges over /vsicurl/.

Structured inspection

pyramids.dataset.cog.inspect

Structured Cloud Optimized GeoTIFF inspection.

Provides :func:cog_info — a GDAL-only, metadata-only inspection of a raster that answers "what compression / predictor / blocksize / overview pyramid does this COG use?" without reading any pixels. It depends only on the GDAL Python bindings pyramids already uses.

The result is a frozen :class:COGInfo dataclass carrying the band/geo profile plus a per-level :class:OverviewLevel list. Validity is delegated to :func:pyramids.dataset.cog.validate.validate so :attr:COGInfo.is_cog agrees with :meth:pyramids.dataset.engines.cog.COG.validate_cog.

OverviewLevel dataclass

One level of a COG's internal overview pyramid.

Attributes:

Name	Type	Description
index	int	Zero-based overview index (0 is the coarsest-to-finest order GDAL reports, i.e. index 0 is the first/largest overview).
width	int	Overview width in pixels.
height	int	Overview height in pixels.
blocksize	tuple[int, int]	(block_x, block_y) tile size of this overview.
decimation	int	Integer shrink factor relative to full resolution, round(full_width / width) (e.g. 2, 4, 8).

Source code in src/pyramids/dataset/cog/inspect.py

@dataclass(frozen=True)
class OverviewLevel:
    """One level of a COG's internal overview pyramid.

    Attributes:
        index: Zero-based overview index (0 is the coarsest-to-finest order
            GDAL reports, i.e. index 0 is the first/largest overview).
        width: Overview width in pixels.
        height: Overview height in pixels.
        blocksize: ``(block_x, block_y)`` tile size of this overview.
        decimation: Integer shrink factor relative to full resolution,
            ``round(full_width / width)`` (e.g. ``2``, ``4``, ``8``).
    """

    index: int
    width: int
    height: int
    blocksize: tuple[int, int]
    decimation: int

COGInfo dataclass

Structured metadata describing a (Cloud Optimized) GeoTIFF.

Attributes:

Name	Type	Description
is_cog	bool	True iff the file validates as a COG (delegated to :func:pyramids.dataset.cog.validate.validate).
driver	str	GDAL driver short name (e.g. "GTiff").
width	int	Full-resolution width in pixels.
height	int	Full-resolution height in pixels.
band_count	int	Number of raster bands.
dtype	str	GDAL data-type name of band 1 (e.g. "Float32").
crs_epsg	int | None	EPSG code of the CRS, or None when unresolved.
bounds	tuple[float, float, float, float]	(min_x, min_y, max_x, max_y) in the raster CRS.
resolution	tuple[float, float]	(pixel_width, pixel_height) (both positive).
compression	str | None	IMAGE_STRUCTURE compression token, or None.
predictor	str | None	IMAGE_STRUCTURE predictor token, or None.
interleave	str | None	IMAGE_STRUCTURE interleave token, or None.
blocksize	tuple[int, int]	(block_x, block_y) tile size of the full-res image.
overviews	list[OverviewLevel]	Per-level overview metadata, finest index first.
band_tags	dict[int, dict[str, Any]]	Per-band metadata dict keyed by 1-based band index.
colormap	bool	True when band 1 carries a colour table.

Source code in src/pyramids/dataset/cog/inspect.py

@dataclass(frozen=True)
class COGInfo:
    """Structured metadata describing a (Cloud Optimized) GeoTIFF.

    Attributes:
        is_cog: ``True`` iff the file validates as a COG (delegated to
            :func:`pyramids.dataset.cog.validate.validate`).
        driver: GDAL driver short name (e.g. ``"GTiff"``).
        width: Full-resolution width in pixels.
        height: Full-resolution height in pixels.
        band_count: Number of raster bands.
        dtype: GDAL data-type name of band 1 (e.g. ``"Float32"``).
        crs_epsg: EPSG code of the CRS, or ``None`` when unresolved.
        bounds: ``(min_x, min_y, max_x, max_y)`` in the raster CRS.
        resolution: ``(pixel_width, pixel_height)`` (both positive).
        compression: ``IMAGE_STRUCTURE`` compression token, or ``None``.
        predictor: ``IMAGE_STRUCTURE`` predictor token, or ``None``.
        interleave: ``IMAGE_STRUCTURE`` interleave token, or ``None``.
        blocksize: ``(block_x, block_y)`` tile size of the full-res image.
        overviews: Per-level overview metadata, finest index first.
        band_tags: Per-band metadata dict keyed by 1-based band index.
        colormap: ``True`` when band 1 carries a colour table.
    """

    is_cog: bool
    driver: str
    width: int
    height: int
    band_count: int
    dtype: str
    crs_epsg: int | None
    bounds: tuple[float, float, float, float]
    resolution: tuple[float, float]
    compression: str | None
    predictor: str | None
    interleave: str | None
    blocksize: tuple[int, int]
    overviews: list[OverviewLevel] = field(default_factory=list)
    band_tags: dict[int, dict[str, Any]] = field(default_factory=dict)
    colormap: bool = False

    @property
    def overview_count(self) -> int:
        """Number of overview levels present.

        Returns:
            int: ``len(self.overviews)``.
        """
        return len(self.overviews)

cog_info(path, config=None)

Inspect a raster and return its structured COG metadata.

Reads only headers/metadata (no pixels), so it is cheap even for very large or remote (/vsicurl/) COGs. Validity is determined by the same validator that backs :meth:pyramids.dataset.engines.cog.COG.validate_cog.

Parameters:

Name	Type	Description	Default
path	str | Path	Local path or /vsi* path to a raster GDAL can open.	required
config	dict[str, str] | None	GDAL config options applied (via gdal.config_options) while opening. When None and path is remote, the :data:~pyramids.dataset.cog.options.COG_READ_DEFAULTS are applied.	None

Returns:

Name	Type	Description
COGInfo	COGInfo	The structured metadata, including the overview pyramid.

Raises:

Type	Description
FileNotFoundError	When path cannot be opened by GDAL.

Examples:

• Inspect a COG and read its compression and overview pyramid:

  >>> from pyramids.dataset.cog import cog_info  # doctest: +SKIP
  >>> info = cog_info("scene_cog.tif")  # doctest: +SKIP
  >>> info.compression  # doctest: +SKIP
  'DEFLATE'
  >>> [o.decimation for o in info.overviews]  # doctest: +SKIP
  [2, 4, 8]
  

• A plain (non-COG) GeoTIFF reports is_cog=False with no overviews:

  >>> info = cog_info("plain.tif")  # doctest: +SKIP
  >>> info.is_cog, info.overview_count  # doctest: +SKIP
  (False, 0)
  

Source code in src/pyramids/dataset/cog/inspect.py

def cog_info(path: str | Path, config: dict[str, str] | None = None) -> COGInfo:
    """Inspect a raster and return its structured COG metadata.

    Reads only headers/metadata (no pixels), so it is cheap even for very large
    or remote (``/vsicurl/``) COGs. Validity is determined by the same validator
    that backs :meth:`pyramids.dataset.engines.cog.COG.validate_cog`.

    Args:
        path: Local path or ``/vsi*`` path to a raster GDAL can open.
        config: GDAL config options applied (via `gdal.config_options`) while
            opening. When `None` and `path` is remote, the
            :data:`~pyramids.dataset.cog.options.COG_READ_DEFAULTS` are applied.

    Returns:
        COGInfo: The structured metadata, including the overview pyramid.

    Raises:
        FileNotFoundError: When ``path`` cannot be opened by GDAL.

    Examples:
        - Inspect a COG and read its compression and overview pyramid:
            ```python
            >>> from pyramids.dataset.cog import cog_info  # doctest: +SKIP
            >>> info = cog_info("scene_cog.tif")  # doctest: +SKIP
            >>> info.compression  # doctest: +SKIP
            'DEFLATE'
            >>> [o.decimation for o in info.overviews]  # doctest: +SKIP
            [2, 4, 8]

            ```
        - A plain (non-COG) GeoTIFF reports ``is_cog=False`` with no overviews:
            ```python
            >>> info = cog_info("plain.tif")  # doctest: +SKIP
            >>> info.is_cog, info.overview_count  # doctest: +SKIP
            (False, 0)

            ```
    """
    p = str(path)
    cfg = _resolve_read_config(p, config)
    with config_context(cfg):
        info = _cog_info_impl(p)
    return info

_cog_info_impl(p)

Build the :class:COGInfo for p (config context already applied).

Parameters:

Name	Type	Description	Default
p	str	Local path or /vsi* path.	required

Returns:

Name	Type	Description
COGInfo	COGInfo	The structured metadata.

Raises:

Type	Description
FileNotFoundError	When p cannot be opened by GDAL.

Source code in src/pyramids/dataset/cog/inspect.py

def _cog_info_impl(p: str) -> COGInfo:
    """Build the :class:`COGInfo` for ``p`` (config context already applied).

    Args:
        p: Local path or ``/vsi*`` path.

    Returns:
        COGInfo: The structured metadata.

    Raises:
        FileNotFoundError: When ``p`` cannot be opened by GDAL.
    """
    try:
        ds = gdal.Open(p)
    except RuntimeError as exc:
        # With gdal.UseExceptions() a missing/unopenable path raises rather
        # than returning None; surface it as FileNotFoundError for callers.
        raise FileNotFoundError(p) from exc
    if ds is None:
        raise FileNotFoundError(p)

    try:
        band0 = ds.GetRasterBand(1)
        struct = ds.GetMetadata("IMAGE_STRUCTURE")
        block_x, block_y = band0.GetBlockSize()
        width, height = ds.RasterXSize, ds.RasterYSize

        gt = ds.GetGeoTransform()
        # Not the top-left/bottom-right pair: on a south-up raster that pair is
        # already swapped, and this reported `min_y` above `max_y`.
        min_x, min_y, max_x, max_y = GeoTransform(*gt).extent(width, height)

        # `epsg_of_crs` is the package's answer to "what EPSG code does this
        # CRS declare". Reading the authority code directly returned the ESRI
        # number for an ESRI-authority CRS as though it were EPSG (#965), and
        # raised on an OGC:CRS84 raster.
        epsg = epsg_of_crs(ds.GetProjectionRef())

        overviews: list[OverviewLevel] = []
        for i in range(band0.GetOverviewCount()):
            ovr = band0.GetOverview(i)
            obx, oby = ovr.GetBlockSize()
            decimation = round(width / ovr.XSize) if ovr.XSize else 0
            overviews.append(
                OverviewLevel(
                    index=i,
                    width=ovr.XSize,
                    height=ovr.YSize,
                    blocksize=(obx, oby),
                    decimation=decimation,
                )
            )

        band_tags = {
            i: dict(ds.GetRasterBand(i).GetMetadata())
            for i in range(1, ds.RasterCount + 1)
        }
        info = COGInfo(
            is_cog=validate(p).is_valid,
            driver=ds.GetDriver().ShortName,
            width=width,
            height=height,
            band_count=ds.RasterCount,
            dtype=gdal.GetDataTypeName(band0.DataType),
            crs_epsg=epsg,
            bounds=(min_x, min_y, max_x, max_y),
            resolution=(abs(gt[1]), abs(gt[5])),
            compression=struct.get("COMPRESSION"),
            predictor=struct.get("PREDICTOR"),
            interleave=struct.get("INTERLEAVE"),
            blocksize=(block_x, block_y),
            overviews=overviews,
            band_tags=band_tags,
            colormap=band0.GetColorTable() is not None,
        )
    finally:
        ds = None
    return info

Overview-decimated reads

pyramids.dataset.engines.cog.COG

Bases: _Engine['Dataset']

Cloud Optimized GeoTIFF read/write/validate operations for Dataset.

Owns the real implementations of to_cog, is_cog (property), and validate_cog. Dataset exposes a same-named facade for each so ds.to_cog(...) and ds.cog.to_cog(...) are equivalent.

to_cog is the single owner of COG write policy: it applies the house defaults, resolves the dtype-aware predictor and overview resampling, and runs the STATISTICS retry. The :func:pyramids.dataset.cog.write_cog facade is a thin delegator that only normalises its input and forwards overrides here, so both entry points produce identical output for identical input. The categorical-raster resampling guardrail (_warn_if_categorical_with_averaging) lives here too.

Source code in src/pyramids/dataset/engines/cog.py

class COG(_Engine["Dataset"]):
    """Cloud Optimized GeoTIFF read/write/validate operations for `Dataset`.

    Owns the real implementations of `to_cog`, `is_cog` (property),
    and `validate_cog`. `Dataset` exposes a same-named facade for each
    so `ds.to_cog(...)` and `ds.cog.to_cog(...)` are equivalent.

    `to_cog` is the **single owner of COG write policy**: it applies the
    house defaults, resolves the dtype-aware predictor and overview
    resampling, and runs the `STATISTICS` retry. The
    :func:`pyramids.dataset.cog.write_cog` facade is a thin delegator that
    only normalises its input and forwards overrides here, so both entry
    points produce identical output for identical input. The
    categorical-raster resampling guardrail
    (`_warn_if_categorical_with_averaging`) lives here too.
    """

    def to_cog(
        self,
        path: str | Path,
        *,
        compression: str | Compression | None = None,
        overviews: Overviews | None = None,
        tiling: Tiling | None = None,
        bands: BandSelection | None = None,
        tags: Tags | None = None,
        layout: Layout | None = None,
        config: dict[str, str] | None = None,
        extra: Mapping[str, Any] | list[str] | None = None,
    ) -> Path:
        """Save the dataset as a Cloud Optimized GeoTIFF.

        The write options are organised into grouped, validated dataclasses
        (:class:`~pyramids.dataset.cog.Compression`,
        :class:`~pyramids.dataset.cog.Overviews`,
        :class:`~pyramids.dataset.cog.Tiling`,
        :class:`~pyramids.dataset.cog.BandSelection`,
        :class:`~pyramids.dataset.cog.Tags`,
        :class:`~pyramids.dataset.cog.Layout`) — see each class for its fields.
        The whole set is accessible under the ``cog`` namespace, e.g.
        ``from pyramids.dataset import cog`` then ``cog.Compression(...)``.

        Args:
            path: Destination path. Parent directory must exist.
            compression: How pixel bytes are compressed. Either a named profile
                string (`deflate`, `zstd`, `lzw`, `packbits`, `jpeg`, `webp`,
                `lerc`, `lerc_deflate`, `lerc_zstd`, `raw`) or a
                :class:`~pyramids.dataset.cog.Compression`. `None` uses the
                house `DEFLATE` default. The `jpeg`/`webp` **profile strings**
                enforce dtype/band constraints (Byte; 1-3 / 3-4 bands) up front;
                a direct `Compression(compress="JPEG")` is passed to GDAL
                unchecked. The predictor auto-resolves per source dtype (`2`
                integer, `3` float) unless set on `Compression`.
            overviews: The internal overview pyramid, as an
                :class:`~pyramids.dataset.cog.Overviews`. `None` builds the
                default pyramid with dtype-aware resampling (`mode` for
                categorical sources, `average` for continuous).
            tiling: Reprojection / web-tiling layout, as a
                :class:`~pyramids.dataset.cog.Tiling` (`target_srs`, `scheme`
                e.g. `"GoogleMapsCompatible"`, warp `resampling`, zoom knobs).
                `None` writes in the source CRS with no tiling scheme.
            bands: Band selection / dtype cast / NoData applied before the
                write, as a :class:`~pyramids.dataset.cog.BandSelection`. Any of
                its fields pre-processes the source through an in-memory
                `gdal.Translate` so the predictor/overview policy sees the
                output bands.
            tags: Metadata and colour table to stamp onto the output, as a
                :class:`~pyramids.dataset.cog.Tags` (`band_tags`, `colormap`,
                `metadata`).
            layout: Physical layout / driver behaviour, as a
                :class:`~pyramids.dataset.cog.Layout` (`blocksize`, `bigtiff`,
                `num_threads`, `add_mask`, `sparse_ok`, `statistics`). `None`
                uses the house defaults.
            config: GDAL config options (e.g. `{"GDAL_NUM_THREADS": "4"}`)
                applied via `gdal.config_options` for the duration of the
                write. `None` (default) applies no extra config.
            extra: Additional GDAL creation options as a mapping or
                legacy `['KEY=VALUE',...]` list. Overrides conflicting
                group fields.

        Returns:
            Path: The resolved destination path.

        Raises:
            ValueError: Invalid blocksize or unknown option key.
            FileNotFoundError: Parent directory does not exist.
            FailedToSaveError: GDAL CreateCopy failed.
            DriverNotExistError: GDAL build lacks the COG driver.

        Warnings:
            UserWarning: When the source looks categorical (integer dtype or a
                color table) and `overviews.resampling` is an averaging method;
                and when both `tiling.scheme` and `tiling.target_srs` are set
                (`scheme` wins, `target_srs` is ignored).

        Note:
            **Larger-than-RAM / parallel writes.** The GDAL COG driver does the
            two-pass overview layout internally and *streams* from the source
            dataset, so a raster bigger than RAM can be COG-encoded as long as
            the source is **on-disk** (or a `/vsi*` file) rather than a fully
            in-RAM array — anchor a MEM dataset with `to_file(path)` first if
            needed. There is no truly dask-parallel COG writer yet:
            `to_file(compute=False)` returns a `dask.delayed` that wraps the
            *synchronous* GDAL write (GeoTIFF writes are serialised by GDAL's
            own file lock), so it defers *scheduling*, not memory or per-tile
            parallelism. For parallel cloud writes use a Zarr-backed output.

        Examples:
            - Write a compressed COG from an in-memory Dataset:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> out = ds.to_cog("out.tif", compression="zstd")  # doctest: +SKIP
                >>> out.name  # doctest: +SKIP
                'out.tif'

                ```
            - Produce a web-optimized COG for a tile server:
                ```python
                >>> from pyramids.dataset import cog  # doctest: +SKIP
                >>> web = ds.to_cog(  # doctest: +SKIP
                ...     "web.tif", tiling=cog.Tiling(scheme="GoogleMapsCompatible"),
                ... )

                ```
            - Select bands and forward additional GDAL options through `extra`:
                ```python
                >>> from pyramids.dataset import cog  # doctest: +SKIP
                >>> _ = ds.to_cog(  # doctest: +SKIP
                ...     "precise.tif",
                ...     compression="lerc",
                ...     bands=cog.BandSelection(indexes=[2, 1, 0]),
                ...     extra={"MAX_Z_ERROR": 0.001},
                ... )

                ```
        """
        # The jpeg/webp dtype/band check is a named-profile convenience (it mirrors
        # the profile presets), so it fires only when the method was selected via a
        # profile string — a direct `Compression(compress="JPEG")` goes straight to
        # GDAL, matching the pre-refactor flat `compress="JPEG"` path (which GDAL
        # accepts for e.g. 4-band Byte).
        compression_from_profile = isinstance(compression, str)
        compression = Compression.coerce(compression)
        if compression is None:
            compression = Compression()
        overviews = overviews or Overviews()
        tiling = tiling or Tiling()
        bands = bands or BandSelection()
        tags = tags or Tags()
        layout = layout or Layout()

        # Build the effective source (PB-4): band-subset/cast/NoData routes through
        # `BandSelection._translate` and tag/colourmap/metadata through
        # `Tags._stamp`, so the predictor/overview policy below — and the COG write
        # itself — see the *output* bands, and the user's dataset is never mutated.
        source_ds, source_band0 = self._effective_source(bands, tags)

        # Single house policy (ARC-1): each option group serializes its own fields
        # — `Compression`/`Overviews` resolve the dtype-aware predictor and overview
        # resampling from the source band — so a direct `ds.to_cog(...)` and the
        # `write_cog(...)` facade produce identical output for identical input.
        # Assembled before the checks below so `Tiling._to_options`'
        # scheme-vs-target_srs conflict warning keeps firing ahead of them (and is
        # not swallowed by a `validate_profile` raise), matching pre-refactor order.
        defaults = {
            **compression._to_options(source_band0),
            **overviews._to_options(source_band0),
            **tiling._to_options(),
            **layout._to_options(),
        }

        # The jpeg/webp dtype/band constraints (PB-5) mirror the named-profile
        # presets, so they are enforced against the *effective* source only for the
        # profile-string path; a direct `Compression(compress="JPEG")` goes straight
        # to GDAL, matching the pre-refactor flat `compress="JPEG"` behaviour.
        if compression_from_profile and compression.compress is not None:
            validate_profile(
                compression.compress.lower(),
                gdal.GetDataTypeName(source_band0.DataType),
                source_ds.RasterCount,
            )

        # Guardrail (ARC-3): warn only when the *caller* explicitly asked for an
        # averaging resampler on categorical data — never for a default resolved
        # inside `Overviews._to_options`.
        if overviews.resampling is not None:
            self._warn_if_categorical_with_averaging(
                overviews.resampling, band=source_band0
            )

        options = merge_options(defaults, extra)
        # If a caller forced a sub-byte-aligned NBITS through `extra`, the
        # promoted-width predictor computed above would make GDAL reject the
        # write — drop the predictor so the caller's explicit width is honoured.
        _reconcile_predictor_with_nbits(options)
        with config_context(config):
            self._translate_with_statistics_retry(path, options, src=source_ds)
        return Path(path)

    def _effective_source(
        self, bands: BandSelection, tags: Tags
    ) -> tuple[gdal.Dataset, Any]:
        """Return the source dataset (optionally pre-processed) and its band 0.

        Orchestrates the two per-group transforms: a band subset / dtype cast /
        NoData override (PB-4) runs the source through
        :meth:`~pyramids.dataset.cog.BandSelection._translate`; band tags /
        colourmap / metadata (PC-2) are applied by
        :meth:`~pyramids.dataset.cog.Tags._stamp` onto a MEM copy so the user's
        open dataset is **never mutated**. With neither, the backing raster is
        returned unchanged.

        Args:
            bands: The band-selection group deciding the ``gdal.Translate``.
            tags: The tag/colourmap/metadata group to stamp.

        Returns:
            A ``(dataset, band1)`` tuple where ``band1`` is GDAL band 1 of the
            returned dataset (used for predictor/resampling resolution).

        Raises:
            FailedToSaveError: When the stamp-only MEM copy fails.
        """
        # The COG write (and the gdal.Translate pre-process) do tiled reads of the source; a
        # NetCDF multidim view can't be window-read by GDAL >= 3.13, so materialise it first (no-op
        # for an ordinary raster).
        self._ds._materialize_md_view()
        needs_translate = bands._needs_translate()
        needs_stamp = tags._has_any()
        if not needs_translate and not needs_stamp:
            ds = self._ds._raster
            return ds, ds.GetRasterBand(1)

        if needs_translate:
            mem = bands._translate(self._ds._raster)
        else:
            # Stamp-only: copy so the user's dataset is not mutated.
            mem = gdal.GetDriverByName("MEM").CreateCopy("", self._ds._raster)
            if mem is None:
                raise FailedToSaveError(
                    "failed to copy the source dataset for COG metadata stamping"
                )
        if needs_stamp:
            tags._stamp(mem)
        return mem, mem.GetRasterBand(1)

    def to_cog_bytes(self, **kwargs: Any) -> bytes:
        """Encode the dataset as a COG and return the file contents as bytes.

        Writes the COG to an in-memory GDAL ``/vsimem/`` file (no temp file on
        disk), reads the bytes back, and unlinks the virtual file. Useful for
        uploading a COG directly to an object store (S3 / GCS / Azure) without
        touching the local filesystem.

        Args:
            **kwargs: Forwarded verbatim to :meth:`to_cog` (e.g.
                ``compression``, ``layout``, ``bands``, ``extra``). The same
                house defaults and dtype-aware resolution apply.

        Returns:
            bytes: The complete COG file contents.

        Raises:
            FailedToSaveError: GDAL failed to encode the COG.

        Examples:
            - Encode an in-memory Dataset to COG bytes and upload them:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("scene.tif")  # doctest: +SKIP
                >>> blob = ds.to_cog_bytes(compression="zstd")  # doctest: +SKIP
                >>> len(blob) > 0  # doctest: +SKIP
                True
                >>> blob[:2] in (b"II", b"MM")  # TIFF byte-order marker  # doctest: +SKIP
                True

                ```
        """
        vsi_path = mint_vsimem("cog")
        try:
            self.to_cog(vsi_path, **kwargs)
            try:
                data = read_vsi_bytes(vsi_path)
            except FileNotFoundError as exc:
                raise FailedToSaveError(
                    f"could not reopen in-memory COG at {vsi_path}"
                ) from exc
        finally:
            # silent_unlink: when to_cog fails before creating the file, a
            # plain gdal.Unlink raises under gdal.UseExceptions() and masks
            # the original exception. Sweep the PAM sidecar too.
            silent_unlink(vsi_path)
            silent_unlink(f"{vsi_path}.aux.xml")
            # Minting registers the path for the exit sweep; reclaiming it here
            # is what makes that entry dead. Saying so keeps the registry from
            # growing one string per call for the life of the process -- a
            # service encoding a COG per request would otherwise never stop.
            unregister_vsimem(vsi_path)
        return data

    def _translate_with_statistics_retry(
        self,
        path: str | Path,
        options: dict[str, Any],
        src: gdal.Dataset | None = None,
    ) -> None:
        """Write the COG, retrying once without STATISTICS on the known failure.

        Some GDAL builds abort the ``STATISTICS=YES`` sampling pass on float
        on-disk sources with "no valid pixels found in sampling". The COG
        itself is fine without embedded statistics, so on that specific error
        we retry once with ``STATISTICS`` dropped. Lives here (ARC-4) rather
        than in the :func:`write_cog` facade so a direct ``ds.to_cog(...)`` is
        equally robust.

        Args:
            path: Destination file path.
            options: Fully-merged COG creation options.
            src: Source :class:`gdal.Dataset` to encode. Defaults to the
                backing raster; a pre-processed in-memory dataset is passed
                when band-subsetting / casting / setting NoData (PB-4).
        """
        if src is None:
            # COG CreateCopy does tiled reads of the source; a NetCDF multidim view can't be
            # window-read by GDAL >= 3.13, so materialise it first (no-op for an ordinary raster).
            self._ds._materialize_md_view()
        source = self._ds._raster if src is None else src

        def _run(opts: dict[str, Any]) -> None:
            dst: gdal.Dataset | None = None
            try:
                dst = translate_to_cog(source, path, opts)
                dst.FlushCache()
            finally:
                dst = None

        try:
            _run(options)
        except (RuntimeError, FailedToSaveError) as exc:
            # translate_to_cog wraps CreateCopy RuntimeErrors into
            # FailedToSaveError; a deferred STATISTICS failure at FlushCache
            # time surfaces as a raw RuntimeError — catch both.
            statistics_on = str(options.get("STATISTICS", "")).upper() in (
                "YES",
                "TRUE",
            )
            if statistics_on and "valid pixels" in str(exc).lower():
                retry = {k: v for k, v in options.items() if k != "STATISTICS"}
                _run(retry)
            else:
                raise

    @property
    @under_gdal_env
    def is_cog(self) -> bool:
        """`True` iff the backing file on disk is a valid COG.

        `False` for MEM datasets, `/vsimem/` paths, and unsaved
        datasets (empty :attr:`file_name`).

        Examples:
            - Check the backing file of a newly-opened COG:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("scene.tif")  # doctest: +SKIP
                >>> ds.is_cog  # doctest: +SKIP
                True

                ```
            - Plain GeoTIFFs and MEM datasets return False:
                ```python
                >>> plain = Dataset.read_file("plain.tif")  # doctest: +SKIP
                >>> plain.is_cog  # doctest: +SKIP
                False

                ```
            - Use in a conditional pipeline:
                ```python
                >>> if not ds.is_cog:  # doctest: +SKIP
                ...     ds.to_cog("fixed.tif")

                ```
        """
        result: bool
        fn = self._on_disk_path()
        if fn is None:
            result = False
        else:
            result = self._is_cog_cheap(fn)
        return result

    @staticmethod
    def _is_cog_cheap(path: str) -> bool:
        """Fast, metadata-only heuristic for "is this file a COG?" (ARC-7).

        Avoids the full COG validator on every `is_cog` access (which reads the
        whole IFD/offset table — costly over `/vsicurl`). Checks: GTiff driver,
        no external `.ovr` sidecar, internally tiled (square blocks or a single
        tile), and internal overviews present when the image is larger than one
        tile. This can FALSE-POSITIVE on a tiled GeoTIFF that is not laid out in
        strict COG order — use :meth:`validate_cog` for the authoritative check.

        Args:
            path: On-disk or remote `/vsi*` path.

        Returns:
            bool: `True` when the file looks like a COG by the cheap heuristic.
        """
        cfg = _resolve_read_config(path, None)
        with config_context(cfg):
            try:
                ds = gdal.Open(path)
            except RuntimeError:
                return False
            if ds is None:
                return False
            try:
                if ds.GetDriver().ShortName != "GTiff":
                    return False
                files = ds.GetFileList() or []
                if any(str(f).lower().endswith(".ovr") for f in files):
                    return False
                band = ds.GetRasterBand(1)
                block_x, block_y = band.GetBlockSize()
                width, height = ds.RasterXSize, ds.RasterYSize
                single_tile = block_x >= width and block_y >= height
                tiled = block_x == block_y or single_tile
                if not tiled:
                    return False
                needs_overviews = max(width, height) > max(block_x, block_y)
                if needs_overviews and band.GetOverviewCount() == 0:
                    return False
                return True
            finally:
                ds = None

    @under_gdal_env
    def validate_cog(
        self, strict: bool = False, config: dict[str, str] | None = None
    ) -> ValidationReport:
        """Validate the backing file as a COG.

        Args:
            strict: If `True`, warnings are treated as errors.
            config: GDAL config options for the read; defaults to the remote
                read tuning for `/vsicurl` paths (see
                :func:`pyramids.dataset.cog.validate.validate`).

        Returns:
            ValidationReport with errors, warnings, and structural details.

        Raises:
            FileNotFoundError: Dataset has no on-disk backing file
                (MEM-only or `/vsimem/`).

        Examples:
            - Validate and branch on the result:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("scene.tif")  # doctest: +SKIP
                >>> report = ds.validate_cog()  # doctest: +SKIP
                >>> bool(report)  # doctest: +SKIP
                True

                ```
            - Strict mode promotes warnings to errors:
                ```python
                >>> strict = ds.validate_cog(strict=True)  # doctest: +SKIP
                >>> if not strict:  # doctest: +SKIP
                ...     for err in strict.errors: print(err)

                ```
            - Inspect structural details from the report:
                ```python
                >>> report.details.get("blocksize")  # doctest: +SKIP
                [512, 512]

                ```
        """
        fn = self._on_disk_path()
        if fn is None:
            raise FileNotFoundError(
                "Dataset has no on-disk backing file to validate "
                "(is this a MEM or /vsimem/ dataset?)"
            )
        return validate(fn, strict=strict, config=config)

    @under_gdal_env
    def info(self, config: dict[str, str] | None = None) -> COGInfo:
        """Return structured COG metadata for the backing file.

        Reads only headers/metadata (no pixels) and reports compression,
        predictor, blocksize, dtype, CRS/bounds/resolution, the overview
        pyramid, per-band tags, and colour-table presence. See
        :class:`pyramids.dataset.cog.inspect.COGInfo`.

        Args:
            config: GDAL config options for the read; defaults to the remote
                read tuning for `/vsicurl` paths.

        Returns:
            COGInfo: The structured metadata for the on-disk file.

        Raises:
            FileNotFoundError: Dataset has no on-disk backing file
                (MEM-only or `/vsimem/`).

        Examples:
            - Inspect a COG's compression and overview pyramid:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
                >>> info = ds.cog_info()  # doctest: +SKIP
                >>> info.compression  # doctest: +SKIP
                'DEFLATE'
                >>> [o.decimation for o in info.overviews]  # doctest: +SKIP
                [2, 4, 8]

                ```
            - Read the tile size and band count:
                ```python
                >>> info.blocksize  # doctest: +SKIP
                (512, 512)
                >>> info.band_count  # doctest: +SKIP
                1

                ```
        """
        fn = self._on_disk_path()
        if fn is None:
            raise FileNotFoundError(
                "Dataset has no on-disk backing file to inspect "
                "(is this a MEM or /vsimem/ dataset?)"
            )
        return cog_info(fn, config=config)

    def _on_disk_path(self) -> str | None:
        """Return the validatable on-disk path of the backing raster, or None.

        A single predicate shared by :attr:`is_cog`, :meth:`validate_cog`, and
        :meth:`info` (ARC-5) so the definition of "has a real backing file to
        validate/inspect" cannot drift between them.

        Returns:
            str | None: The file path when the dataset is backed by a real
            on-disk (or remote `/vsi*`, but not in-memory `/vsimem/`) file;
            `None` for MEM datasets, `/vsimem/` paths, and unsaved datasets.
        """
        # `_vsimem_path` first, and the driver: a MEM container reports a
        # placeholder `file_name` (a bare driver name such as "netcdf"), which
        # is neither empty nor a /vsimem/ path, so the old test accepted it and
        # `is_cog` went on to validate a file that does not exist.
        fn = getattr(self._ds, "_vsimem_path", "") or self._ds.file_name
        if not fn or fn.startswith("/vsimem/"):
            return None
        # `driver_type` calls `_require_open()`, so it is asked only once the
        # cheap tests have passed and only while the handle is open -- a closed
        # dataset keeps returning None rather than raising out of `is_cog`.
        if self._ds._raster is not None and self._ds.driver_type == "memory":
            return None
        return fn

    @staticmethod
    def _reject_nonpositive_output_sizes(
        dst_width: int | None, dst_height: int | None
    ) -> None:
        """Reject a zero or negative explicit output size before any read.

        Guards `read_part` up front: a non-positive `dst_width`/`dst_height`
        otherwise surfaced as a bare `ZeroDivisionError` -- both when the output
        geotransform divides by the size and, once that transform became
        unconditional, on the plain-array path too.

        Args:
            dst_width: Requested output width in pixels, or `None` for native.
            dst_height: Requested output height in pixels, or `None` for native.

        Raises:
            ValueError: Either size is given (not `None`) and is `<= 0`.
        """
        for name, value in (("dst_width", dst_width), ("dst_height", dst_height)):
            if value is not None and value <= 0:
                raise ValueError(f"{name} must be a positive pixel count, not {value}")

    @overload
    def read_part(
        self,
        bbox: tuple[float, float, float, float],
        *,
        dst_width: int | None = ...,
        dst_height: int | None = ...,
        bbox_crs: int | str | None = ...,
        resampling: str = ...,
        band: int | None = ...,
        return_transform: Literal[False] = ...,
    ) -> np.typing.NDArray: ...

    @overload
    def read_part(
        self,
        bbox: tuple[float, float, float, float],
        *,
        dst_width: int | None = ...,
        dst_height: int | None = ...,
        bbox_crs: int | str | None = ...,
        resampling: str = ...,
        band: int | None = ...,
        return_transform: Literal[True],
    ) -> tuple[np.typing.NDArray, tuple[float, float, float, float, float, float]]: ...

    @overload
    def read_part(
        self,
        bbox: tuple[float, float, float, float],
        *,
        dst_width: int | None = ...,
        dst_height: int | None = ...,
        bbox_crs: int | str | None = ...,
        resampling: str = ...,
        band: int | None = ...,
        return_transform: bool,
    ) -> (
        np.typing.NDArray
        | tuple[np.typing.NDArray, tuple[float, float, float, float, float, float]]
    ): ...

    @under_gdal_env
    def read_part(
        self,
        bbox: tuple[float, float, float, float],
        *,
        dst_width: int | None = None,
        dst_height: int | None = None,
        bbox_crs: int | str | None = None,
        resampling: str = "bilinear",
        band: int | None = None,
        return_transform: bool = False,
    ) -> (
        np.typing.NDArray
        | tuple[np.typing.NDArray, tuple[float, float, float, float, float, float]]
    ):
        """Read a geographic window, decimated from the nearest overview.

        Requesting a `dst_width`/`dst_height` smaller than the source window
        makes GDAL serve the data from the nearest overview level, so for a COG
        over `/vsicurl/` only the relevant byte ranges are fetched — the
        cloud-native partial-read pattern.

        The window answers in physical units, like `read_array`: a CF-packed band
        (`scale_factor` / `add_offset`) is unpacked and comes back `float64`, whether
        the window lies inside the raster or only partly overlaps it (see the note
        below).

        Args:
            bbox: `(min_x, min_y, max_x, max_y)` window in `bbox_crs`.
            dst_width: Output width in pixels. Defaults to the source window
                width (no decimation).
            dst_height: Output height in pixels. Defaults to the source window
                height.
            bbox_crs: CRS of `bbox`, reprojected to the dataset CRS when it
                differs. Defaults to `None`, meaning the bbox is already in the
                raster's own coordinates, so nothing is transformed.
            resampling: Resampling method, case-insensitive. One of `nearest`,
                `bilinear`, `cubic`, `cubicspline` (alias `cubic_spline`),
                `lanczos`, `average`, `mode`, plus `gauss` and `rms` when the
                GDAL build provides them.
            band: 0-based band index. `None` reads all bands.
            return_transform: When `True`, return `(array, geotransform)`
                instead of the bare array, so a caller can place the result
                without reproducing the pixel-snapping arithmetic. The
                geotransform is the six-tuple `GetGeoTransform` uses, describing
                the returned buffer at its output resolution: it maps output
                pixel coordinates to the **dataset's** CRS, over the window
                actually read (the requested bbox snapped outward to whole
                source pixels). Its origin, cell size and orientation place the
                buffer's outer extent exactly, and that extent always matches the
                snapped window. Per-cell placement -- each cell's data sitting on
                its own transform centre -- is exact only for a read at native
                resolution (`dst_width`/`dst_height` omitted, so the output keeps
                the snapped window's size). A resampled read can shift a cell's
                data off its transform centre: by a fraction of a cell in the
                interior, and by more near the raster edge, where a window that
                overruns the edge has its out-of-raster remainder padded before
                resampling -- reaching several output cells when such a window is
                upsampled. The outer extent stays exact in every case; read at
                native resolution when you need cell-exact data. Defaults to
                `False`, which keeps the bare-array return unchanged.

        Returns:
            numpy.ndarray: `(rows, cols)` for a single band, or
            `(bands, rows, cols)` for all bands; always sized
            `dst_height x dst_width` (the requested output size). Pixel values
            only — no transform, bounds, or CRS is attached, unless
            `return_transform=True`, in which case a
            `(array, geotransform)` tuple is returned instead (see
            `return_transform`).

        Raises:
            CRSError: An explicit `bbox_crs` was given but the raster has no CRS
                to transform into. Omit it to read in the raster's own
                coordinates (ARC-26).
            TypeError: `resampling` is not a string.
            ValueError: Unknown `resampling`, or `dst_width` / `dst_height` is
                given as a non-positive pixel count.
            OutOfBoundsError: The window does not intersect the raster at all.

        Note:
            A window that only **partially** overlaps the raster is **not**
            stretched to fill the output: the intersection is read and placed
            at its correct offset inside a `dst_height x dst_width` buffer
            whose out-of-raster remainder is filled with NoData (the band's
            NoData value, else NaN for float / `0` for integer — see
            :meth:`_nodata_fill`). A fully-inside window is returned without
            padding. This keeps the result aligned to the requested window,
            which matters for edge tiles served by :meth:`read_tile`.

            The padded partial-overlap read is assembled in stored units --
            the stored values, padded with the stored NoData -- and then
            unpacked once, like a fully-inside read. So on a CF-packed band it
            is in the same physical units and dtype as any other window, and
            on a band that declares NoData the padding reads as its gap does in
            a default read (the sentinel unpacked with the data, e.g. `-9999`
            at `scale=0.5` as `-4999.5`).

        Examples:
            - Read a 256x256 decimated thumbnail of a bbox:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
                >>> arr = ds.read_part(  # doctest: +SKIP
                ...     (12.4, 41.8, 12.6, 42.0), dst_width=256, dst_height=256,
                ... )
                >>> arr.shape[-2:]  # doctest: +SKIP
                (256, 256)

                ```
            - A window inside a CF-packed raster comes back in physical units:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset, GeoReference
                >>> packed = Dataset.from_array(
                ...     np.full((4, 4), 100, dtype="int16"),
                ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
                ... )
                >>> packed.scale = [0.01]
                >>> part = packed.read_part((0, -2, 2, 0), band=0)
                >>> part.shape, part.dtype, float(part.max())
                ((2, 2), dtype('float64'), 1.0)

                ```
            - A window straddling the raster's edge is padded with the NoData and
              unpacked with the data, so it is in the same units as one inside:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset, GeoReference
                >>> packed = Dataset.from_array(
                ...     np.full((4, 4), 100, dtype="int16"),
                ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
                ...     no_data_value=-9999,
                ... )
                >>> packed.scale = [0.5]
                >>> edge = packed.read_part((-1, -2, 1, 0), band=0)
                >>> edge.dtype, edge.tolist()
                (dtype('float64'), [[-4999.5, 50.0], [-4999.5, 50.0]])

                ```
            - Ask for the geotransform of the window actually read. The bbox is
              snapped outward to whole source pixels, so the transform, not the
              requested bbox, is what places the cells:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset, GeoReference
                >>> grid = Dataset.from_array(
                ...     np.tile(np.arange(8, dtype="float64"), (8, 1)),
                ...     geo_ref=GeoReference(top_left_corner=(0, 8), cell_size=1.0, epsg=3857),
                ... )
                >>> arr, gt = grid.read_part(
                ...     (1.5, 1.5, 4.5, 4.5), dst_width=3, dst_height=3, return_transform=True,
                ... )
                >>> arr.shape
                (3, 3)
                >>> gt  # origin snapped to x=1, cells 4/3 wide
                (1.0, 1.3333333333333333, 0.0, 5.0, 0.0, -1.3333333333333333)

                ```
        """
        alg = _resolve_read_resampling(resampling)
        self._reject_nonpositive_output_sizes(dst_width, dst_height)
        # This serves a decimated window from the source; a NetCDF multidim view can't be window-read
        # by GDAL >= 3.13, so materialise it first (no-op for an ordinary raster).
        self._ds._materialize_md_view()
        ds = self._ds._raster
        min_x, min_y, max_x, max_y = self._reproject_bbox(bbox, bbox_crs)
        geotransform = ds.GetGeoTransform()
        px_tl, py_tl = world_to_pixel(geotransform, min_x, max_y)
        px_br, py_br = world_to_pixel(geotransform, max_x, min_y)

        # The full requested window, in source pixel coordinates (may extend
        # beyond the raster on any side).
        req_xoff = int(math.floor(min(px_tl, px_br)))
        req_yoff = int(math.floor(min(py_tl, py_br)))
        req_xsize = int(math.ceil(max(px_tl, px_br))) - req_xoff
        req_ysize = int(math.ceil(max(py_tl, py_br))) - req_yoff
        if req_xsize <= 0 or req_ysize <= 0:
            raise OutOfBoundsError(
                f"bbox {bbox} (crs {bbox_crs}) has zero pixel extent"
            )

        # Intersection of the requested window with the raster.
        ix0 = max(0, req_xoff)
        iy0 = max(0, req_yoff)
        ix1 = min(ds.RasterXSize, req_xoff + req_xsize)
        iy1 = min(ds.RasterYSize, req_yoff + req_ysize)
        if ix1 - ix0 <= 0 or iy1 - iy0 <= 0:
            raise OutOfBoundsError(
                f"bbox {bbox} (crs {bbox_crs}) does not intersect the raster"
            )

        out_w = dst_width if dst_width is not None else req_xsize
        out_h = dst_height if dst_height is not None else req_ysize
        # Describes the buffer both arms below return: the snapped window at the
        # output resolution, in the dataset's CRS. Computed here, before the read,
        # because it depends only on the window and the output size, not the pixels.
        out_gt = self._output_geotransform(
            geotransform, req_xoff, req_yoff, req_xsize, req_ysize, out_w, out_h
        )
        source = ds if band is None else ds.GetRasterBand(band + 1)

        fully_inside = (ix0, iy0, ix1, iy1) == (
            req_xoff,
            req_yoff,
            req_xoff + req_xsize,
            req_yoff + req_ysize,
        )
        if fully_inside:
            # Physical units, like `read_array`. These are overview-decimated reads
            # of the same band, and a raster that renders at one scale through
            # `plot()` and another through `plot(overview=True)` is a defect a
            # reader can see.
            array = np.asarray(
                self._ds.io._apply_scale_offset(
                    np.asarray(
                        source.ReadAsArray(
                            ix0,
                            iy0,
                            ix1 - ix0,
                            iy1 - iy0,
                            buf_xsize=out_w,
                            buf_ysize=out_h,
                            resample_alg=alg,
                        )
                    ),
                    band,
                )
            )
            return (array, out_gt) if return_transform else array

        # Partial overlap: read only the intersection, then place it at its
        # correct offset inside a full-size output buffer padded with NoData,
        # so the returned array stays aligned to the requested window.
        scale_x = out_w / req_xsize
        scale_y = out_h / req_ysize
        ox0 = max(0, min(out_w, int(round((ix0 - req_xoff) * scale_x))))
        oy0 = max(0, min(out_h, int(round((iy0 - req_yoff) * scale_y))))
        ox1 = max(ox0 + 1, min(out_w, int(round((ix1 - req_xoff) * scale_x))))
        oy1 = max(oy0 + 1, min(out_h, int(round((iy1 - req_yoff) * scale_y))))
        sub = np.asarray(
            source.ReadAsArray(
                ix0,
                iy0,
                ix1 - ix0,
                iy1 - iy0,
                buf_xsize=ox1 - ox0,
                buf_ysize=oy1 - oy0,
                resample_alg=alg,
            )
        )
        fill = self._nodata_fill(ds.GetRasterBand(1))
        if sub.ndim == 3:
            out = np.full((sub.shape[0], out_h, out_w), fill, dtype=sub.dtype)
            out[:, oy0:oy1, ox0:ox1] = sub
        else:
            out = np.full((out_h, out_w), fill, dtype=sub.dtype)
            out[oy0:oy1, ox0:ox1] = sub
        # Assembled in stored units -- the counts and the stored sentinel padding them
        # -- and unpacked once, so the padding reads as a gap exactly the way the
        # fully-inside arm above reads its own. Returning the stored array made a window
        # straddling the raster's edge come back as `int16` counts while a window just
        # inside it came back as physical `float64`, and every edge `read_tile` took
        # this arm.
        array = np.asarray(self._ds.io._apply_scale_offset(out, band))
        return (array, out_gt) if return_transform else array

    @staticmethod
    def _output_geotransform(
        source_gt: tuple[float, float, float, float, float, float],
        req_xoff: int,
        req_yoff: int,
        req_xsize: int,
        req_ysize: int,
        out_w: int,
        out_h: int,
    ) -> tuple[float, float, float, float, float, float]:
        """Geotransform of a `read_part` buffer, from the source grid and the window.

        The buffer covers source pixels `[req_xoff, req_xoff + req_xsize)` by
        `[req_yoff, req_yoff + req_ysize)` -- the requested bbox snapped outward
        to whole pixels -- resampled to `out_w` by `out_h`. Composing the source
        affine with that offset-and-scale gives an affine in the dataset's CRS
        that maps an output pixel to the same world point the source grid does,
        so it carries a rotated or south-up source through unchanged (the pixel
        step keeps its sign) rather than assuming a north-up grid.

        Args:
            source_gt: The source `GetGeoTransform` six-tuple.
            req_xoff: Left edge of the window, in source pixels.
            req_yoff: Top edge of the window, in source pixels.
            req_xsize: Window width, in source pixels.
            req_ysize: Window height, in source pixels.
            out_w: Output width, in pixels.
            out_h: Output height, in pixels.

        Returns:
            tuple: The `(origin_x, x_size, x_skew, origin_y, y_skew, y_size)`
            geotransform of the output buffer.

        Examples:
            - A north-up window snapped to whole pixels and decimated 4 source
              pixels into 3 output cells (the same window as `read_part`'s
              transform example): the origin moves to the snapped corner and the
              cell size grows to 4/3:
                ```python
                >>> from pyramids.dataset.engines.cog import COG
                >>> COG._output_geotransform((0.0, 1.0, 0.0, 8.0, 0.0, -1.0), 1, 3, 4, 4, 3, 3)
                (1.0, 1.3333333333333333, 0.0, 5.0, 0.0, -1.3333333333333333)

                ```
            - A south-up source (a positive y-step) keeps its sign instead of
              being forced north-up:
                ```python
                >>> from pyramids.dataset.engines.cog import COG
                >>> gt = COG._output_geotransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0), 0, 0, 8, 8, 4, 4)
                >>> gt[5]  # y-step stays positive
                2.0
                >>> gt
                (0.0, 2.0, 0.0, 0.0, 0.0, 2.0)

                ```
            - A rotated (skewed) source carries its skew terms through the
              compose, so the output stays aligned to the same rotated grid:
                ```python
                >>> from pyramids.dataset.engines.cog import COG
                >>> COG._output_geotransform((0.0, 2.0, 1.0, 0.0, 1.0, -2.0), 1, 1, 4, 4, 2, 2)
                (3.0, 4.0, 2.0, -1.0, 2.0, -4.0)

                ```
        """
        gt0, gt1, gt2, gt3, gt4, gt5 = source_gt
        return (
            gt0 + req_xoff * gt1 + req_yoff * gt2,
            gt1 * req_xsize / out_w,
            gt2 * req_ysize / out_h,
            gt3 + req_xoff * gt4 + req_yoff * gt5,
            gt4 * req_xsize / out_w,
            gt5 * req_ysize / out_h,
        )

    @staticmethod
    def _nodata_fill(band: Any) -> float:
        """Pick a fill value for padding partial reads.

        Args:
            band: The GDAL band whose NoData value (if any) to use.

        Returns:
            float: The band's NoData value, else NaN for floating-point bands
            and ``0`` for integer bands.
        """
        nodata = band.GetNoDataValue()
        if nodata is not None:
            return cast(float, nodata)
        return 0 if is_integer_gdal_dtype(band.DataType) else float("nan")

    @under_gdal_env
    def preview(
        self,
        *,
        max_size: int = 1024,
        resampling: str = "bilinear",
        band: int | None = None,
    ) -> np.typing.NDArray:
        """Read a whole-image thumbnail downsampled to `max_size` on the long edge.

        Pulls from a coarse overview when one exists, so previewing a huge COG
        is cheap. The thumbnail is the same band at a coarser sampling, so it
        answers in the same units as `read_array`: a CF-packed band
        (`scale_factor` / `add_offset`) is unpacked to physical values.

        Args:
            max_size: Maximum pixels on the longer edge. Defaults to 1024.
            resampling: Resampling method (see :meth:`read_part`).
            band: 0-based band index. `None` reads all bands.

        Returns:
            numpy.ndarray: The downsampled array, `(rows, cols)` or
            `(bands, rows, cols)`, `float64` when a band that was read is packed. Pixel
            values only — no transform, bounds, or CRS is attached to the
            returned array.

        Raises:
            TypeError: `resampling` is not a string.
            ValueError: Unknown `resampling`.

        Examples:
            - Build a 128px thumbnail of a single band:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
                >>> thumb = ds.preview(max_size=128, band=0)  # doctest: +SKIP
                >>> max(thumb.shape)  # doctest: +SKIP
                128

                ```
            - A thumbnail of a CF-packed raster is in physical units:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset, GeoReference
                >>> packed = Dataset.from_array(
                ...     np.full((4, 4), 100, dtype="int16"),
                ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
                ... )
                >>> packed.scale = [0.01]
                >>> thumb = packed.preview(max_size=2, band=0)
                >>> thumb.shape, float(thumb[0, 0])
                ((2, 2), 1.0)

                ```
        """
        alg = _resolve_read_resampling(resampling)
        width, height = self._ds.columns, self._ds.rows
        scale = max(width, height) / max_size
        if scale <= 1:
            out_w, out_h = width, height
        else:
            out_w, out_h = max(1, round(width / scale)), max(1, round(height / scale))
        ds = self._ds._raster
        source = ds if band is None else ds.GetRasterBand(band + 1)
        # Physical units, like `read_array`: `preview` is the same band at a coarser
        # sampling, not a different quantity.
        return np.asarray(
            self._ds.io._apply_scale_offset(
                np.asarray(
                    source.ReadAsArray(
                        buf_xsize=out_w, buf_ysize=out_h, resample_alg=alg
                    )
                ),
                band,
            )
        )

    @under_gdal_env
    def point(
        self,
        x: float,
        y: float,
        *,
        point_crs: int | str | None = None,
        band: int | None = None,
    ) -> np.typing.NDArray:
        """Sample band value(s) at a single coordinate.

        The sample answers in the same units as `read_array`: a CF-packed band
        (`scale_factor` / `add_offset`) is unpacked, so `point` and
        `read_array()[row, col]` agree. A no-data cell answers its stored
        sentinel transformed like any other value; compare against
        `no_data_value` on an `unpack=False` read instead.

        Args:
            x: X / longitude / easting in `point_crs`.
            y: Y / latitude / northing in `point_crs`.
            point_crs: CRS of `(x, y)`, reprojected to the dataset CRS when it
                differs. Defaults to `None`, meaning the coordinates are already
                in the raster's own CRS, so nothing is transformed.
            band: 0-based band index. `None` samples all bands.

        Returns:
            numpy.ndarray: A scalar 0-d array for a single band, or a
            `(bands,)` array when `band` is `None`, `float64` when a band that was read
            is packed. Pixel values only — no coordinate metadata is attached.

        Raises:
            CRSError: An explicit `point_crs` was given but the raster has no
                CRS to transform into. Omit it to read in the raster's own
                coordinates (ARC-26).
            OutOfBoundsError: The point falls outside the raster extent.

        Examples:
            - Sample all bands at a lon/lat coordinate:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
                >>> ds.point(12.5, 41.9)  # doctest: +SKIP
                array([1234.], dtype=float32)

                ```
            - On a CF-packed raster the sample matches the default `read_array`:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset, GeoReference
                >>> packed = Dataset.from_array(
                ...     np.array([[100, 200], [300, 400]], dtype="int16"),
                ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
                ... )
                >>> packed.scale = [0.01]
                >>> float(packed.point(1.5, -0.5, band=0)), float(packed.read_array()[0, 1])
                (2.0, 2.0)

                ```
        """
        col, row = self._world_to_pixel(x, y, point_crs)
        if not (0 <= col < self._ds.columns and 0 <= row < self._ds.rows):
            raise OutOfBoundsError(
                f"point ({x}, {y}) in crs {point_crs} is outside the raster extent"
            )
        ds = self._ds._raster
        source = ds if band is None else ds.GetRasterBand(band + 1)
        # Unpacked while the read still has its `(bands, 1, 1)` shape, so the
        # per-band broadcast works; the reshape to a plain value comes after. A
        # `point` that answered 100 where `read_array()[row, col]` answered 2.5 was
        # the same raster sampled two ways.
        arr = np.asarray(
            self._ds.io._apply_scale_offset(
                np.asarray(source.ReadAsArray(col, row, 1, 1)), band
            )
        )
        return arr.reshape(-1) if band is None else arr.reshape(())

    @under_gdal_env
    def read_tile(
        self,
        z: int,
        x: int,
        y: int,
        *,
        tilesize: int = 256,
        resampling: str = "bilinear",
        band: int | None = None,
    ) -> np.typing.NDArray:
        """Read a Web-Mercator XYZ/slippy-map tile.

        Computes the EPSG:3857 bounds of tile `(z, x, y)` from the closed-form
        Web-Mercator formula and delegates to :meth:`read_part` at `tilesize`
        resolution — no extra tiling dependency needed.

        Args:
            z: Zoom level.
            x: Tile column index.
            y: Tile row index (origin top-left / north-west).
            tilesize: Output tile size in pixels (square). Defaults to 256.
            resampling: Resampling method (see :meth:`read_part`).
            band: 0-based band index. `None` reads all bands.

        Returns:
            numpy.ndarray: A `(tilesize, tilesize)` or
            `(bands, tilesize, tilesize)` array, in physical units like
            `read_array` (`float64` when a band read is CF-packed). Pixel values
            only — the tile's georeferencing is defined by its `(z, x, y)`, not
            attached to the array; edge tiles are NoData-padded and unpacked in
            the same units as inner ones (see :meth:`read_part`).

        Raises:
            OutOfBoundsError: The tile does not intersect the raster.

        Examples:
            - Read the zoom-0 world tile of a global COG:
                ```python
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> ds = Dataset.read_file("global_cog.tif")  # doctest: +SKIP
                >>> tile = ds.read_tile(0, 0, 0)  # doctest: +SKIP
                >>> tile.shape[-2:]  # doctest: +SKIP
                (256, 256)

                ```
        """
        bounds = _xyz_bounds_3857(z, x, y)
        return self.read_part(
            bounds,
            dst_width=tilesize,
            dst_height=tilesize,
            bbox_crs=3857,
            resampling=resampling,
            band=band,
        )

    def _reproject_bbox(
        self, bbox: tuple[float, float, float, float], bbox_crs: int | str | None
    ) -> tuple[float, float, float, float]:
        """Reproject a bbox into the dataset CRS, returning its envelope.

        Args:
            bbox: `(min_x, min_y, max_x, max_y)` in `bbox_crs`.
            bbox_crs: CRS of `bbox`. `None` (default) means the bbox is already
                in the raster's own coordinates.

        Returns:
            `(min_x, min_y, max_x, max_y)` in the dataset CRS. When
            `bbox_crs` already matches the dataset EPSG the bbox is
            returned unchanged.
        """
        min_x, min_y, max_x, max_y = bbox
        envelope = bbox
        # `None` means the caller named no CRS, so the bbox is already in the
        # raster's own coordinates and there is nothing to transform (ARC-26).
        # An *explicit* bbox_crs goes through `require_crs_spec`, so a raster
        # with no CRS reports the mismatch rather than silently ignoring the
        # argument.
        if bbox_crs is not None:
            target = require_crs_spec(
                self._ds.epsg, self._ds.crs, "read a bbox window in another CRS"
            )
            if not crs_equal(bbox_crs, target):
                transformer = _cached_transformer(bbox_crs, target)
                corners = [
                    transformer.transform(min_x, min_y),
                    transformer.transform(min_x, max_y),
                    transformer.transform(max_x, min_y),
                    transformer.transform(max_x, max_y),
                ]
                xs = [c[0] for c in corners]
                ys = [c[1] for c in corners]
                envelope = (min(xs), min(ys), max(xs), max(ys))
        return envelope

    def _world_to_pixel(
        self, x: float, y: float, point_crs: int | str | None
    ) -> tuple[int, int]:
        """Convert a world coordinate to integer `(col, row)` pixel indices.

        Args:
            x: X / longitude in `point_crs`.
            y: Y / latitude in `point_crs`.
            point_crs: CRS of `(x, y)`. `None` (the default) means the point is
                already in the raster's own coordinates.

        Returns:
            `(col, row)` integer pixel indices (floored).
        """
        # Mirrors the bbox path: `None` means the point is already in the
        # raster's own coordinates, so there is nothing to transform. Without
        # this early-out the default builds a transformer FROM None and fails on
        # every georeferenced raster. Past it the caller HAS named a CRS, so the
        # raster must have one to transform into -- `require_crs_spec` rather
        # than a silent skip, which ignored the argument that was passed.
        if point_crs is not None:
            target = require_crs_spec(
                self._ds.epsg, self._ds.crs, "sample a point given in another CRS"
            )
            if not crs_equal(point_crs, target):
                transformer = _cached_transformer(point_crs, target)
                x, y = transformer.transform(x, y)
        col, row = world_to_pixel(self._ds._raster.GetGeoTransform(), x, y)
        return int(math.floor(col)), int(math.floor(row))

    def _warn_if_categorical_with_averaging(
        self, overview_resampling: str, band: Any | None = None
    ) -> None:
        """Emit a `UserWarning` if an averaging resampler is used on categorical data.

        Args:
            overview_resampling: The resampling method requested by the
                caller. Case-insensitive. Only averaging-family methods
                (`average`, `bilinear`, `cubic`, `cubicspline`,
                `lanczos`) trigger the check.
            band: GDAL band whose dtype/colour-table decides "categorical".
                Defaults to band 1 of the backing raster; a pre-processed
                (cast/subset) band is passed when those options are used so
                the check reflects the *output* dtype (PB-4).

        Warns:
            UserWarning: When `overview_resampling` is an averaging
                method and the source has a color table OR integer
                dtype — both strong signals of categorical data.

        Note:
            Silent when `overview_resampling` is `nearest` or
            `mode` (both category-safe) or when the source is
            floating-point and has no color table (continuous data).

        Examples:
            - Integer dataset + averaging method emits a warning:
                ```python
                >>> import warnings  # doctest: +SKIP
                >>> with warnings.catch_warnings(record=True) as caught:  # doctest: +SKIP
                ...     warnings.simplefilter("always")
                ...     byte_ds.cog._warn_if_categorical_with_averaging("average")
                ...     [str(w.message) for w in caught if issubclass(w.category, UserWarning)]
                ['overview_resampling=\\'average\\' averages pixel values, ...']

                ```
            - Nearest resampling is always silent:
                ```python
                >>> with warnings.catch_warnings(record=True) as caught:  # doctest: +SKIP
                ...     warnings.simplefilter("always")
                ...     byte_ds.cog._warn_if_categorical_with_averaging("nearest")
                ...     len(caught)
                0

                ```
        """
        if overview_resampling.lower() not in _AVERAGING_RESAMPLERS:
            return
        first_band = band if band is not None else self._ds._raster.GetRasterBand(1)
        has_color_table = first_band.GetColorTable() is not None
        is_integer = is_integer_gdal_dtype(first_band.DataType)
        if has_color_table or is_integer:
            warnings.warn(
                f"overview_resampling={overview_resampling!r} averages pixel "
                "values, which corrupts categorical rasters (land cover, IDs). "
                "Use overview_resampling='nearest' or 'mode' instead.",
                UserWarning,
                stacklevel=3,
            )

is_cog property

True iff the backing file on disk is a valid COG.

False for MEM datasets, /vsimem/ paths, and unsaved datasets (empty :attr:file_name).

Examples:

• Check the backing file of a newly-opened COG:

  >>> from pyramids.dataset import Dataset  # doctest: +SKIP
  >>> ds = Dataset.read_file("scene.tif")  # doctest: +SKIP
  >>> ds.is_cog  # doctest: +SKIP
  True
  

• Plain GeoTIFFs and MEM datasets return False:

  >>> plain = Dataset.read_file("plain.tif")  # doctest: +SKIP
  >>> plain.is_cog  # doctest: +SKIP
  False
  

• Use in a conditional pipeline:

  >>> if not ds.is_cog:  # doctest: +SKIP
  ...     ds.to_cog("fixed.tif")
  

_is_cog_cheap(path) staticmethod

Fast, metadata-only heuristic for "is this file a COG?" (ARC-7).

Avoids the full COG validator on every is_cog access (which reads the whole IFD/offset table — costly over /vsicurl). Checks: GTiff driver, no external .ovr sidecar, internally tiled (square blocks or a single tile), and internal overviews present when the image is larger than one tile. This can FALSE-POSITIVE on a tiled GeoTIFF that is not laid out in strict COG order — use :meth:validate_cog for the authoritative check.

Parameters:

Name	Type	Description	Default
path	str	On-disk or remote /vsi* path.	required

Returns:

Name	Type	Description
bool	bool	True when the file looks like a COG by the cheap heuristic.

Source code in src/pyramids/dataset/engines/cog.py

@staticmethod
def _is_cog_cheap(path: str) -> bool:
    """Fast, metadata-only heuristic for "is this file a COG?" (ARC-7).

    Avoids the full COG validator on every `is_cog` access (which reads the
    whole IFD/offset table — costly over `/vsicurl`). Checks: GTiff driver,
    no external `.ovr` sidecar, internally tiled (square blocks or a single
    tile), and internal overviews present when the image is larger than one
    tile. This can FALSE-POSITIVE on a tiled GeoTIFF that is not laid out in
    strict COG order — use :meth:`validate_cog` for the authoritative check.

    Args:
        path: On-disk or remote `/vsi*` path.

    Returns:
        bool: `True` when the file looks like a COG by the cheap heuristic.
    """
    cfg = _resolve_read_config(path, None)
    with config_context(cfg):
        try:
            ds = gdal.Open(path)
        except RuntimeError:
            return False
        if ds is None:
            return False
        try:
            if ds.GetDriver().ShortName != "GTiff":
                return False
            files = ds.GetFileList() or []
            if any(str(f).lower().endswith(".ovr") for f in files):
                return False
            band = ds.GetRasterBand(1)
            block_x, block_y = band.GetBlockSize()
            width, height = ds.RasterXSize, ds.RasterYSize
            single_tile = block_x >= width and block_y >= height
            tiled = block_x == block_y or single_tile
            if not tiled:
                return False
            needs_overviews = max(width, height) > max(block_x, block_y)
            if needs_overviews and band.GetOverviewCount() == 0:
                return False
            return True
        finally:
            ds = None

validate_cog(strict=False, config=None)

Validate the backing file as a COG.

Parameters:

Name	Type	Description	Default
strict	bool	If True, warnings are treated as errors.	False
config	dict[str, str] | None	GDAL config options for the read; defaults to the remote read tuning for /vsicurl paths (see :func:pyramids.dataset.cog.validate.validate).	None

Returns:

Type	Description
ValidationReport	ValidationReport with errors, warnings, and structural details.

Raises:

Type	Description
FileNotFoundError	Dataset has no on-disk backing file (MEM-only or /vsimem/).

Examples:

• Validate and branch on the result:

  >>> from pyramids.dataset import Dataset  # doctest: +SKIP
  >>> ds = Dataset.read_file("scene.tif")  # doctest: +SKIP
  >>> report = ds.validate_cog()  # doctest: +SKIP
  >>> bool(report)  # doctest: +SKIP
  True
  

• Strict mode promotes warnings to errors:

  >>> strict = ds.validate_cog(strict=True)  # doctest: +SKIP
  >>> if not strict:  # doctest: +SKIP
  ...     for err in strict.errors: print(err)
  

• Inspect structural details from the report:

  >>> report.details.get("blocksize")  # doctest: +SKIP
  [512, 512]
  

Source code in src/pyramids/dataset/engines/cog.py

@under_gdal_env
def validate_cog(
    self, strict: bool = False, config: dict[str, str] | None = None
) -> ValidationReport:
    """Validate the backing file as a COG.

    Args:
        strict: If `True`, warnings are treated as errors.
        config: GDAL config options for the read; defaults to the remote
            read tuning for `/vsicurl` paths (see
            :func:`pyramids.dataset.cog.validate.validate`).

    Returns:
        ValidationReport with errors, warnings, and structural details.

    Raises:
        FileNotFoundError: Dataset has no on-disk backing file
            (MEM-only or `/vsimem/`).

    Examples:
        - Validate and branch on the result:
            ```python
            >>> from pyramids.dataset import Dataset  # doctest: +SKIP
            >>> ds = Dataset.read_file("scene.tif")  # doctest: +SKIP
            >>> report = ds.validate_cog()  # doctest: +SKIP
            >>> bool(report)  # doctest: +SKIP
            True

            ```
        - Strict mode promotes warnings to errors:
            ```python
            >>> strict = ds.validate_cog(strict=True)  # doctest: +SKIP
            >>> if not strict:  # doctest: +SKIP
            ...     for err in strict.errors: print(err)

            ```
        - Inspect structural details from the report:
            ```python
            >>> report.details.get("blocksize")  # doctest: +SKIP
            [512, 512]

            ```
    """
    fn = self._on_disk_path()
    if fn is None:
        raise FileNotFoundError(
            "Dataset has no on-disk backing file to validate "
            "(is this a MEM or /vsimem/ dataset?)"
        )
    return validate(fn, strict=strict, config=config)

info(config=None)

Return structured COG metadata for the backing file.

Reads only headers/metadata (no pixels) and reports compression, predictor, blocksize, dtype, CRS/bounds/resolution, the overview pyramid, per-band tags, and colour-table presence. See :class:pyramids.dataset.cog.inspect.COGInfo.

Parameters:

Name	Type	Description	Default
config	dict[str, str] | None	GDAL config options for the read; defaults to the remote read tuning for /vsicurl paths.	None

Returns:

Name	Type	Description
COGInfo	COGInfo	The structured metadata for the on-disk file.

Raises:

Type	Description
FileNotFoundError	Dataset has no on-disk backing file (MEM-only or /vsimem/).

Examples:

• Inspect a COG's compression and overview pyramid:

  >>> from pyramids.dataset import Dataset  # doctest: +SKIP
  >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
  >>> info = ds.cog_info()  # doctest: +SKIP
  >>> info.compression  # doctest: +SKIP
  'DEFLATE'
  >>> [o.decimation for o in info.overviews]  # doctest: +SKIP
  [2, 4, 8]
  

• Read the tile size and band count:

  >>> info.blocksize  # doctest: +SKIP
  (512, 512)
  >>> info.band_count  # doctest: +SKIP
  1
  

Source code in src/pyramids/dataset/engines/cog.py

@under_gdal_env
def info(self, config: dict[str, str] | None = None) -> COGInfo:
    """Return structured COG metadata for the backing file.

    Reads only headers/metadata (no pixels) and reports compression,
    predictor, blocksize, dtype, CRS/bounds/resolution, the overview
    pyramid, per-band tags, and colour-table presence. See
    :class:`pyramids.dataset.cog.inspect.COGInfo`.

    Args:
        config: GDAL config options for the read; defaults to the remote
            read tuning for `/vsicurl` paths.

    Returns:
        COGInfo: The structured metadata for the on-disk file.

    Raises:
        FileNotFoundError: Dataset has no on-disk backing file
            (MEM-only or `/vsimem/`).

    Examples:
        - Inspect a COG's compression and overview pyramid:
            ```python
            >>> from pyramids.dataset import Dataset  # doctest: +SKIP
            >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
            >>> info = ds.cog_info()  # doctest: +SKIP
            >>> info.compression  # doctest: +SKIP
            'DEFLATE'
            >>> [o.decimation for o in info.overviews]  # doctest: +SKIP
            [2, 4, 8]

            ```
        - Read the tile size and band count:
            ```python
            >>> info.blocksize  # doctest: +SKIP
            (512, 512)
            >>> info.band_count  # doctest: +SKIP
            1

            ```
    """
    fn = self._on_disk_path()
    if fn is None:
        raise FileNotFoundError(
            "Dataset has no on-disk backing file to inspect "
            "(is this a MEM or /vsimem/ dataset?)"
        )
    return cog_info(fn, config=config)

read_part(bbox, *, dst_width=None, dst_height=None, bbox_crs=None, resampling='bilinear', band=None, return_transform=False)

read_part(bbox: tuple[float, float, float, float], *, dst_width: int | None = ..., dst_height: int | None = ..., bbox_crs: int | str | None = ..., resampling: str = ..., band: int | None = ..., return_transform: Literal[False] = ...) -> np.typing.NDArray

read_part(bbox: tuple[float, float, float, float], *, dst_width: int | None = ..., dst_height: int | None = ..., bbox_crs: int | str | None = ..., resampling: str = ..., band: int | None = ..., return_transform: Literal[True]) -> tuple[np.typing.NDArray, tuple[float, float, float, float, float, float]]

read_part(bbox: tuple[float, float, float, float], *, dst_width: int | None = ..., dst_height: int | None = ..., bbox_crs: int | str | None = ..., resampling: str = ..., band: int | None = ..., return_transform: bool) -> np.typing.NDArray | tuple[np.typing.NDArray, tuple[float, float, float, float, float, float]]

Read a geographic window, decimated from the nearest overview.

Requesting a dst_width/dst_height smaller than the source window makes GDAL serve the data from the nearest overview level, so for a COG over /vsicurl/ only the relevant byte ranges are fetched — the cloud-native partial-read pattern.

The window answers in physical units, like read_array: a CF-packed band (scale_factor / add_offset) is unpacked and comes back float64, whether the window lies inside the raster or only partly overlaps it (see the note below).

Parameters:

Name	Type	Description	Default
bbox	tuple[float, float, float, float]	(min_x, min_y, max_x, max_y) window in bbox_crs.	required
dst_width	int | None	Output width in pixels. Defaults to the source window width (no decimation).	None
dst_height	int | None	Output height in pixels. Defaults to the source window height.	None
bbox_crs	int | str | None	CRS of bbox, reprojected to the dataset CRS when it differs. Defaults to None, meaning the bbox is already in the raster's own coordinates, so nothing is transformed.	None
resampling	str	Resampling method, case-insensitive. One of nearest, bilinear, cubic, cubicspline (alias cubic_spline), lanczos, average, mode, plus gauss and rms when the GDAL build provides them.	'bilinear'
band	int | None	0-based band index. None reads all bands.	None
return_transform	bool	When True, return (array, geotransform) instead of the bare array, so a caller can place the result without reproducing the pixel-snapping arithmetic. The geotransform is the six-tuple GetGeoTransform uses, describing the returned buffer at its output resolution: it maps output pixel coordinates to the dataset's CRS, over the window actually read (the requested bbox snapped outward to whole source pixels). Its origin, cell size and orientation place the buffer's outer extent exactly, and that extent always matches the snapped window. Per-cell placement -- each cell's data sitting on its own transform centre -- is exact only for a read at native resolution (dst_width/dst_height omitted, so the output keeps the snapped window's size). A resampled read can shift a cell's data off its transform centre: by a fraction of a cell in the interior, and by more near the raster edge, where a window that overruns the edge has its out-of-raster remainder padded before resampling -- reaching several output cells when such a window is upsampled. The outer extent stays exact in every case; read at native resolution when you need cell-exact data. Defaults to False, which keeps the bare-array return unchanged.	False

Returns:

Type	Description
NDArray | tuple[NDArray, tuple[float, float, float, float, float, float]]	numpy.ndarray: (rows, cols) for a single band, or
NDArray | tuple[NDArray, tuple[float, float, float, float, float, float]]	(bands, rows, cols) for all bands; always sized
NDArray | tuple[NDArray, tuple[float, float, float, float, float, float]]	dst_height x dst_width (the requested output size). Pixel values
NDArray | tuple[NDArray, tuple[float, float, float, float, float, float]]	only — no transform, bounds, or CRS is attached, unless
NDArray | tuple[NDArray, tuple[float, float, float, float, float, float]]	return_transform=True, in which case a
NDArray | tuple[NDArray, tuple[float, float, float, float, float, float]]	(array, geotransform) tuple is returned instead (see
NDArray | tuple[NDArray, tuple[float, float, float, float, float, float]]	return_transform).

Raises:

Type	Description
CRSError	An explicit bbox_crs was given but the raster has no CRS to transform into. Omit it to read in the raster's own coordinates (ARC-26).
TypeError	resampling is not a string.
ValueError	Unknown resampling, or dst_width / dst_height is given as a non-positive pixel count.
OutOfBoundsError	The window does not intersect the raster at all.

Note

A window that only partially overlaps the raster is not stretched to fill the output: the intersection is read and placed at its correct offset inside a dst_height x dst_width buffer whose out-of-raster remainder is filled with NoData (the band's NoData value, else NaN for float / 0 for integer — see :meth:_nodata_fill). A fully-inside window is returned without padding. This keeps the result aligned to the requested window, which matters for edge tiles served by :meth:read_tile.

The padded partial-overlap read is assembled in stored units -- the stored values, padded with the stored NoData -- and then unpacked once, like a fully-inside read. So on a CF-packed band it is in the same physical units and dtype as any other window, and on a band that declares NoData the padding reads as its gap does in a default read (the sentinel unpacked with the data, e.g. -9999 at scale=0.5 as -4999.5).

Examples:

• Read a 256x256 decimated thumbnail of a bbox:

  >>> from pyramids.dataset import Dataset  # doctest: +SKIP
  >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
  >>> arr = ds.read_part(  # doctest: +SKIP
  ...     (12.4, 41.8, 12.6, 42.0), dst_width=256, dst_height=256,
  ... )
  >>> arr.shape[-2:]  # doctest: +SKIP
  (256, 256)
  

• A window inside a CF-packed raster comes back in physical units:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset, GeoReference
  >>> packed = Dataset.from_array(
  ...     np.full((4, 4), 100, dtype="int16"),
  ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
  ... )
  >>> packed.scale = [0.01]
  >>> part = packed.read_part((0, -2, 2, 0), band=0)
  >>> part.shape, part.dtype, float(part.max())
  ((2, 2), dtype('float64'), 1.0)
  

• A window straddling the raster's edge is padded with the NoData and unpacked with the data, so it is in the same units as one inside:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset, GeoReference
  >>> packed = Dataset.from_array(
  ...     np.full((4, 4), 100, dtype="int16"),
  ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
  ...     no_data_value=-9999,
  ... )
  >>> packed.scale = [0.5]
  >>> edge = packed.read_part((-1, -2, 1, 0), band=0)
  >>> edge.dtype, edge.tolist()
  (dtype('float64'), [[-4999.5, 50.0], [-4999.5, 50.0]])
  

• Ask for the geotransform of the window actually read. The bbox is snapped outward to whole source pixels, so the transform, not the requested bbox, is what places the cells:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset, GeoReference
  >>> grid = Dataset.from_array(
  ...     np.tile(np.arange(8, dtype="float64"), (8, 1)),
  ...     geo_ref=GeoReference(top_left_corner=(0, 8), cell_size=1.0, epsg=3857),
  ... )
  >>> arr, gt = grid.read_part(
  ...     (1.5, 1.5, 4.5, 4.5), dst_width=3, dst_height=3, return_transform=True,
  ... )
  >>> arr.shape
  (3, 3)
  >>> gt  # origin snapped to x=1, cells 4/3 wide
  (1.0, 1.3333333333333333, 0.0, 5.0, 0.0, -1.3333333333333333)
  

Source code in src/pyramids/dataset/engines/cog.py

@under_gdal_env
def read_part(
    self,
    bbox: tuple[float, float, float, float],
    *,
    dst_width: int | None = None,
    dst_height: int | None = None,
    bbox_crs: int | str | None = None,
    resampling: str = "bilinear",
    band: int | None = None,
    return_transform: bool = False,
) -> (
    np.typing.NDArray
    | tuple[np.typing.NDArray, tuple[float, float, float, float, float, float]]
):
    """Read a geographic window, decimated from the nearest overview.

    Requesting a `dst_width`/`dst_height` smaller than the source window
    makes GDAL serve the data from the nearest overview level, so for a COG
    over `/vsicurl/` only the relevant byte ranges are fetched — the
    cloud-native partial-read pattern.

    The window answers in physical units, like `read_array`: a CF-packed band
    (`scale_factor` / `add_offset`) is unpacked and comes back `float64`, whether
    the window lies inside the raster or only partly overlaps it (see the note
    below).

    Args:
        bbox: `(min_x, min_y, max_x, max_y)` window in `bbox_crs`.
        dst_width: Output width in pixels. Defaults to the source window
            width (no decimation).
        dst_height: Output height in pixels. Defaults to the source window
            height.
        bbox_crs: CRS of `bbox`, reprojected to the dataset CRS when it
            differs. Defaults to `None`, meaning the bbox is already in the
            raster's own coordinates, so nothing is transformed.
        resampling: Resampling method, case-insensitive. One of `nearest`,
            `bilinear`, `cubic`, `cubicspline` (alias `cubic_spline`),
            `lanczos`, `average`, `mode`, plus `gauss` and `rms` when the
            GDAL build provides them.
        band: 0-based band index. `None` reads all bands.
        return_transform: When `True`, return `(array, geotransform)`
            instead of the bare array, so a caller can place the result
            without reproducing the pixel-snapping arithmetic. The
            geotransform is the six-tuple `GetGeoTransform` uses, describing
            the returned buffer at its output resolution: it maps output
            pixel coordinates to the **dataset's** CRS, over the window
            actually read (the requested bbox snapped outward to whole
            source pixels). Its origin, cell size and orientation place the
            buffer's outer extent exactly, and that extent always matches the
            snapped window. Per-cell placement -- each cell's data sitting on
            its own transform centre -- is exact only for a read at native
            resolution (`dst_width`/`dst_height` omitted, so the output keeps
            the snapped window's size). A resampled read can shift a cell's
            data off its transform centre: by a fraction of a cell in the
            interior, and by more near the raster edge, where a window that
            overruns the edge has its out-of-raster remainder padded before
            resampling -- reaching several output cells when such a window is
            upsampled. The outer extent stays exact in every case; read at
            native resolution when you need cell-exact data. Defaults to
            `False`, which keeps the bare-array return unchanged.

    Returns:
        numpy.ndarray: `(rows, cols)` for a single band, or
        `(bands, rows, cols)` for all bands; always sized
        `dst_height x dst_width` (the requested output size). Pixel values
        only — no transform, bounds, or CRS is attached, unless
        `return_transform=True`, in which case a
        `(array, geotransform)` tuple is returned instead (see
        `return_transform`).

    Raises:
        CRSError: An explicit `bbox_crs` was given but the raster has no CRS
            to transform into. Omit it to read in the raster's own
            coordinates (ARC-26).
        TypeError: `resampling` is not a string.
        ValueError: Unknown `resampling`, or `dst_width` / `dst_height` is
            given as a non-positive pixel count.
        OutOfBoundsError: The window does not intersect the raster at all.

    Note:
        A window that only **partially** overlaps the raster is **not**
        stretched to fill the output: the intersection is read and placed
        at its correct offset inside a `dst_height x dst_width` buffer
        whose out-of-raster remainder is filled with NoData (the band's
        NoData value, else NaN for float / `0` for integer — see
        :meth:`_nodata_fill`). A fully-inside window is returned without
        padding. This keeps the result aligned to the requested window,
        which matters for edge tiles served by :meth:`read_tile`.

        The padded partial-overlap read is assembled in stored units --
        the stored values, padded with the stored NoData -- and then
        unpacked once, like a fully-inside read. So on a CF-packed band it
        is in the same physical units and dtype as any other window, and
        on a band that declares NoData the padding reads as its gap does in
        a default read (the sentinel unpacked with the data, e.g. `-9999`
        at `scale=0.5` as `-4999.5`).

    Examples:
        - Read a 256x256 decimated thumbnail of a bbox:
            ```python
            >>> from pyramids.dataset import Dataset  # doctest: +SKIP
            >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
            >>> arr = ds.read_part(  # doctest: +SKIP
            ...     (12.4, 41.8, 12.6, 42.0), dst_width=256, dst_height=256,
            ... )
            >>> arr.shape[-2:]  # doctest: +SKIP
            (256, 256)

            ```
        - A window inside a CF-packed raster comes back in physical units:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset, GeoReference
            >>> packed = Dataset.from_array(
            ...     np.full((4, 4), 100, dtype="int16"),
            ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
            ... )
            >>> packed.scale = [0.01]
            >>> part = packed.read_part((0, -2, 2, 0), band=0)
            >>> part.shape, part.dtype, float(part.max())
            ((2, 2), dtype('float64'), 1.0)

            ```
        - A window straddling the raster's edge is padded with the NoData and
          unpacked with the data, so it is in the same units as one inside:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset, GeoReference
            >>> packed = Dataset.from_array(
            ...     np.full((4, 4), 100, dtype="int16"),
            ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
            ...     no_data_value=-9999,
            ... )
            >>> packed.scale = [0.5]
            >>> edge = packed.read_part((-1, -2, 1, 0), band=0)
            >>> edge.dtype, edge.tolist()
            (dtype('float64'), [[-4999.5, 50.0], [-4999.5, 50.0]])

            ```
        - Ask for the geotransform of the window actually read. The bbox is
          snapped outward to whole source pixels, so the transform, not the
          requested bbox, is what places the cells:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset, GeoReference
            >>> grid = Dataset.from_array(
            ...     np.tile(np.arange(8, dtype="float64"), (8, 1)),
            ...     geo_ref=GeoReference(top_left_corner=(0, 8), cell_size=1.0, epsg=3857),
            ... )
            >>> arr, gt = grid.read_part(
            ...     (1.5, 1.5, 4.5, 4.5), dst_width=3, dst_height=3, return_transform=True,
            ... )
            >>> arr.shape
            (3, 3)
            >>> gt  # origin snapped to x=1, cells 4/3 wide
            (1.0, 1.3333333333333333, 0.0, 5.0, 0.0, -1.3333333333333333)

            ```
    """
    alg = _resolve_read_resampling(resampling)
    self._reject_nonpositive_output_sizes(dst_width, dst_height)
    # This serves a decimated window from the source; a NetCDF multidim view can't be window-read
    # by GDAL >= 3.13, so materialise it first (no-op for an ordinary raster).
    self._ds._materialize_md_view()
    ds = self._ds._raster
    min_x, min_y, max_x, max_y = self._reproject_bbox(bbox, bbox_crs)
    geotransform = ds.GetGeoTransform()
    px_tl, py_tl = world_to_pixel(geotransform, min_x, max_y)
    px_br, py_br = world_to_pixel(geotransform, max_x, min_y)

    # The full requested window, in source pixel coordinates (may extend
    # beyond the raster on any side).
    req_xoff = int(math.floor(min(px_tl, px_br)))
    req_yoff = int(math.floor(min(py_tl, py_br)))
    req_xsize = int(math.ceil(max(px_tl, px_br))) - req_xoff
    req_ysize = int(math.ceil(max(py_tl, py_br))) - req_yoff
    if req_xsize <= 0 or req_ysize <= 0:
        raise OutOfBoundsError(
            f"bbox {bbox} (crs {bbox_crs}) has zero pixel extent"
        )

    # Intersection of the requested window with the raster.
    ix0 = max(0, req_xoff)
    iy0 = max(0, req_yoff)
    ix1 = min(ds.RasterXSize, req_xoff + req_xsize)
    iy1 = min(ds.RasterYSize, req_yoff + req_ysize)
    if ix1 - ix0 <= 0 or iy1 - iy0 <= 0:
        raise OutOfBoundsError(
            f"bbox {bbox} (crs {bbox_crs}) does not intersect the raster"
        )

    out_w = dst_width if dst_width is not None else req_xsize
    out_h = dst_height if dst_height is not None else req_ysize
    # Describes the buffer both arms below return: the snapped window at the
    # output resolution, in the dataset's CRS. Computed here, before the read,
    # because it depends only on the window and the output size, not the pixels.
    out_gt = self._output_geotransform(
        geotransform, req_xoff, req_yoff, req_xsize, req_ysize, out_w, out_h
    )
    source = ds if band is None else ds.GetRasterBand(band + 1)

    fully_inside = (ix0, iy0, ix1, iy1) == (
        req_xoff,
        req_yoff,
        req_xoff + req_xsize,
        req_yoff + req_ysize,
    )
    if fully_inside:
        # Physical units, like `read_array`. These are overview-decimated reads
        # of the same band, and a raster that renders at one scale through
        # `plot()` and another through `plot(overview=True)` is a defect a
        # reader can see.
        array = np.asarray(
            self._ds.io._apply_scale_offset(
                np.asarray(
                    source.ReadAsArray(
                        ix0,
                        iy0,
                        ix1 - ix0,
                        iy1 - iy0,
                        buf_xsize=out_w,
                        buf_ysize=out_h,
                        resample_alg=alg,
                    )
                ),
                band,
            )
        )
        return (array, out_gt) if return_transform else array

    # Partial overlap: read only the intersection, then place it at its
    # correct offset inside a full-size output buffer padded with NoData,
    # so the returned array stays aligned to the requested window.
    scale_x = out_w / req_xsize
    scale_y = out_h / req_ysize
    ox0 = max(0, min(out_w, int(round((ix0 - req_xoff) * scale_x))))
    oy0 = max(0, min(out_h, int(round((iy0 - req_yoff) * scale_y))))
    ox1 = max(ox0 + 1, min(out_w, int(round((ix1 - req_xoff) * scale_x))))
    oy1 = max(oy0 + 1, min(out_h, int(round((iy1 - req_yoff) * scale_y))))
    sub = np.asarray(
        source.ReadAsArray(
            ix0,
            iy0,
            ix1 - ix0,
            iy1 - iy0,
            buf_xsize=ox1 - ox0,
            buf_ysize=oy1 - oy0,
            resample_alg=alg,
        )
    )
    fill = self._nodata_fill(ds.GetRasterBand(1))
    if sub.ndim == 3:
        out = np.full((sub.shape[0], out_h, out_w), fill, dtype=sub.dtype)
        out[:, oy0:oy1, ox0:ox1] = sub
    else:
        out = np.full((out_h, out_w), fill, dtype=sub.dtype)
        out[oy0:oy1, ox0:ox1] = sub
    # Assembled in stored units -- the counts and the stored sentinel padding them
    # -- and unpacked once, so the padding reads as a gap exactly the way the
    # fully-inside arm above reads its own. Returning the stored array made a window
    # straddling the raster's edge come back as `int16` counts while a window just
    # inside it came back as physical `float64`, and every edge `read_tile` took
    # this arm.
    array = np.asarray(self._ds.io._apply_scale_offset(out, band))
    return (array, out_gt) if return_transform else array

preview(*, max_size=1024, resampling='bilinear', band=None)

Read a whole-image thumbnail downsampled to max_size on the long edge.

Pulls from a coarse overview when one exists, so previewing a huge COG is cheap. The thumbnail is the same band at a coarser sampling, so it answers in the same units as read_array: a CF-packed band (scale_factor / add_offset) is unpacked to physical values.

Parameters:

Name	Type	Description	Default
max_size	int	Maximum pixels on the longer edge. Defaults to 1024.	1024
resampling	str	Resampling method (see :meth:read_part).	'bilinear'
band	int | None	0-based band index. None reads all bands.	None

Returns:

Type	Description
NDArray	numpy.ndarray: The downsampled array, (rows, cols) or
NDArray	(bands, rows, cols), float64 when a band that was read is packed. Pixel
NDArray	values only — no transform, bounds, or CRS is attached to the
NDArray	returned array.

Raises:

Type	Description
TypeError	resampling is not a string.
ValueError	Unknown resampling.

Examples:

• Build a 128px thumbnail of a single band:

  >>> from pyramids.dataset import Dataset  # doctest: +SKIP
  >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
  >>> thumb = ds.preview(max_size=128, band=0)  # doctest: +SKIP
  >>> max(thumb.shape)  # doctest: +SKIP
  128
  

• A thumbnail of a CF-packed raster is in physical units:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset, GeoReference
  >>> packed = Dataset.from_array(
  ...     np.full((4, 4), 100, dtype="int16"),
  ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
  ... )
  >>> packed.scale = [0.01]
  >>> thumb = packed.preview(max_size=2, band=0)
  >>> thumb.shape, float(thumb[0, 0])
  ((2, 2), 1.0)
  

Source code in src/pyramids/dataset/engines/cog.py

@under_gdal_env
def preview(
    self,
    *,
    max_size: int = 1024,
    resampling: str = "bilinear",
    band: int | None = None,
) -> np.typing.NDArray:
    """Read a whole-image thumbnail downsampled to `max_size` on the long edge.

    Pulls from a coarse overview when one exists, so previewing a huge COG
    is cheap. The thumbnail is the same band at a coarser sampling, so it
    answers in the same units as `read_array`: a CF-packed band
    (`scale_factor` / `add_offset`) is unpacked to physical values.

    Args:
        max_size: Maximum pixels on the longer edge. Defaults to 1024.
        resampling: Resampling method (see :meth:`read_part`).
        band: 0-based band index. `None` reads all bands.

    Returns:
        numpy.ndarray: The downsampled array, `(rows, cols)` or
        `(bands, rows, cols)`, `float64` when a band that was read is packed. Pixel
        values only — no transform, bounds, or CRS is attached to the
        returned array.

    Raises:
        TypeError: `resampling` is not a string.
        ValueError: Unknown `resampling`.

    Examples:
        - Build a 128px thumbnail of a single band:
            ```python
            >>> from pyramids.dataset import Dataset  # doctest: +SKIP
            >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
            >>> thumb = ds.preview(max_size=128, band=0)  # doctest: +SKIP
            >>> max(thumb.shape)  # doctest: +SKIP
            128

            ```
        - A thumbnail of a CF-packed raster is in physical units:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset, GeoReference
            >>> packed = Dataset.from_array(
            ...     np.full((4, 4), 100, dtype="int16"),
            ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
            ... )
            >>> packed.scale = [0.01]
            >>> thumb = packed.preview(max_size=2, band=0)
            >>> thumb.shape, float(thumb[0, 0])
            ((2, 2), 1.0)

            ```
    """
    alg = _resolve_read_resampling(resampling)
    width, height = self._ds.columns, self._ds.rows
    scale = max(width, height) / max_size
    if scale <= 1:
        out_w, out_h = width, height
    else:
        out_w, out_h = max(1, round(width / scale)), max(1, round(height / scale))
    ds = self._ds._raster
    source = ds if band is None else ds.GetRasterBand(band + 1)
    # Physical units, like `read_array`: `preview` is the same band at a coarser
    # sampling, not a different quantity.
    return np.asarray(
        self._ds.io._apply_scale_offset(
            np.asarray(
                source.ReadAsArray(
                    buf_xsize=out_w, buf_ysize=out_h, resample_alg=alg
                )
            ),
            band,
        )
    )

point(x, y, *, point_crs=None, band=None)

Sample band value(s) at a single coordinate.

The sample answers in the same units as read_array: a CF-packed band (scale_factor / add_offset) is unpacked, so point and read_array()[row, col] agree. A no-data cell answers its stored sentinel transformed like any other value; compare against no_data_value on an unpack=False read instead.

Parameters:

Name	Type	Description	Default
x	float	X / longitude / easting in point_crs.	required
y	float	Y / latitude / northing in point_crs.	required
point_crs	int | str | None	CRS of (x, y), reprojected to the dataset CRS when it differs. Defaults to None, meaning the coordinates are already in the raster's own CRS, so nothing is transformed.	None
band	int | None	0-based band index. None samples all bands.	None

Returns:

Type	Description
NDArray	numpy.ndarray: A scalar 0-d array for a single band, or a
NDArray	(bands,) array when band is None, float64 when a band that was read
NDArray	is packed. Pixel values only — no coordinate metadata is attached.

Raises:

Type	Description
CRSError	An explicit point_crs was given but the raster has no CRS to transform into. Omit it to read in the raster's own coordinates (ARC-26).
OutOfBoundsError	The point falls outside the raster extent.

Examples:

• Sample all bands at a lon/lat coordinate:

  >>> from pyramids.dataset import Dataset  # doctest: +SKIP
  >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
  >>> ds.point(12.5, 41.9)  # doctest: +SKIP
  array([1234.], dtype=float32)
  

• On a CF-packed raster the sample matches the default read_array:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset, GeoReference
  >>> packed = Dataset.from_array(
  ...     np.array([[100, 200], [300, 400]], dtype="int16"),
  ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
  ... )
  >>> packed.scale = [0.01]
  >>> float(packed.point(1.5, -0.5, band=0)), float(packed.read_array()[0, 1])
  (2.0, 2.0)
  

Source code in src/pyramids/dataset/engines/cog.py

@under_gdal_env
def point(
    self,
    x: float,
    y: float,
    *,
    point_crs: int | str | None = None,
    band: int | None = None,
) -> np.typing.NDArray:
    """Sample band value(s) at a single coordinate.

    The sample answers in the same units as `read_array`: a CF-packed band
    (`scale_factor` / `add_offset`) is unpacked, so `point` and
    `read_array()[row, col]` agree. A no-data cell answers its stored
    sentinel transformed like any other value; compare against
    `no_data_value` on an `unpack=False` read instead.

    Args:
        x: X / longitude / easting in `point_crs`.
        y: Y / latitude / northing in `point_crs`.
        point_crs: CRS of `(x, y)`, reprojected to the dataset CRS when it
            differs. Defaults to `None`, meaning the coordinates are already
            in the raster's own CRS, so nothing is transformed.
        band: 0-based band index. `None` samples all bands.

    Returns:
        numpy.ndarray: A scalar 0-d array for a single band, or a
        `(bands,)` array when `band` is `None`, `float64` when a band that was read
        is packed. Pixel values only — no coordinate metadata is attached.

    Raises:
        CRSError: An explicit `point_crs` was given but the raster has no
            CRS to transform into. Omit it to read in the raster's own
            coordinates (ARC-26).
        OutOfBoundsError: The point falls outside the raster extent.

    Examples:
        - Sample all bands at a lon/lat coordinate:
            ```python
            >>> from pyramids.dataset import Dataset  # doctest: +SKIP
            >>> ds = Dataset.read_file("scene_cog.tif")  # doctest: +SKIP
            >>> ds.point(12.5, 41.9)  # doctest: +SKIP
            array([1234.], dtype=float32)

            ```
        - On a CF-packed raster the sample matches the default `read_array`:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset, GeoReference
            >>> packed = Dataset.from_array(
            ...     np.array([[100, 200], [300, 400]], dtype="int16"),
            ...     geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=1.0, epsg=4326),
            ... )
            >>> packed.scale = [0.01]
            >>> float(packed.point(1.5, -0.5, band=0)), float(packed.read_array()[0, 1])
            (2.0, 2.0)

            ```
    """
    col, row = self._world_to_pixel(x, y, point_crs)
    if not (0 <= col < self._ds.columns and 0 <= row < self._ds.rows):
        raise OutOfBoundsError(
            f"point ({x}, {y}) in crs {point_crs} is outside the raster extent"
        )
    ds = self._ds._raster
    source = ds if band is None else ds.GetRasterBand(band + 1)
    # Unpacked while the read still has its `(bands, 1, 1)` shape, so the
    # per-band broadcast works; the reshape to a plain value comes after. A
    # `point` that answered 100 where `read_array()[row, col]` answered 2.5 was
    # the same raster sampled two ways.
    arr = np.asarray(
        self._ds.io._apply_scale_offset(
            np.asarray(source.ReadAsArray(col, row, 1, 1)), band
        )
    )
    return arr.reshape(-1) if band is None else arr.reshape(())

read_tile(z, x, y, *, tilesize=256, resampling='bilinear', band=None)

Read a Web-Mercator XYZ/slippy-map tile.

Computes the EPSG:3857 bounds of tile (z, x, y) from the closed-form Web-Mercator formula and delegates to :meth:read_part at tilesize resolution — no extra tiling dependency needed.

Parameters:

Name	Type	Description	Default
z	int	Zoom level.	required
x	int	Tile column index.	required
y	int	Tile row index (origin top-left / north-west).	required
tilesize	int	Output tile size in pixels (square). Defaults to 256.	256
resampling	str	Resampling method (see :meth:read_part).	'bilinear'
band	int | None	0-based band index. None reads all bands.	None

Returns:

Type	Description
NDArray	numpy.ndarray: A (tilesize, tilesize) or
NDArray	(bands, tilesize, tilesize) array, in physical units like
NDArray	read_array (float64 when a band read is CF-packed). Pixel values
NDArray	only — the tile's georeferencing is defined by its (z, x, y), not
NDArray	attached to the array; edge tiles are NoData-padded and unpacked in
NDArray	the same units as inner ones (see :meth:read_part).

Raises:

Type	Description
OutOfBoundsError	The tile does not intersect the raster.

Examples:

• Read the zoom-0 world tile of a global COG:

  >>> from pyramids.dataset import Dataset  # doctest: +SKIP
  >>> ds = Dataset.read_file("global_cog.tif")  # doctest: +SKIP
  >>> tile = ds.read_tile(0, 0, 0)  # doctest: +SKIP
  >>> tile.shape[-2:]  # doctest: +SKIP
  (256, 256)
  

Source code in src/pyramids/dataset/engines/cog.py

@under_gdal_env
def read_tile(
    self,
    z: int,
    x: int,
    y: int,
    *,
    tilesize: int = 256,
    resampling: str = "bilinear",
    band: int | None = None,
) -> np.typing.NDArray:
    """Read a Web-Mercator XYZ/slippy-map tile.

    Computes the EPSG:3857 bounds of tile `(z, x, y)` from the closed-form
    Web-Mercator formula and delegates to :meth:`read_part` at `tilesize`
    resolution — no extra tiling dependency needed.

    Args:
        z: Zoom level.
        x: Tile column index.
        y: Tile row index (origin top-left / north-west).
        tilesize: Output tile size in pixels (square). Defaults to 256.
        resampling: Resampling method (see :meth:`read_part`).
        band: 0-based band index. `None` reads all bands.

    Returns:
        numpy.ndarray: A `(tilesize, tilesize)` or
        `(bands, tilesize, tilesize)` array, in physical units like
        `read_array` (`float64` when a band read is CF-packed). Pixel values
        only — the tile's georeferencing is defined by its `(z, x, y)`, not
        attached to the array; edge tiles are NoData-padded and unpacked in
        the same units as inner ones (see :meth:`read_part`).

    Raises:
        OutOfBoundsError: The tile does not intersect the raster.

    Examples:
        - Read the zoom-0 world tile of a global COG:
            ```python
            >>> from pyramids.dataset import Dataset  # doctest: +SKIP
            >>> ds = Dataset.read_file("global_cog.tif")  # doctest: +SKIP
            >>> tile = ds.read_tile(0, 0, 0)  # doctest: +SKIP
            >>> tile.shape[-2:]  # doctest: +SKIP
            (256, 256)

            ```
    """
    bounds = _xyz_bounds_3857(z, x, y)
    return self.read_part(
        bounds,
        dst_width=tilesize,
        dst_height=tilesize,
        bbox_crs=3857,
        resampling=resampling,
        band=band,
    )

Validation

pyramids.dataset.cog.validate

COG validation wrapping osgeo_utils sample validator.

Provides :func:validate — a thin wrapper over osgeo_utils.samples.validate_cloud_optimized_geotiff.validate (GDAL ships it as a "sample"; the signature has drifted between GDAL 3.4 / 3.6 / 3.8 / 3.12, so we defensively probe the return shape). If the import fails entirely, a minimal in-house fallback checks that the file is tiled and has overviews.

Returns a :class:ValidationReport — a frozen dataclass usable as a :class:bool (is_valid) with errors, warnings, and details fields for richer reporting.

ValidationReport dataclass

Outcome of validating whether a file is a Cloud Optimized GeoTIFF.

Attributes:

Name	Type	Description
is_valid	bool	True iff :attr:errors is empty (and, under strict=True, no warnings either).
errors	list[str]	Error messages (empty when valid).
warnings	list[str]	Non-fatal warnings (e.g., "no overviews").
details	dict[str, Any]	Structural metadata from the validator — typically ifd_offsets, data_offsets, and, in the fallback path, blocksize and overview_count.

Source code in src/pyramids/dataset/cog/validate.py

@dataclass(frozen=True)
class ValidationReport:
    """Outcome of validating whether a file is a Cloud Optimized GeoTIFF.

    Attributes:
        is_valid: `True` iff :attr:`errors` is empty (and, under
            `strict=True`, no warnings either).
        errors: Error messages (empty when valid).
        warnings: Non-fatal warnings (e.g., "no overviews").
        details: Structural metadata from the validator — typically
            `ifd_offsets`, `data_offsets`, and, in the fallback
            path, `blocksize` and `overview_count`.
    """

    is_valid: bool
    errors: list[str] = field(default_factory=list)
    warnings: list[str] = field(default_factory=list)
    details: dict[str, Any] = field(default_factory=dict)

    def __bool__(self) -> bool:
        """Truthy iff the file validates as a COG.

        Examples:
            - A valid report is truthy:
                ```python
                >>> bool(ValidationReport(is_valid=True))
                True

                ```
            - An invalid report (with errors) is falsy:
                ```python
                >>> bool(ValidationReport(is_valid=False, errors=["bad"]))
                False

                ```
            - The report is usable directly in conditionals:
                ```python
                >>> report = ValidationReport(is_valid=True, details={"blocksize": [512, 512]})
                >>> "OK" if report else "bad"
                'OK'
                >>> report.details["blocksize"]
                [512, 512]

                ```
        """
        return self.is_valid

_osgeo_validate(path)

Invoke the osgeo_utils sample validator; return (errors, warnings, details).

The sample validator's signature has drifted across GDAL versions. We probe defensively: GDAL 3.6+ returns (warnings, errors, details) while older builds may return just (warnings, errors).

Parameters:

Name	Type	Description	Default
path	str	File path or /vsi* path.	required

Returns:

Type	Description
list[str]	Tuple of (errors, warnings, details) — errors listed first
list[str]	to match this module's public convention.

Raises:

Type	Description
ImportError	The osgeo_utils sample module is unavailable.
FileNotFoundError	The underlying file cannot be opened (raised via ValidateCloudOptimizedGeoTIFFException).

Source code in src/pyramids/dataset/cog/validate.py

def _osgeo_validate(
    path: str,
) -> tuple[list[str], list[str], dict[str, Any]]:
    """Invoke the osgeo_utils sample validator; return `(errors, warnings, details)`.

    The sample validator's signature has drifted across GDAL versions.
    We probe defensively: GDAL 3.6+ returns
    `(warnings, errors, details)` while older builds may return just
    `(warnings, errors)`.

    Args:
        path: File path or `/vsi*` path.

    Returns:
        Tuple of `(errors, warnings, details)` — errors listed first
        to match this module's public convention.

    Raises:
        ImportError: The `osgeo_utils` sample module is unavailable.
        FileNotFoundError: The underlying file cannot be opened
            (raised via `ValidateCloudOptimizedGeoTIFFException`).
    """
    from osgeo_utils.samples import validate_cloud_optimized_geotiff as v

    # Structural pre-check before invoking the validator — avoids
    # depending on GDAL's error-message phrasing (which varies by
    # version and locale) to detect "file not found".
    _raise_if_missing(path)

    try:
        result = v.validate(path, full_check=True)
    except v.ValidateCloudOptimizedGeoTIFFException as exc:
        # If a ValidateCloudOptimizedGeoTIFFException escapes despite
        # the pre-check, it's not about a missing file — surface it
        # as a validation error rather than letting it propagate.
        return [str(exc)], [], {}
    except RuntimeError as exc:
        # Same rationale as above for RuntimeErrors from gdal.Open
        # inside the sample validator (locale-independent fallback).
        return [str(exc)], [], {}

    errors: list[str]
    warnings: list[str]
    details: dict[str, Any]
    if len(result) == 3:
        warnings, errors, details = result
    else:  # pragma: no cover — defensive; older GDAL
        warnings, errors = result
        details = {}
    return list(errors), list(warnings), dict(details)

_fallback_validate(path)

Minimal in-house validator used when the sample module is unavailable.

Checks: file opens; image is tiled (block dimensions smaller than full extent); at least one overview present. Does NOT check the IFD-before-data layout; recommends upgrading GDAL if used.

Heuristic limitations

The "is stripped" check compares the block shape reported by :func:GetBlockSize — stripped TIFFs typically return (width, small_N) (e.g. (512, 4)) while tiled files return (tile, tile). The rule used is by!= bx and by * 4 < bx, which:

• Correctly flags standard stripped layouts ((W, 1), (W, 4), (W, 8)).
• Correctly passes square-tiled COGs ((256, 256), (512, 512)).
• Can FALSE-NEGATIVE on pathological cases such as near-square strips (by == bx) — extremely rare in practice.
• Can FALSE-POSITIVE on legitimately non-square TIFF tiles (e.g. (512, 128) used for tall elongated rasters) — also rare; the GTiff driver requires square tiles for COG.

The authoritative check is the TIFF TILEWIDTH / STRIPBYTECOUNTS tag, but reading it requires either :mod:tifffile or a direct libtiff binding. We accept the heuristic because this fallback runs only when :mod:osgeo_utils.samples.validate_cloud_optimized_geotiff is unavailable — which, in practice, is never on GDAL >= 3.4.

Parameters:

Name	Type	Description	Default
path	str	File path or /vsi* path.	required

Returns:

Type	Description
list[str]	(errors, warnings, details) — same convention as
list[str]	func:_osgeo_validate.

Source code in src/pyramids/dataset/cog/validate.py

def _fallback_validate(
    path: str,
) -> tuple[list[str], list[str], dict[str, Any]]:
    """Minimal in-house validator used when the sample module is unavailable.

    Checks: file opens; image is tiled (block dimensions smaller than
    full extent); at least one overview present. Does NOT check the
    IFD-before-data layout; recommends upgrading GDAL if used.

    Heuristic limitations:
        The "is stripped" check compares the block shape reported by
        :func:`GetBlockSize` — stripped TIFFs typically return
        `(width, small_N)` (e.g. `(512, 4)`) while tiled files
        return `(tile, tile)`. The rule used is `by!= bx and
        by * 4 < bx`, which:

        - Correctly flags standard stripped layouts (`(W, 1)`,
          `(W, 4)`, `(W, 8)`).
        - Correctly passes square-tiled COGs (`(256, 256)`,
          `(512, 512)`).
        - Can FALSE-NEGATIVE on pathological cases such as
          near-square strips (`by == bx`) — extremely rare in
          practice.
        - Can FALSE-POSITIVE on legitimately non-square TIFF tiles
          (e.g. `(512, 128)` used for tall elongated rasters) —
          also rare; the GTiff driver requires square tiles for COG.

        The authoritative check is the TIFF `TILEWIDTH` /
        `STRIPBYTECOUNTS` tag, but reading it requires either
        :mod:`tifffile` or a direct `libtiff` binding. We accept
        the heuristic because this fallback runs only when
        :mod:`osgeo_utils.samples.validate_cloud_optimized_geotiff`
        is unavailable — which, in practice, is never on GDAL >= 3.4.

    Args:
        path: File path or `/vsi*` path.

    Returns:
        `(errors, warnings, details)` — same convention as
        :func:`_osgeo_validate`.
    """
    # Locale-independent missing-file pre-check (ARC-6): each validator path
    # owns this so validate() needs no redundant outer call.
    _raise_if_missing(path)
    errors: list[str] = []
    warnings: list[str] = ["using fallback validator; osgeo_utils sample unavailable"]
    details: dict[str, Any] = {}
    ds = gdal.Open(path)
    if ds is None:
        errors.append(f"cannot open {path}")
    else:
        band = ds.GetRasterBand(1)
        bx, by = band.GetBlockSize()
        details["blocksize"] = [bx, by]
        details["overview_count"] = band.GetOverviewCount()
        # See the "Heuristic limitations" note in the docstring.
        is_stripped = by != bx and by * 4 < bx
        if is_stripped:
            errors.append("not tiled (stripped layout)")
        if band.GetOverviewCount() == 0:
            warnings.append("no overviews present")
        ds = None
    return errors, warnings, details

validate(path, strict=False, config=None)

Validate that the file at path is a valid Cloud Optimized GeoTIFF.

Delegates to osgeo_utils.samples.validate_cloud_optimized_geotiff when available (GDAL ≥ 3.4). Falls back to a minimal in-house check (tiled + overviews) when the import fails.

Parameters:

Name	Type	Description	Default
path	str | Path	Local path or /vsi* VSI path.	required
strict	bool	If True, warnings are promoted to errors.	False
config	dict[str, str] | None	GDAL config options applied (via gdal.config_options) for the duration of the validation. When None and path is a remote//vsicurl path, the :data:~pyramids.dataset.cog.options.COG_READ_DEFAULTS are applied so remote reads avoid a directory-listing round-trip.	None

Returns:

Name	Type	Description
ValidationReport	ValidationReport	Includes is_valid, error/warning lists,
	ValidationReport	and a details dict with structural metadata. Usable as a
	ValidationReport	boolean (bool(report) == report.is_valid).

Raises:

Type	Description
FileNotFoundError	When a local path does not exist. VSI paths are passed through to GDAL, which reports the error through the normal validator surface.

Examples:

• Validate a local COG and inspect the report:

  >>> from pyramids.dataset.cog import validate  # doctest: +SKIP
  >>> report = validate("scene.tif")  # doctest: +SKIP
  >>> bool(report)  # doctest: +SKIP
  True
  >>> report.details.get("blocksize")  # doctest: +SKIP
  [512, 512]
  

• Strict mode promotes warnings (e.g. "no overviews") to errors:

  >>> strict = validate("scene.tif", strict=True)  # doctest: +SKIP
  >>> if not strict:  # doctest: +SKIP
  ...     for err in strict.errors: print(err)
  

• Validate a cloud-hosted COG via VSI path:

  >>> validate("/vsis3/public-bucket/scene.tif").is_valid  # doctest: +SKIP
  True
  

Source code in src/pyramids/dataset/cog/validate.py

def validate(
    path: str | Path,
    strict: bool = False,
    config: dict[str, str] | None = None,
) -> ValidationReport:
    """Validate that the file at `path` is a valid Cloud Optimized GeoTIFF.

    Delegates to `osgeo_utils.samples.validate_cloud_optimized_geotiff`
    when available (GDAL ≥ 3.4). Falls back to a minimal in-house check
    (tiled + overviews) when the import fails.

    Args:
        path: Local path or `/vsi*` VSI path.
        strict: If `True`, warnings are promoted to errors.
        config: GDAL config options applied (via `gdal.config_options`) for
            the duration of the validation. When `None` and `path` is a
            remote/`/vsicurl` path, the
            :data:`~pyramids.dataset.cog.options.COG_READ_DEFAULTS` are
            applied so remote reads avoid a directory-listing round-trip.

    Returns:
        ValidationReport: Includes `is_valid`, error/warning lists,
        and a `details` dict with structural metadata. Usable as a
        boolean (`bool(report) == report.is_valid`).

    Raises:
        FileNotFoundError: When a *local* `path` does not exist. VSI
            paths are passed through to GDAL, which reports the error
            through the normal validator surface.

    Examples:
        - Validate a local COG and inspect the report:
            ```python
            >>> from pyramids.dataset.cog import validate  # doctest: +SKIP
            >>> report = validate("scene.tif")  # doctest: +SKIP
            >>> bool(report)  # doctest: +SKIP
            True
            >>> report.details.get("blocksize")  # doctest: +SKIP
            [512, 512]

            ```
        - Strict mode promotes warnings (e.g. "no overviews") to errors:
            ```python
            >>> strict = validate("scene.tif", strict=True)  # doctest: +SKIP
            >>> if not strict:  # doctest: +SKIP
            ...     for err in strict.errors: print(err)

            ```
        - Validate a cloud-hosted COG via VSI path:
            ```python
            >>> validate("/vsis3/public-bucket/scene.tif").is_valid  # doctest: +SKIP
            True

            ```
    """
    p = str(path)
    cfg = _resolve_read_config(p, config)
    with config_context(cfg):
        # The missing-file pre-check lives in each validator path
        # (_osgeo_validate / _fallback_validate), so no redundant outer call
        # here (ARC-6).
        try:
            errors, warnings, details = _osgeo_validate(p)
        except ImportError:  # pragma: no cover — osgeo_utils is a hard dep of GDAL
            errors, warnings, details = _fallback_validate(p)
        # Enrich with the per-overview layout (PC-4) so a report is
        # self-sufficient for debugging "why is this slow / not a COG" without
        # a second cog_info() call. Existing validator keys win on conflict.
        details = {**_probe_overview_layout(p), **dict(details)}

    if strict:
        errors = list(errors) + list(warnings)
        warnings = []

    return ValidationReport(
        is_valid=not errors,
        errors=list(errors),
        warnings=list(warnings),
        details=dict(details),
    )