Reading, inspecting & validating

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

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

Structured inspection

pyramids.dataset.cog.inspect

Structured Cloud Optimized GeoTIFF inspection.

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

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

OverviewLevel dataclass

One level of a COG's internal overview pyramid.

Attributes:

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

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

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

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

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

COGInfo dataclass

Structured metadata describing a (Cloud Optimized) GeoTIFF.

Attributes:

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

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

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

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

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

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

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

cog_info(path, config=None)

Inspect a raster and return its structured COG metadata.

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

Parameters:

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

Returns:

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

Raises:

Type	Description
FileNotFoundError	When path cannot be opened by GDAL.

Examples:

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

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

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

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

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

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

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

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

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

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

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

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

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

_cog_info_impl(p)

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

Parameters:

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

Returns:

Name	Type	Description
COGInfo	COGInfo	The structured metadata.

Raises:

Type	Description
FileNotFoundError	When p cannot be opened by GDAL.

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

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

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

    Returns:
        COGInfo: The structured metadata.

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

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

        gt = ds.GetGeoTransform()
        min_x, max_y = gt[0], gt[3]
        max_x = min_x + gt[1] * width
        min_y = max_y + gt[5] * height

        srs = ds.GetSpatialRef()
        epsg: int | None = None
        if srs is not None:
            code = srs.GetAuthorityCode(None)
            epsg = int(code) if code is not None else None

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

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

Overview-decimated reads

pyramids.dataset.engines.cog.COG

