Cloud Optimized GeoTIFF (COG)

The COG surface lets you write, inspect, validate, and partially read Cloud Optimized GeoTIFFs. It is GDAL-native — no extra third-party dependency — and is split between a small subpackage of helpers (pyramids.dataset.cog) and the COG engine that owns the user-facing Dataset methods.

Module layout

Hold "Ctrl" to enable pan & zoom

classDiagram
    class Dataset {
        +to_cog(...)
        +to_cog_bytes(...)
        +cog_info()
        +is_cog
        +validate_cog(...)
        +read_part(...) preview(...) point(...) read_tile(...)
    }
    class COG {
        <<engine>>
    }
    class write_cog {
        <<facade>>
    }
    Dataset --> COG : ds.cog
    write_cog ..> COG : delegates

Dataset.to_cog(...) is the single owner of write policy; the multi-input write_cog(...) facade (NumPy array / a labeled-array DataArray / gdal.Dataset / path / Dataset) normalises its input and delegates to it, so both produce identical output.

Public API at a glance

Concern	Symbol	Page
Write a COG (typed kwargs)	Dataset.to_cog	this page
Write from array/DataArray/path	write_cog	this page
Encode to in-memory bytes	Dataset.to_cog_bytes	this page
Named compression profiles	PROFILES · profile_options · validate_profile	this page
Creation-option helpers	merge_options · to_gdal_options · validate_blocksize	this page
Structured inspection	cog_info · COGInfo · OverviewLevel	Read & inspect
Validate	validate · ValidationReport · Dataset.is_cog	Read & inspect
Overview-decimated reads	read_part · preview · point · read_tile	Read & inspect
Command line	pyramids cog create\|validate\|info	CLI

See the COG cookbook for an end-to-end walkthrough and the COG basics notebook for a runnable, offline example.

Defaults

to_cog resolves the two pixel-affecting options per source dtype: the predictor (2 for integer, 3 for float) and the overview resampling (mode for categorical sources — integer dtype or a colour table — and average for continuous float). Both can be overridden explicitly.

Writing

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,
            )

to_cog(path, *, compression=None, overviews=None, tiling=None, bands=None, tags=None, layout=None, config=None, extra=None)

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(...).

Parameters:

Name	Type	Description	Default
path	str | Path	Destination path. Parent directory must exist.	required
compression	str | Compression | None	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.	None
overviews	Overviews | None	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).	None
tiling	Tiling | None	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.	None
bands	BandSelection | None	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.	None
tags	Tags | None	Metadata and colour table to stamp onto the output, as a :class:~pyramids.dataset.cog.Tags (band_tags, colormap, metadata).	None
layout	Layout | None	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.	None
config	dict[str, str] | None	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.	None
extra	Mapping[str, Any] | list[str] | None	Additional GDAL creation options as a mapping or legacy ['KEY=VALUE',...] list. Overrides conflicting group fields.	None

Returns:

Name	Type	Description
Path	Path	The resolved destination path.

Raises:

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

Warns:

Type	Description
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:

  >>> 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:

  >>> 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:

  >>> 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},
  ... )
  

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

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)

to_cog_bytes(**kwargs)

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.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Forwarded verbatim to :meth:to_cog (e.g. compression, layout, bands, extra). The same house defaults and dtype-aware resolution apply.	{}

Returns:

Name	Type	Description
bytes	bytes	The complete COG file contents.

Raises:

Type	Description
FailedToSaveError	GDAL failed to encode the COG.

Examples:

• Encode an in-memory Dataset to COG bytes and upload them:

  >>> 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
  

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

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

The write_cog facade

pyramids.dataset.cog.facade.write_cog(data, output, *, crs=None, transform=None, nodata=None, options=None, validate=True, strict=False)

Write raster data to disk as a Cloud Optimized GeoTIFF.

Thin convenience facade over :meth:pyramids.dataset.engines.cog.COG.to_cog. It accepts a wider range of inputs (NumPy array, a labeled DataArray, gdal.Dataset, path, or :class:~pyramids.dataset.Dataset), normalises them into a :class:~pyramids.dataset.Dataset, then delegates the entire write to to_cog — which owns all COG policy: the house defaults (DEFLATE, 512px tiles, BIGTIFF=IF_SAFER, NUM_THREADS=ALL_CPUS, embedded statistics), the dtype-aware predictor (2 for integer, 3 for float), the dtype-aware default overview resampling (mode for categorical, average for continuous), and the STATISTICS retry. Because policy lives in one place, write_cog and a direct ds.to_cog(...) produce identical output for identical input. Any caller-supplied options are forwarded as extra and override the defaults. By default the result is round-tripped through :func:pyramids.dataset.cog.validate.validate.

Parameters:

Name	Type	Description	Default
data	Any	Source raster. Accepted forms: :class:~pyramids.dataset.Dataset — used directly. :class:osgeo.gdal.Dataset — wrapped in a :class:Dataset. :class:numpy.ndarray — requires crs and transform. a labeled DataArray — geotransform is derived from the spatial coordinates; CRS from crs / DataArray metadata. :class:str / :class:~pathlib.Path — an existing raster.	required
output	str | Path	Destination path. The parent directory must exist.	required
crs	Any | None	CRS for the NumPy-array form (EPSG int, "EPSG:XXXX", WKT, or PROJ string). Also used as a fallback for DataArrays.	None
transform	tuple[float, ...] | None	6-tuple GDAL geotransform; required for the array form.	None
nodata	float | int | None	NoData scalar. For an array input it is passed to array construction; for a pre-built input it is applied before the write — on an in-memory copy when the source is read-only, so the source is never mutated.	None
options	CreationOptions | None	Caller overrides merged on top of :data:PYRAMIDS_COG_DEFAULTS. Keys are GDAL COG driver options (validated downstream). When PREDICTOR is absent it is auto-resolved from the raster dtype.	None
validate	bool	When True (default), validate the written file and raise :class:RuntimeError if it is not a valid COG.	True
strict	bool	Promote validation warnings to errors.	False

Returns:

Type	Description
Path	A (output_path, report) tuple. report is a
ValidationReport | None	class:~pyramids.dataset.cog.validate.ValidationReport when
tuple[Path, ValidationReport | None]	validate is True, otherwise None.

Raises:

Type	Description
ValueError	Required crs/transform missing for an array.
TypeError	data is an unsupported type.
RuntimeError	validate is True and the file failed COG validation.

Examples:

• Write a COG from a NumPy array (predictor auto-resolves to 3 for float):

  >>> import numpy as np  # doctest: +SKIP
  >>> arr = np.random.rand(256, 256).astype("float32")  # doctest: +SKIP
  >>> path, report = write_cog(  # doctest: +SKIP
  ...     arr, "out.tif", crs=4326,
  ...     transform=(0.0, 0.01, 0.0, 10.0, 0.0, -0.01),
  ... )
  >>> report.is_valid  # doctest: +SKIP
  True
  

• Re-encode an existing raster with overrides and skip validation:

  >>> path, report = write_cog(  # doctest: +SKIP
  ...     "plain.tif", "scene_cog.tif",
  ...     options={"COMPRESS": "ZSTD", "LEVEL": 18},
  ...     validate=False,
  ... )
  >>> report is None  # doctest: +SKIP
  True
  

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