Bases: _Engine

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):
    """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,
        *,
        profile: str | None = None,
        compress: str | None = None,
        level: int | None = None,
        quality: int | None = None,
        blocksize: int = 512,
        predictor: str | int | None = None,
        bigtiff: str = "IF_SAFER",
        num_threads: int | str = "ALL_CPUS",
        overview_resampling: str | None = None,
        overview_count: int | None = None,
        overview_compress: str | None = None,
        tiling_scheme: str | None = None,
        zoom_level: int | None = None,
        zoom_level_strategy: str = "auto",
        aligned_levels: int | None = None,
        resampling: str = "nearest",
        add_mask: bool = False,
        sparse_ok: bool = False,
        target_srs: int | str | None = None,
        statistics: bool = True,
        indexes: list[int] | None = None,
        out_dtype: str | None = None,
        nodata: float | int | None = None,
        band_tags: dict[int, dict[str, Any]] | None = None,
        colormap: dict[int, tuple[int, int, int, int]] | None = None,
        metadata: dict[str, Any] | 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.

        Args:
            path: Destination path. Parent directory must exist.
            profile: Named compression preset (case-insensitive) — one of
                `deflate`, `zstd`, `lzw`, `packbits`, `jpeg`, `webp`,
                `lerc`, `lerc_deflate`, `lerc_zstd`, `raw`. Seeds the
                compression options; explicit `compress`/`level`/`quality`
                and `extra` override it. `jpeg`/`webp` enforce dtype/band
                constraints (Byte; 1-3 / 3-4 bands).
            compress: Compression method. `DEFLATE`, `LZW`, and
                `NONE` are guaranteed by every GDAL build. `JPEG`
                is almost always available. `ZSTD`, `WEBP`,
                `LERC`, `LERC_DEFLATE`, and `LERC_ZSTD` require
                the GDAL build to have been compiled with the
                corresponding library (libzstd / libwebp / LERC); on
                a GDAL build lacking them, the COG driver will raise
                at write time. To probe what your GDAL supports:

                ```python
                from osgeo import gdal
                meta = gdal.GetDriverByName("GTiff").GetMetadataItem(
                    "DMD_CREATIONOPTIONLIST"
                )
                print("ZSTD" in (meta or ""))
                ```
            level: Compression level (e.g., 1-12 for DEFLATE, 1-22 ZSTD).
            quality: Lossy-compression quality 1-100 (JPEG/WEBP).
            blocksize: Internal tile size; power of 2 in [64, 4096].
            predictor: `"YES"`/`"STANDARD"`/`"FLOATING_POINT"` or 1/2/3.
                Defaults to `None`, which auto-resolves per the source
                dtype: `2` (horizontal differencing) for integer rasters,
                `3` (floating-point predictor) for float rasters. Pass an
                explicit value to override.
            bigtiff: `"IF_SAFER"` (default), `"YES"`, `"NO"`,
                `"IF_NEEDED"`.
            num_threads: Worker threads; `"ALL_CPUS"` or an int.
            overview_resampling: `nearest`, `average`, `bilinear`,
                `cubic`, `cubicspline`, `lanczos`, `mode`,
                `rms`, `gauss`. Defaults to `None`, which auto-resolves
                per the source dtype: `mode` for categorical sources
                (integer dtype or a colour table) and `average` for
                continuous (float) sources. The categorical guardrail
                warns only when *you* explicitly pass an averaging method
                on categorical data — never for this auto-resolved default.
            overview_count: Number of overview levels (default: auto).
            overview_compress: Compression for overview IFDs.
            tiling_scheme: e.g., `"GoogleMapsCompatible"` for a
                web-optimized COG (EPSG:3857).
            zoom_level: Advanced tiling-scheme knob: pin the maximum zoom level.
            zoom_level_strategy: Advanced tiling-scheme knob: `auto` (default),
                `lower`, or `upper` zoom-level selection.
            aligned_levels: Advanced tiling-scheme knob: number of overview
                levels aligned to the tiling scheme.
            resampling: Warp resampling when `tiling_scheme` or
                `target_srs` reprojects.
            add_mask: Add an alpha band for transparency.
            sparse_ok: Allow sparse (unfilled) tiles.
            target_srs: Reproject before write. Int for EPSG or a WKT
                / PROJ string.
            statistics: Compute and embed band statistics.
            indexes: 0-based band indices to keep, in order (e.g. `[3, 2, 1]`
                to select and reorder bands). `None` keeps all bands. When
                set, the source is pre-processed through an in-memory
                `gdal.Translate` before the COG write.
            out_dtype: Output NumPy dtype name to cast to (e.g. `"uint8"`,
                `"int16"`). `None` keeps the source dtype. The dtype-aware
                predictor is resolved from the *post-cast* dtype.
            nodata: NoData value to set on the output. `None` keeps the
                source NoData.
            band_tags: Per-band metadata to stamp onto the output, keyed by
                0-based band index, e.g. `{0: {"name": "NDVI"}}`. Useful when
                the source is a bare array/DataArray that carries no band
                descriptions.
            colormap: Palette to attach to band 1, mapping pixel value to an
                `(R, G, B, A)` tuple, e.g. `{0: (0, 0, 0, 255), 1: (255, 0, 0, 255)}`.
            metadata: Dataset-level metadata items to stamp onto the output.
            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 kwargs.

        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 has a color table) and `overview_resampling`
                is an averaging method.

        Note:
            Setting `tiling_scheme` (e.g., `GoogleMapsCompatible`)
            implies a specific SRS — `target_srs` is ignored in that
            case. A `UserWarning` is emitted if both are provided.

        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
                >>> import numpy as np  # doctest: +SKIP
                >>> from pyramids.dataset import Dataset  # doctest: +SKIP
                >>> arr = np.random.rand(256, 256).astype("float32")  # doctest: +SKIP
                >>> ds = Dataset.create_from_array(  # doctest: +SKIP
                ...     arr, top_left_corner=(0, 0), cell_size=0.001, epsg=4326,
                ... )
                >>> out = ds.to_cog("out.tif", compress="ZSTD")  # doctest: +SKIP
                >>> out.name  # doctest: +SKIP
                'out.tif'

                ```
            - Produce a web-optimized COG for a tile server:
                ```python
                >>> web = ds.to_cog("web.tif", tiling_scheme="GoogleMapsCompatible")  # doctest: +SKIP
                >>> reopened = Dataset.read_file(web)  # doctest: +SKIP
                >>> reopened.epsg  # doctest: +SKIP
                3857

                ```
            - Forward additional GDAL options through `extra`:
                ```python
                >>> _ = ds.to_cog(  # doctest: +SKIP
                ...     "precise.tif",
                ...     compress="LERC",
                ...     extra={"MAX_Z_ERROR": 0.001},
                ... )

                ```
        """
        validate_blocksize(blocksize)
        if tiling_scheme is not None and target_srs is not None:
            warnings.warn(
                "Both tiling_scheme and target_srs provided; "
                "tiling_scheme wins and target_srs is ignored.",
                UserWarning,
                stacklevel=2,
            )
            target_srs = None

        # Build the effective source (PB-4): when band-subsetting, casting the
        # dtype, or (re)setting NoData, pre-process through an in-memory
        # gdal.Translate so the predictor/overview policy below — and the COG
        # write itself — see the *output* bands, not the original source.
        source_ds, source_band0 = self._effective_source(
            indexes, out_dtype, nodata, band_tags, colormap, metadata
        )

        # Resolve a named profile (PB-5): it seeds the compression options;
        # explicit kwargs and `extra` override it. jpeg/webp enforce dtype/band
        # constraints against the *effective* source.
        profile_opts: dict[str, Any] = {}
        if profile is not None:
            validate_profile(
                profile,
                gdal.GetDataTypeName(source_band0.DataType),
                source_ds.RasterCount,
            )
            profile_opts = profile_options(profile)
        eff_compress = (
            compress
            if compress is not None
            else profile_opts.get("COMPRESS", "DEFLATE")
        )
        eff_level = level if level is not None else profile_opts.get("LEVEL")
        eff_quality = quality if quality is not None else profile_opts.get("QUALITY")
        profile_extra = {
            k: v
            for k, v in profile_opts.items()
            if k not in ("COMPRESS", "LEVEL", "QUALITY")
        }

        # Single house policy lives here (ARC-1): `to_cog` resolves the
        # dtype-dependent defaults so a direct `ds.to_cog(...)` and the
        # `write_cog(...)` facade — which now just delegates here — produce
        # identical output for identical input.
        if predictor is None:
            # Per-dtype predictor (ARC-2): 2 for integer, 3 for float. GeoTIFF
            # bands share a dtype, so band 0 decides for the whole file. Pass an
            # explicit `predictor=` to override for an (atypical) mixed source.
            predictor = resolve_cog_predictor(source_band0.DataType)
        caller_chose_resampling = overview_resampling is not None
        if overview_resampling is None:
            # Category-safe default (ARC-3): `mode` for integer/colour-table
            # sources, `average` for continuous. Chosen so the default never
            # corrupts categorical rasters and never trips the guardrail below.
            overview_resampling = default_cog_overview_resampling(
                source_band0.DataType, source_band0.GetColorTable() is not None
            )
        if caller_chose_resampling:
            # Only warn when the *caller* explicitly asked for an averaging
            # resampler on categorical data — never for a default we picked.
            self._warn_if_categorical_with_averaging(
                overview_resampling, band=source_band0
            )

        num_threads_str = (
            num_threads if isinstance(num_threads, str) else str(num_threads)
        )
        defaults: dict[str, Any] = {
            "COMPRESS": eff_compress,
            "LEVEL": eff_level,
            "QUALITY": eff_quality,
            **profile_extra,
            "BLOCKSIZE": blocksize,
            "PREDICTOR": predictor,
            "BIGTIFF": bigtiff,
            "NUM_THREADS": num_threads_str,
            "OVERVIEW_RESAMPLING": overview_resampling,
            "OVERVIEW_COUNT": overview_count,
            "OVERVIEW_COMPRESS": overview_compress,
            "TILING_SCHEME": tiling_scheme,
            "ZOOM_LEVEL": zoom_level,
            "ZOOM_LEVEL_STRATEGY": zoom_level_strategy,
            "ALIGNED_LEVELS": aligned_levels,
            "WARP_RESAMPLING": (resampling if (tiling_scheme or target_srs) else None),
            "ADD_ALPHA": True if add_mask else None,
            "SPARSE_OK": True if sparse_ok else None,
            "STATISTICS": "YES" if statistics else None,
        }
        if target_srs is not None:
            defaults["TARGET_SRS"] = (
                f"EPSG:{target_srs}" if isinstance(target_srs, int) else target_srs
            )

        options = merge_options(defaults, extra)
        with config_context(config):
            self._translate_with_statistics_retry(path, options, src=source_ds)
        return Path(path)

    def _effective_source(
        self,
        indexes: list[int] | None,
        out_dtype: str | None,
        nodata: float | int | None,
        band_tags: dict[int, dict[str, Any]] | None = None,
        colormap: dict[int, tuple[int, int, int, int]] | None = None,
        metadata: dict[str, Any] | None = None,
    ) -> tuple[gdal.Dataset, Any]:
        """Return the source dataset (optionally pre-processed) and its band 0.

        When band-subsetting / dtype-casting / setting NoData (PB-4) the backing
        raster is run through an in-memory ``gdal.Translate``; when stamping
        band tags / colourmap / metadata (PC-2) it is copied to a MEM dataset
        first so the user's open dataset is **never mutated**. With neither, the
        backing raster is returned unchanged.

        Args:
            indexes: 0-based band indices to keep/reorder, or ``None``.
            out_dtype: Output NumPy dtype name to cast to, or ``None``.
            nodata: NoData value to set, or ``None``.
            band_tags: Per-band metadata keyed by 0-based band index, or ``None``.
            colormap: Palette for band 1 (value -> RGBA), or ``None``.
            metadata: Dataset-level metadata, or ``None``.

        Returns:
            A ``(dataset, band1)`` tuple where ``band1`` is GDAL band 1 of the
            returned dataset (used for predictor/resampling resolution).
        """
        needs_translate = (
            indexes is not None or out_dtype is not None or nodata is not None
        )
        needs_stamp = bool(band_tags or colormap or metadata)
        if not needs_translate and not needs_stamp:
            ds = self._ds._raster
            return ds, ds.GetRasterBand(1)

        if needs_translate:
            translate_kwargs: dict[str, Any] = {}
            if indexes is not None:
                # pyramids band indices are 0-based; GDAL bandList is 1-based.
                translate_kwargs["bandList"] = [i + 1 for i in indexes]
            if out_dtype is not None:
                translate_kwargs["outputType"] = numpy_to_gdal_dtype(out_dtype)
            if nodata is not None:
                translate_kwargs["noData"] = nodata
            mem = gdal.Translate(
                "", self._ds._raster, format="MEM", **translate_kwargs
            )
        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 build the pre-processed COG source "
                f"(indexes={indexes}, out_dtype={out_dtype}, nodata={nodata})"
            )
        if needs_stamp:
            self._stamp_metadata(mem, band_tags, colormap, metadata)
        return mem, mem.GetRasterBand(1)

    @staticmethod
    def _stamp_metadata(
        ds: gdal.Dataset,
        band_tags: dict[int, dict[str, Any]] | None,
        colormap: dict[int, tuple[int, int, int, int]] | None,
        metadata: dict[str, Any] | None,
    ) -> None:
        """Stamp band tags / colourmap / dataset metadata onto a MEM dataset.

        Args:
            ds: The (copied) dataset to mutate.
            band_tags: Per-band metadata keyed by 0-based band index.
            colormap: Palette for band 1 (value -> RGBA tuple). GeoTIFF only
                supports a colour table on a single-band `Byte` / `UInt16`
                raster; a `ValueError` is raised up-front for other dtypes
                (GDAL would otherwise fail deep in `CreateCopy`).
            metadata: Dataset-level metadata items.

        Raises:
            ValueError: When `colormap` is applied to a band whose dtype is
                not `Byte`/`UInt16`.
        """
        if metadata:
            ds.SetMetadata({str(k): str(v) for k, v in metadata.items()})
        if colormap:
            band = ds.GetRasterBand(1)
            if band.DataType not in _PALETTE_GDAL_DTYPES:
                raise ValueError(
                    f"colormap is only supported on Byte/UInt16 rasters; got "
                    f"{gdal.GetDataTypeName(band.DataType)}. Cast first with "
                    f"to_cog(..., out_dtype='uint8'), or drop the colormap."
                )
            color_table = gdal.ColorTable()
            for value, rgba in colormap.items():
                color_table.SetColorEntry(int(value), tuple(rgba))
            band.SetColorTable(color_table)
            band.SetColorInterpretation(gdal.GCI_PaletteIndex)
        if band_tags:
            for index, tags in band_tags.items():
                # 0-based index -> GDAL 1-based band number.
                ds.GetRasterBand(index + 1).SetMetadata(
                    {str(k): str(v) for k, v in tags.items()}
                )

    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. ``compress``,
                ``blocksize``, ``predictor``, ``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(compress="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 = f"/vsimem/{uuid.uuid4().hex}.tif"
        try:
            self.to_cog(vsi_path, **kwargs)
            handle = gdal.VSIFOpenL(vsi_path, "rb")
            if handle is None:
                raise FailedToSaveError(
                    f"could not reopen in-memory COG at {vsi_path}"
                )
            try:
                gdal.VSIFSeekL(handle, 0, 2)  # SEEK_END
                size = gdal.VSIFTellL(handle)
                gdal.VSIFSeekL(handle, 0, 0)  # SEEK_SET
                data = gdal.VSIFReadL(1, size, handle)
            finally:
                gdal.VSIFCloseL(handle)
        finally:
            gdal.Unlink(vsi_path)
        return bytes(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).
        """
        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
    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

    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)

    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.
        """
        fn = self._ds.file_name
        if not fn or fn.startswith("/vsimem/"):
            return None
        return fn

    def read_part(
        self,
        bbox: tuple[float, float, float, float],
        *,
        dst_width: int | None = None,
        dst_height: int | None = None,
        bbox_crs: int = 4326,
        resampling: str = "bilinear",
        band: int | None = None,
    ) -> np.ndarray:
        """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.

        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: EPSG code of `bbox`. Reprojected to the dataset CRS
                when different. Defaults to 4326 (WGS84 lon/lat).
            resampling: One of `nearest`, `bilinear`, `cubic`,
                `cubicspline`, `lanczos`, `average`, `mode`.
            band: 0-based band index. `None` reads all bands.

        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.

        Raises:
            ValueError: Unknown `resampling`.
            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`.

        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)

                ```
        """
        if resampling not in _RESAMPLING_ALG:
            raise ValueError(
                f"unknown resampling {resampling!r}; "
                f"choose from {sorted(_RESAMPLING_ALG)}"
            )
        ds = self._ds._raster
        min_x, min_y, max_x, max_y = self._reproject_bbox(bbox, bbox_crs)
        inv = gdal.InvGeoTransform(ds.GetGeoTransform())
        px_tl, py_tl = gdal.ApplyGeoTransform(inv, min_x, max_y)
        px_br, py_br = gdal.ApplyGeoTransform(inv, 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
        alg = _RESAMPLING_ALG[resampling]
        source = ds if band is None else ds.GetRasterBand(band + 1)

        fully_inside = (
            ix0 == req_xoff
            and iy0 == req_yoff
            and ix1 == req_xoff + req_xsize
            and iy1 == req_yoff + req_ysize
        )
        if fully_inside:
            return np.asarray(
                source.ReadAsArray(
                    ix0,
                    iy0,
                    ix1 - ix0,
                    iy1 - iy0,
                    buf_xsize=out_w,
                    buf_ysize=out_h,
                    resample_alg=alg,
                )
            )

        # 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
        return out

    @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 nodata
        return 0 if is_integer_gdal_dtype(band.DataType) else float("nan")

    def preview(
        self,
        *,
        max_size: int = 1024,
        resampling: str = "bilinear",
        band: int | None = None,
    ) -> np.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.

        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)`. Pixel values only — no transform, bounds,
            or CRS is attached to the returned array.

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

                ```
        """
        if resampling not in _RESAMPLING_ALG:
            raise ValueError(
                f"unknown resampling {resampling!r}; "
                f"choose from {sorted(_RESAMPLING_ALG)}"
            )
        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))
        alg = _RESAMPLING_ALG[resampling]
        ds = self._ds._raster
        source = ds if band is None else ds.GetRasterBand(band + 1)
        return np.asarray(
            source.ReadAsArray(buf_xsize=out_w, buf_ysize=out_h, resample_alg=alg)
        )

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

        Args:
            x: X / longitude / easting in `point_crs`.
            y: Y / latitude / northing in `point_crs`.
            point_crs: EPSG code of `(x, y)`. Reprojected to the dataset CRS
                when different. Defaults to 4326.
            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`. Pixel values only — no
            coordinate metadata is attached.

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

                ```
        """
        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)
        arr = np.asarray(source.ReadAsArray(col, row, 1, 1))
        return arr.reshape(-1) if band is None else arr.reshape(())

    def read_tile(
        self,
        z: int,
        x: int,
        y: int,
        *,
        tilesize: int = 256,
        resampling: str = "bilinear",
        band: int | None = None,
    ) -> np.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. Pixel values only — the tile's
            georeferencing is defined by its `(z, x, y)`, not attached to the
            array; edge tiles are NoData-padded (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
    ) -> 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: EPSG code of `bbox`.

        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
        if self._ds.epsg == bbox_crs:
            return min_x, min_y, max_x, max_y
        transformer = Transformer.from_crs(bbox_crs, self._ds.epsg, always_xy=True)
        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]
        return min(xs), min(ys), max(xs), max(ys)

    def _world_to_pixel(self, x: float, y: float, point_crs: int) -> 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: EPSG code of `(x, y)`.

        Returns:
            `(col, row)` integer pixel indices (floored).
        """
        if self._ds.epsg != point_crs:
            transformer = Transformer.from_crs(point_crs, self._ds.epsg, always_xy=True)
            x, y = transformer.transform(x, y)
        inv = gdal.InvGeoTransform(self._ds._raster.GetGeoTransform())
        col, row = gdal.ApplyGeoTransform(inv, x, y)
        return int(math.floor(col)), int(math.floor(row))

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

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

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

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

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

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

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

is_cog property

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

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

Examples:

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

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

• Plain GeoTIFFs and MEM datasets return False:

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

• Use in a conditional pipeline:

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

_is_cog_cheap(path) staticmethod

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

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

Parameters:

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

Returns:

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

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

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

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

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

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

validate_cog(strict=False, config=None)

Validate the backing file as a COG.

Parameters:

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

Returns:

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

Raises:

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

Examples:

• Validate and branch on the result:

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

• Strict mode promotes warnings to errors:

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

• Inspect structural details from the report:

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

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

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

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

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

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

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

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

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

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

info(config=None)

Return structured COG metadata for the backing file.

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

Parameters:

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

Returns:

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

Raises:

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

Examples:

• Inspect a COG's compression and overview pyramid:

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

• Read the tile size and band count:

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

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

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

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

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

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

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

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

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

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

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

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.

Parameters:

Name	Type	Description	Default
bbox	tuple[float, float, float, float]	(min_x, min_y, max_x, max_y) window in bbox_crs.	required
dst_width	int | None	Output width in pixels. Defaults to the source window width (no decimation).	None
dst_height	int | None	Output height in pixels. Defaults to the source window height.	None
bbox_crs	int	EPSG code of bbox. Reprojected to the dataset CRS when different. Defaults to 4326 (WGS84 lon/lat).	4326
resampling	str	One of nearest, bilinear, cubic, cubicspline, lanczos, average, mode.	'bilinear'
band	int | None	0-based band index. None reads all bands.	None

Returns:

Type	Description
ndarray	numpy.ndarray: (rows, cols) for a single band, or
ndarray	(bands, rows, cols) for all bands; always sized
ndarray	dst_height x dst_width (the requested output size). Pixel values
ndarray	only — no transform, bounds, or CRS is attached.

Raises:

Type	Description
ValueError	Unknown resampling.
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.

Examples:

• Read a 256x256 decimated thumbnail of a bbox:

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

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

def read_part(
    self,
    bbox: tuple[float, float, float, float],
    *,
    dst_width: int | None = None,
    dst_height: int | None = None,
    bbox_crs: int = 4326,
    resampling: str = "bilinear",
    band: int | None = None,
) -> np.ndarray:
    """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.

    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: EPSG code of `bbox`. Reprojected to the dataset CRS
            when different. Defaults to 4326 (WGS84 lon/lat).
        resampling: One of `nearest`, `bilinear`, `cubic`,
            `cubicspline`, `lanczos`, `average`, `mode`.
        band: 0-based band index. `None` reads all bands.

    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.

    Raises:
        ValueError: Unknown `resampling`.
        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`.

    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)

            ```
    """
    if resampling not in _RESAMPLING_ALG:
        raise ValueError(
            f"unknown resampling {resampling!r}; "
            f"choose from {sorted(_RESAMPLING_ALG)}"
        )
    ds = self._ds._raster
    min_x, min_y, max_x, max_y = self._reproject_bbox(bbox, bbox_crs)
    inv = gdal.InvGeoTransform(ds.GetGeoTransform())
    px_tl, py_tl = gdal.ApplyGeoTransform(inv, min_x, max_y)
    px_br, py_br = gdal.ApplyGeoTransform(inv, 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
    alg = _RESAMPLING_ALG[resampling]
    source = ds if band is None else ds.GetRasterBand(band + 1)

    fully_inside = (
        ix0 == req_xoff
        and iy0 == req_yoff
        and ix1 == req_xoff + req_xsize
        and iy1 == req_yoff + req_ysize
    )
    if fully_inside:
        return np.asarray(
            source.ReadAsArray(
                ix0,
                iy0,
                ix1 - ix0,
                iy1 - iy0,
                buf_xsize=out_w,
                buf_ysize=out_h,
                resample_alg=alg,
            )
        )

    # 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
    return out

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

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

Pulls from a coarse overview when one exists, so previewing a huge COG is cheap.

Parameters:

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

Returns:

Type	Description
ndarray	numpy.ndarray: The downsampled array, (rows, cols) or
ndarray	(bands, rows, cols). Pixel values only — no transform, bounds,
ndarray	or CRS is attached to the returned array.

Raises:

Type	Description
ValueError	Unknown resampling.

Examples:

• Build a 128px thumbnail of a single band:

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

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

def preview(
    self,
    *,
    max_size: int = 1024,
    resampling: str = "bilinear",
    band: int | None = None,
) -> np.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.

    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)`. Pixel values only — no transform, bounds,
        or CRS is attached to the returned array.

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

            ```
    """
    if resampling not in _RESAMPLING_ALG:
        raise ValueError(
            f"unknown resampling {resampling!r}; "
            f"choose from {sorted(_RESAMPLING_ALG)}"
        )
    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))
    alg = _RESAMPLING_ALG[resampling]
    ds = self._ds._raster
    source = ds if band is None else ds.GetRasterBand(band + 1)
    return np.asarray(
        source.ReadAsArray(buf_xsize=out_w, buf_ysize=out_h, resample_alg=alg)
    )

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