def write_cog(
    data: Any,
    output: str | Path,
    *,
    crs: Any | None = None,
    transform: tuple[float, ...] | None = None,
    nodata: float | int | None = None,
    options: CreationOptions | None = None,
    validate: bool = True,
    strict: bool = False,
) -> tuple[Path, ValidationReport | None]:
    """Write raster data to disk as a Cloud Optimized GeoTIFF.

    Thin convenience facade over
    :meth:`pyramids.dataset.engines.cog.COG.to_cog`. It accepts a wider
    range of inputs (NumPy array, a labeled ``DataArray``, ``gdal.Dataset``,
    path, or :class:`~pyramids.dataset.Dataset`), normalises them into a
    :class:`~pyramids.dataset.Dataset`, then **delegates the entire write
    to ``to_cog``** — which owns all COG policy: the house defaults
    (DEFLATE, 512px tiles, ``BIGTIFF=IF_SAFER``, ``NUM_THREADS=ALL_CPUS``,
    embedded statistics), the dtype-aware predictor (``2`` for integer,
    ``3`` for float), the dtype-aware default overview resampling
    (``mode`` for categorical, ``average`` for continuous), and the
    ``STATISTICS`` retry. Because policy lives in one place, ``write_cog``
    and a direct ``ds.to_cog(...)`` produce identical output for identical
    input. Any caller-supplied ``options`` are forwarded as ``extra`` and
    override the defaults. By default the result is round-tripped through
    :func:`pyramids.dataset.cog.validate.validate`.

    Args:
        data: Source raster. Accepted forms:

            - :class:`~pyramids.dataset.Dataset` — used directly.
            - :class:`osgeo.gdal.Dataset` — wrapped in a :class:`Dataset`.
            - :class:`numpy.ndarray` — requires ``crs`` and ``transform``.
            - a labeled ``DataArray`` — geotransform is derived from the
              spatial coordinates; CRS from ``crs`` / DataArray metadata.
            - :class:`str` / :class:`~pathlib.Path` — an existing raster.
        output: Destination path. The parent directory must exist.
        crs: CRS for the NumPy-array form (EPSG int, ``"EPSG:XXXX"``, WKT,
            or PROJ string). Also used as a fallback for DataArrays.
        transform: 6-tuple GDAL geotransform; required for the array form.
        nodata: NoData scalar. For an array input it is passed to array
            construction; for a pre-built input it is applied before the write —
            on an in-memory copy when the source is read-only, so the source is
            never mutated.
        options: Caller overrides merged on top of
            :data:`PYRAMIDS_COG_DEFAULTS`. Keys are GDAL COG driver
            options (validated downstream). When ``PREDICTOR`` is absent it
            is auto-resolved from the raster dtype.
        validate: When ``True`` (default), validate the written file and
            raise :class:`RuntimeError` if it is not a valid COG.
        strict: Promote validation warnings to errors.

    Returns:
        A ``(output_path, report)`` tuple. ``report`` is a
        :class:`~pyramids.dataset.cog.validate.ValidationReport` when
        ``validate`` is ``True``, otherwise ``None``.

    Raises:
        ValueError: Required ``crs``/``transform`` missing for an array.
        TypeError: ``data`` is an unsupported type.
        RuntimeError: ``validate`` is ``True`` and the file failed COG
            validation.

    Examples:
        - Write a COG from a NumPy array (predictor auto-resolves to 3 for
          float):
            ```python
            >>> import numpy as np  # doctest: +SKIP
            >>> arr = np.random.rand(256, 256).astype("float32")  # doctest: +SKIP
            >>> path, report = write_cog(  # doctest: +SKIP
            ...     arr, "out.tif", crs=4326,
            ...     transform=(0.0, 0.01, 0.0, 10.0, 0.0, -0.01),
            ... )
            >>> report.is_valid  # doctest: +SKIP
            True

            ```
        - Re-encode an existing raster with overrides and skip validation:
            ```python
            >>> path, report = write_cog(  # doctest: +SKIP
            ...     "plain.tif", "scene_cog.tif",
            ...     options={"COMPRESS": "ZSTD", "LEVEL": 18},
            ...     validate=False,
            ... )
            >>> report is None  # doctest: +SKIP
            True

            ```
    """
    ds = _normalize_to_dataset(data, crs, transform, nodata)

    # Single write policy lives in COG.to_cog (ARC-1): house defaults, the
    # dtype-aware predictor (ARC-2), the category-safe default overview
    # resampling (ARC-3), and the STATISTICS retry (ARC-4) are all applied
    # there. write_cog only normalises the input and forwards the caller's
    # overrides as `extra`, so write_cog and a direct ds.to_cog(...) produce
    # identical output for identical input.
    output_path = ds.to_cog(output, extra=options or None)

    report: ValidationReport | None = None
    if validate:
        report = _validate_file(output_path, strict=strict)
        if not report.is_valid:
            raise RuntimeError(
                f"write_cog produced an invalid COG at {output_path}: {report.errors}"
            )
    return output_path, report

Creation options & profiles

pyramids.dataset.cog.options

COG creation-option types, serialization, and validation.

Provides the :data:CreationOptions alias (a Mapping[str, Any]), the named :data:PROFILES, and helpers used by :mod:pyramids.dataset.cog.write and :class:pyramids.dataset.engines.cog.COG:

• :func:to_gdal_options — serialize a mapping into GDAL's ['KEY=VALUE',...] list form.
• :func:merge_options — merge defaults with user-supplied extras (dict or legacy list[str]).
• :func:validate_blocksize — enforce the COG driver's power-of-2-in-[64, 4096] constraint.
• :func:validate_option_keys — gate unknown keys against :data:COG_DRIVER_OPTIONS.
• :func:profile_options / :func:validate_profile — named compression presets.

Each grouped option dataclass owns the logic that turns its own fields into the effective GDAL options via a _to_options method (Layout/Tiling serialize from their fields alone; Compression/Overviews also take the source band so the dtype-aware predictor / overview-resampling default can be resolved). BandSelection and Tags instead transform the source raster — _translate runs the in-memory band-subset/cast/NoData gdal.Translate and _stamp applies the colour table / metadata — so the write call site is a thin orchestrator over these methods rather than the home of the mapping logic.

These _-prefixed methods are private by convention, but they are the COG engine's internal contract: :class:pyramids.dataset.engines.cog.COG is their only intended caller (a sibling module in this subpackage), and they are not part of the public API. Rename or re-signature them in lockstep with that engine.

COG_DRIVER_OPTIONS = frozenset({'COMPRESS', 'LEVEL', 'QUALITY', 'NUM_THREADS', 'BLOCKSIZE', 'BIGTIFF', 'RESAMPLING', 'OVERVIEW_RESAMPLING', 'OVERVIEW_COUNT', 'OVERVIEW_COMPRESS', 'OVERVIEW_QUALITY', 'WARP_RESAMPLING', 'OVERVIEW_PREDICTOR', 'PREDICTOR', 'NBITS', 'TARGET_SRS', 'RES', 'EXTENT', 'ALIGNED_LEVELS', 'ADD_ALPHA', 'TILING_SCHEME', 'ZOOM_LEVEL', 'ZOOM_LEVEL_STRATEGY', 'MAX_Z_ERROR', 'STATISTICS', 'GEOTIFF_VERSION', 'SPARSE_OK', 'COPY_SRC_MDD', 'SRC_MDD'}) module-attribute

Whitelist of GDAL COG driver option keys (uppercased).

COG_READ_DEFAULTS = {'GDAL_DISABLE_READDIR_ON_OPEN': 'EMPTY_DIR', 'GDAL_HTTP_MERGE_CONSECUTIVE_RANGES': 'YES', 'GDAL_HTTP_MULTIRANGE': 'YES', 'GDAL_HTTP_MAX_RETRY': '3', 'GDAL_HTTP_RETRY_DELAY': '0.5', 'VSI_CACHE': 'TRUE'} module-attribute

GDAL config options that make remote /vsicurl/ COG reads efficient.

Without GDAL_DISABLE_READDIR_ON_OPEN a remote open issues a directory listing — often the single biggest latency hit, and the reason this preset exists. GDAL_HTTP_MULTIRANGE lets GDAL issue the scattered tile ranges a COG read produces as one multi-range request (GDAL_HTTP_MERGE_CONSECUTIVE_RANGES then coalesces the adjacent ones), and the retry pair rides out a transient 5xx from object storage instead of failing the whole read. That budget is deliberately larger than the OGC discovery one (:data:pyramids.base._ogc_api.GDAL_HTTP_MAX_RETRY): a COG read issues many range requests, so one flaky range should not lose the whole read, whereas a discovery pre-check is a single request in front of work that has not started yet.

Deliberately not included: CPL_VSIL_CURL_ALLOWED_EXTENSIONS. Restricting /vsicurl/ to .tif,.tiff makes GDAL refuse any URL whose path does not end in one — extensionless object keys, presigned S3 links carrying a query string, Planetary Computer and STAC asset hrefs — so a preset meant to speed reads up instead made them impossible. The readdir-skip above already delivers the latency win on its own. Applied by :func:pyramids.dataset.cog.validate.validate and :func:pyramids.dataset.cog.inspect.cog_info for remote paths when the caller passes no explicit config. Pure strings (no GDAL dependency here).