Sample band value(s) at a single coordinate.

Parameters:

Name	Type	Description	Default
x	float	X / longitude / easting in point_crs.	required
y	float	Y / latitude / northing in point_crs.	required
point_crs	int	EPSG code of (x, y). Reprojected to the dataset CRS when different. Defaults to 4326.	4326
band	int | None	0-based band index. None samples all bands.	None

Returns:

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

Raises:

Type	Description
OutOfBoundsError	The point falls outside the raster extent.

Examples:

• Sample all bands at a lon/lat coordinate:

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

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

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

    Args:
        x: X / longitude / easting in `point_crs`.
        y: Y / latitude / northing in `point_crs`.
        point_crs: EPSG code of `(x, y)`. Reprojected to the dataset CRS
            when different. Defaults to 4326.
        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`. Pixel values only — no
        coordinate metadata is attached.

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

            ```
    """
    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)
    arr = np.asarray(source.ReadAsArray(col, row, 1, 1))
    return arr.reshape(-1) if band is None else arr.reshape(())

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

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

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

Parameters:

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

Returns:

Type	Description
ndarray	numpy.ndarray: A (tilesize, tilesize) or
ndarray	(bands, tilesize, tilesize) array. Pixel values only — the tile's
ndarray	georeferencing is defined by its (z, x, y), not attached to the
ndarray	array; edge tiles are NoData-padded (see :meth:read_part).

Raises:

Type	Description
OutOfBoundsError	The tile does not intersect the raster.

Examples:

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

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

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

def read_tile(
    self,
    z: int,
    x: int,
    y: int,
    *,
    tilesize: int = 256,
    resampling: str = "bilinear",
    band: int | None = None,
) -> np.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. Pixel values only — the tile's
        georeferencing is defined by its `(z, x, y)`, not attached to the
        array; edge tiles are NoData-padded (see :meth:`read_part`).

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

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

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

Validation

pyramids.dataset.cog.validate

COG validation wrapping osgeo_utils sample validator.

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

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

ValidationReport dataclass

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

Attributes:

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

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

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

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

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

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

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

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

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

                ```
        """
        return self.is_valid

_osgeo_validate(path)

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

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

Parameters:

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

Returns:

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

Raises:

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

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

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

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

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

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

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

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

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

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

_fallback_validate(path)

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

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

Heuristic limitations

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

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

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

Parameters:

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

Returns:

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

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

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

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

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

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

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

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

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

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

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

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

Parameters:

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

Returns:

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

Raises:

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

Examples:

• Validate a local COG and inspect the report:

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

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

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

• Validate a cloud-hosted COG via VSI path:

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

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

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

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

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

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

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

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

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

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

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

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

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