PROFILES = {'deflate': {'COMPRESS': 'DEFLATE', 'LEVEL': 9}, 'zstd': {'COMPRESS': 'ZSTD', 'LEVEL': 9}, 'lzw': {'COMPRESS': 'LZW'}, 'packbits': {'COMPRESS': 'PACKBITS'}, 'jpeg': {'COMPRESS': 'JPEG', 'QUALITY': 85}, 'webp': {'COMPRESS': 'WEBP', 'QUALITY': 85}, 'lerc': {'COMPRESS': 'LERC', 'MAX_Z_ERROR': 0.0}, 'lerc_deflate': {'COMPRESS': 'LERC_DEFLATE', 'MAX_Z_ERROR': 0.0}, 'lerc_zstd': {'COMPRESS': 'LERC_ZSTD', 'MAX_Z_ERROR': 0.0}, 'raw': {'COMPRESS': 'NONE'}} module-attribute

Named COG creation profiles (compression presets).

Each profile seeds only the compression-related options; the predictor and overview resampling are still resolved per-dtype by :meth:pyramids.dataset.engines.cog.COG.to_cog, and explicit kwargs / extra override the profile. Limited to keys in :data:COG_DRIVER_OPTIONS.

profile_options(name)

Return a copy of the named profile's creation options.

Parameters:

Name	Type	Description	Default
name	str	Profile name (case-insensitive), e.g. "deflate", "zstd", "jpeg".	required

Returns:

Type	Description
dict[str, Any]	A new dict of the profile's options.

Raises:

Type	Description
ValueError	When name is not a known profile.

Examples:

• Look up the zstd preset:

  >>> profile_options("zstd")
  {'COMPRESS': 'ZSTD', 'LEVEL': 9}
  

• Names are case-insensitive:

  >>> profile_options("LZW")
  {'COMPRESS': 'LZW'}
  

• Unknown names are rejected:

  >>> profile_options("bogus")  # doctest: +IGNORE_EXCEPTION_DETAIL
  Traceback (most recent call last):
  ...
  ValueError: unknown COG profile 'bogus'...
  

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

def profile_options(name: str) -> dict[str, Any]:
    """Return a copy of the named profile's creation options.

    Args:
        name: Profile name (case-insensitive), e.g. ``"deflate"``, ``"zstd"``,
            ``"jpeg"``.

    Returns:
        A new dict of the profile's options.

    Raises:
        ValueError: When ``name`` is not a known profile.

    Examples:
        - Look up the zstd preset:
            ```python
            >>> profile_options("zstd")
            {'COMPRESS': 'ZSTD', 'LEVEL': 9}

            ```
        - Names are case-insensitive:
            ```python
            >>> profile_options("LZW")
            {'COMPRESS': 'LZW'}

            ```
        - Unknown names are rejected:
            ```python
            >>> profile_options("bogus")  # doctest: +IGNORE_EXCEPTION_DETAIL
            Traceback (most recent call last):
            ...
            ValueError: unknown COG profile 'bogus'...

            ```
    """
    key = name.lower()
    if key not in PROFILES:
        raise ValueError(
            f"unknown COG profile {name!r}; choose from {sorted(PROFILES)}"
        )
    return dict(PROFILES[key])

validate_profile(name, dtype_name, band_count)

Raise :class:ValueError if a source violates a profile's constraints.

Some profiles only accept specific dtypes / band counts (JPEG: Byte with 1-3 bands; WEBP: Byte with 3-4 bands). Other profiles are unconstrained and pass silently.

Parameters:

Name	Type	Description	Default
name	str	Profile name (case-insensitive).	required
dtype_name	str	GDAL dtype name of the source (e.g. "Byte", "Float32").	required
band_count	int	Number of bands in the source.	required

Raises:

Type	Description
ValueError	When the source dtype or band count is incompatible.

Examples:

• An unconstrained profile always passes:

  >>> validate_profile("deflate", "Float32", 4)
  

• JPEG requires Byte with <= 3 bands:

  >>> validate_profile("jpeg", "Float32", 1)  # doctest: +IGNORE_EXCEPTION_DETAIL
  Traceback (most recent call last):
  ...
  ValueError: jpeg profile requires dtype in...
  

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

def validate_profile(name: str, dtype_name: str, band_count: int) -> None:
    """Raise :class:`ValueError` if a source violates a profile's constraints.

    Some profiles only accept specific dtypes / band counts (JPEG: ``Byte``
    with 1-3 bands; WEBP: ``Byte`` with 3-4 bands). Other profiles are
    unconstrained and pass silently.

    Args:
        name: Profile name (case-insensitive).
        dtype_name: GDAL dtype name of the source (e.g. ``"Byte"``,
            ``"Float32"``).
        band_count: Number of bands in the source.

    Raises:
        ValueError: When the source dtype or band count is incompatible.

    Examples:
        - An unconstrained profile always passes:
            ```python
            >>> validate_profile("deflate", "Float32", 4)

            ```
        - JPEG requires Byte with <= 3 bands:
            ```python
            >>> validate_profile("jpeg", "Float32", 1)  # doctest: +IGNORE_EXCEPTION_DETAIL
            Traceback (most recent call last):
            ...
            ValueError: jpeg profile requires dtype in...

            ```
    """
    key = name.lower()
    if key not in _PROFILE_DTYPE_CONSTRAINTS:
        return
    allowed_dtypes, (min_bands, max_bands) = _PROFILE_DTYPE_CONSTRAINTS[key]
    if dtype_name not in allowed_dtypes:
        raise ValueError(
            f"{key} profile requires dtype in {sorted(allowed_dtypes)}; "
            f"got {dtype_name}. Use a different profile (e.g. 'deflate')."
        )
    if not (min_bands <= band_count <= max_bands):
        raise ValueError(
            f"{key} profile requires {min_bands}-{max_bands} bands; got {band_count}."
        )

to_gdal_options(opts)

Serialize a mapping into GDAL's ['KEY=VALUE',...] list form.

Keys are uppercased; values are stringified via :func:_stringify (booleans become "YES"/"NO"). None values are skipped so callers can pass optional kwargs through unchanged.

Parameters:

Name	Type	Description	Default
opts	CreationOptions | None	Mapping of option names to values, or None.	required

Returns:

Type	Description
list[str]	List of "KEY=VALUE" strings. Empty list when opts is None.

Examples:

• Serialize a compression config:

  >>> to_gdal_options({"COMPRESS": "DEFLATE", "LEVEL": 9})
  ['COMPRESS=DEFLATE', 'LEVEL=9']
  

• Booleans become GDAL's YES/NO convention:

  >>> to_gdal_options({"STATISTICS": True, "SPARSE_OK": False})
  ['STATISTICS=YES', 'SPARSE_OK=NO']
  

• None values are dropped so optional kwargs flow through unchanged:

  >>> to_gdal_options({"COMPRESS": "LZW", "LEVEL": None})
  ['COMPRESS=LZW']
  >>> to_gdal_options(None)
  []
  

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

def to_gdal_options(opts: CreationOptions | None) -> list[str]:
    """Serialize a mapping into GDAL's `['KEY=VALUE',...]` list form.

    Keys are uppercased; values are stringified via :func:`_stringify`
    (booleans become `"YES"`/`"NO"`). `None` values are skipped so
    callers can pass optional kwargs through unchanged.

    Args:
        opts: Mapping of option names to values, or `None`.

    Returns:
        List of `"KEY=VALUE"` strings. Empty list when `opts` is `None`.

    Examples:
        - Serialize a compression config:
            ```python
            >>> to_gdal_options({"COMPRESS": "DEFLATE", "LEVEL": 9})
            ['COMPRESS=DEFLATE', 'LEVEL=9']

            ```
        - Booleans become GDAL's YES/NO convention:
            ```python
            >>> to_gdal_options({"STATISTICS": True, "SPARSE_OK": False})
            ['STATISTICS=YES', 'SPARSE_OK=NO']

            ```
        - None values are dropped so optional kwargs flow through unchanged:
            ```python
            >>> to_gdal_options({"COMPRESS": "LZW", "LEVEL": None})
            ['COMPRESS=LZW']
            >>> to_gdal_options(None)
            []

            ```
    """
    result: list[str]
    if opts is None:
        result = []
    else:
        result = [
            f"{str(k).upper()}={_stringify(v)}"
            for k, v in opts.items()
            if v is not None
        ]
    return result

merge_options(defaults, extra)

Merge default options with user-supplied extras; extras win.

Accepts extra as either a mapping {'KEY': value} or the legacy list form ['KEY=VALUE',...] used by :meth:pyramids.dataset.Dataset.to_file. All keys in the returned dict are uppercased; None values from either source are dropped.

Parameters:

Name	Type	Description	Default
defaults	CreationOptions	Baseline options (typically derived from kwargs in :meth:pyramids.dataset.engines.COG.to_cog).	required
extra	CreationOptions | list[str] | None	User-provided overrides as a mapping, list[str], or None.	required

Returns:

Name	Type	Description
New	dict[str, Any]	class:dict with all keys uppercased and None values
	dict[str, Any]	removed; extra entries override defaults on conflict.

Raises:

Type	Description
ValueError	When a legacy list-form entry lacks =.

Examples:

• Dict extras override defaults on conflict:

  >>> merge_options({"COMPRESS": "DEFLATE"}, {"COMPRESS": "ZSTD"})
  {'COMPRESS': 'ZSTD'}
  

• Legacy list-of-string form is also accepted for back-compat:

  >>> merge_options({"COMPRESS": "DEFLATE"}, ["LEVEL=9"])
  {'COMPRESS': 'DEFLATE', 'LEVEL': '9'}
  

• None extras returns a copy of the defaults:

  >>> merge_options({"COMPRESS": "DEFLATE"}, None)
  {'COMPRESS': 'DEFLATE'}
  

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

def merge_options(
    defaults: CreationOptions,
    extra: CreationOptions | list[str] | None,
) -> dict[str, Any]:
    """Merge default options with user-supplied extras; extras win.

    Accepts `extra` as either a mapping `{'KEY': value}` or the legacy
    list form `['KEY=VALUE',...]` used by
    :meth:`pyramids.dataset.Dataset.to_file`. All keys in the returned
    dict are uppercased; `None` values from either source are dropped.

    Args:
        defaults: Baseline options (typically derived from kwargs in
            :meth:`pyramids.dataset.engines.COG.to_cog`).
        extra: User-provided overrides as a mapping, `list[str]`, or
            `None`.

    Returns:
        New :class:`dict` with all keys uppercased and `None` values
        removed; `extra` entries override `defaults` on conflict.

    Raises:
        ValueError: When a legacy list-form entry lacks `=`.

    Examples:
        - Dict extras override defaults on conflict:
            ```python
            >>> merge_options({"COMPRESS": "DEFLATE"}, {"COMPRESS": "ZSTD"})
            {'COMPRESS': 'ZSTD'}

            ```
        - Legacy list-of-string form is also accepted for back-compat:
            ```python
            >>> merge_options({"COMPRESS": "DEFLATE"}, ["LEVEL=9"])
            {'COMPRESS': 'DEFLATE', 'LEVEL': '9'}

            ```
        - None extras returns a copy of the defaults:
            ```python
            >>> merge_options({"COMPRESS": "DEFLATE"}, None)
            {'COMPRESS': 'DEFLATE'}

            ```
    """
    merged: dict[str, Any] = {
        str(k).upper(): v for k, v in defaults.items() if v is not None
    }
    # `extra is None` -> nothing to merge, keep the defaults.
    if isinstance(extra, list):
        merged.update(_parse_list_extra(extra))
    elif extra is not None:
        merged.update({str(k).upper(): v for k, v in extra.items() if v is not None})
    return merged

validate_blocksize(value)

Raise :class:ValueError if value is not a valid COG tile size.

The GDAL COG driver requires BLOCKSIZE to be a power of 2 in the closed range [64, 4096].

Parameters:

Name	Type	Description	Default
value	int	Proposed blocksize.	required

Raises:

Type	Description
ValueError	If value is outside the allowed set.

Examples:

• Valid power-of-2 blocksizes return silently:

  >>> validate_blocksize(512)
  >>> validate_blocksize(256)
  

• Non-power-of-2 is rejected:

  >>> validate_blocksize(500) # doctest: +IGNORE_EXCEPTION_DETAIL
  Traceback (most recent call last):
  ...
  ValueError: blocksize must be a power of 2 in [64, 4096]; got 500...
  

• Out-of-range values are rejected:

  >>> validate_blocksize(32) # doctest: +IGNORE_EXCEPTION_DETAIL
  Traceback (most recent call last):
  ...
  ValueError: blocksize must be a power of 2 in [64, 4096]; got 32...
  

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

def validate_blocksize(value: int) -> None:
    """Raise :class:`ValueError` if `value` is not a valid COG tile size.

    The GDAL COG driver requires `BLOCKSIZE` to be a power of 2 in
    the closed range [64, 4096].

    Args:
        value: Proposed blocksize.

    Raises:
        ValueError: If `value` is outside the allowed set.

    Examples:
        - Valid power-of-2 blocksizes return silently:
            ```python
            >>> validate_blocksize(512)
            >>> validate_blocksize(256)

            ```
        - Non-power-of-2 is rejected:
            ```python
            >>> validate_blocksize(500) # doctest: +IGNORE_EXCEPTION_DETAIL
            Traceback (most recent call last):
            ...
            ValueError: blocksize must be a power of 2 in [64, 4096]; got 500...

            ```
        - Out-of-range values are rejected:
            ```python
            >>> validate_blocksize(32) # doctest: +IGNORE_EXCEPTION_DETAIL
            Traceback (most recent call last):
            ...
            ValueError: blocksize must be a power of 2 in [64, 4096]; got 32...

            ```
    """
    if value not in _VALID_BLOCKSIZES:
        raise ValueError(
            f"blocksize must be a power of 2 in [64, 4096]; got {value}. "
            f"Valid values: {sorted(_VALID_BLOCKSIZES)}"
        )

validate_option_keys(opts)

Raise :class:ValueError for any key not in :data:COG_DRIVER_OPTIONS.

Keys are compared case-insensitively.

Parameters:

Name	Type	Description	Default
opts	CreationOptions	Mapping of option names to values.	required

Raises:

Type	Description
ValueError	If any key is not a recognized COG driver option.

Examples:

• Known keys return silently:

  >>> validate_option_keys({"COMPRESS": "DEFLATE"})
  >>> validate_option_keys({"BLOCKSIZE": 512, "BIGTIFF": "IF_SAFER"})
  

• Unknown keys raise ValueError naming the offender:

  >>> validate_option_keys({"NONSENSE": "x"}) # doctest: +IGNORE_EXCEPTION_DETAIL
  Traceback (most recent call last):
  ...
  ValueError: Unknown COG driver option(s): ['NONSENSE']...
  

• Empty mapping is accepted:

  >>> validate_option_keys({})
  

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

def validate_option_keys(opts: CreationOptions) -> None:
    """Raise :class:`ValueError` for any key not in :data:`COG_DRIVER_OPTIONS`.

    Keys are compared case-insensitively.

    Args:
        opts: Mapping of option names to values.

    Raises:
        ValueError: If any key is not a recognized COG driver option.

    Examples:
        - Known keys return silently:
            ```python
            >>> validate_option_keys({"COMPRESS": "DEFLATE"})
            >>> validate_option_keys({"BLOCKSIZE": 512, "BIGTIFF": "IF_SAFER"})

            ```
        - Unknown keys raise ValueError naming the offender:
            ```python
            >>> validate_option_keys({"NONSENSE": "x"}) # doctest: +IGNORE_EXCEPTION_DETAIL
            Traceback (most recent call last):
            ...
            ValueError: Unknown COG driver option(s): ['NONSENSE']...

            ```
        - Empty mapping is accepted:
            ```python
            >>> validate_option_keys({})

            ```
    """
    unknown = {str(k).upper() for k in opts.keys()} - COG_DRIVER_OPTIONS
    if unknown:
        raise ValueError(
            f"Unknown COG driver option(s): {sorted(unknown)}. "
            f"Valid options: {sorted(COG_DRIVER_OPTIONS)}"
        )