Metadata-Version: 2.5
Name: raster-mosaic
Version: 0.1.0
Summary: Generic reproject + mosaic tool for GDAL-readable rasters, respecting each input's own internal valid-data mask band, run on PBS Pro or SLURM
Project-URL: Repository, https://gitlab.com/jrsrp/sys/raster_mosaic
Project-URL: Documentation, https://gitlab.com/jrsrp/sys/raster_mosaic/-/blob/main/docs/running_on_hpc.md
Project-URL: Issues, https://gitlab.com/jrsrp/sys/raster_mosaic/-/issues
Project-URL: Changelog, https://gitlab.com/jrsrp/sys/raster_mosaic/-/blob/main/CHANGELOG.md
Author-email: rdenham <robert.denham@uq.edu.au>
License-Expression: MIT
License-File: LICENSE
Keywords: cog,gdal,hpc,mosaic,pbs,raster,reproject,slurm
Classifier: Development Status :: 4 - Beta
Classifier: Environment :: Console
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: GIS
Requires-Python: >=3.12
Requires-Dist: click>=8.3.1
Requires-Dist: numpy>=2.0
Requires-Dist: rasterio>=1.4.3
Provides-Extra: flow
Requires-Dist: dask-jobqueue>=0.9.0; extra == 'flow'
Requires-Dist: prefect-dask>=0.3.5; extra == 'flow'
Requires-Dist: prefect>=3.4.1; extra == 'flow'
Requires-Dist: tomlkit>=0.13; extra == 'flow'
Description-Content-Type: text/markdown

# raster_mosaic

Reprojects and mosaics an arbitrary set of GDAL-readable rasters onto one common CRS/pixel grid,
respecting each input's own validity -- an internal valid-data mask band (`gdalinfo`'s "Mask Flags:
PER_DATASET") or a NoData value, whichever the inputs actually carry -- and choosing which of the
two the output uses accordingly (see "Output validity: mask band vs NoData value" below), rather
than forcing every run through a single reserved NoData value regardless of what the inputs have.
Runs locally for testing, or as a PBS Pro / SLURM job array + final job on an HPC cluster.

To run it on a cluster -- with the bundled Prefect + Dask flow, or with plain PBS/SLURM job
scripts -- see [docs/running_on_hpc.md](docs/running_on_hpc.md).

Features:

- **Respects each input's validity** -- internal valid-data mask bands and NoData values are both
  carried through reprojection and mosaicking, and the output uses one or the other (default
  `"auto"`: decided from the inputs; set `output_validity = "nodata"` to always get a NoData
  value). See "Output validity: mask band vs NoData value".
- **Target-aligned pixels** -- every input is warped with `-tap`, so the output origin is an exact
  multiple of the resolution in the target CRS (EPSG:3577 by default).
- **Policy-conforming COG** -- 512×512 blocks, ZSTD level 9, product-type predictor and overview
  resampling, `PROCESSING_*` provenance tags. See "COG policy compliance".
- **Optional statistics and histograms** embedded in the file (off by default: `compute_stats`,
  `compute_histogram`). See "Band statistics" and "Histograms".
- **Parallel on PBS Pro / SLURM** -- reprojection runs as a job array or a Dask worker fan-out,
  with up to `ncpus` tiles warped concurrently in each job.
- **Containerised** -- cluster jobs run inside the published apptainer image, so the host needs no
  GDAL install. See "HPC image".

Works with any band count and dtype. A single-band uint8 classification/palette raster is a common
case (resampling defaults to nearest-neighbour for exactly this reason -- see
`MosaicConfig.resampling`), but nothing here is specific to that.

## Install

From PyPI -- no clone needed:

```bash
uvx raster-mosaic --help                             # run it without installing
uvx --from 'raster-mosaic[flow]' raster-mosaic run   # the Prefect + Dask flow
uv add raster-mosaic                                 # or add it to your own project
```

The local `reproject-tile` / `mosaic` commands need the GDAL command-line tools (`gdalwarp`,
`gdalbuildvrt`, `gdal_translate`) on `PATH`. On a cluster, the jobs run inside the container
instead (see "HPC image" below), so the host only submits.

For development:

```bash
uv sync --group dev --group test
```

## Commands

```bash
# 1. Discover inputs (cfg.input_globs in mosaic.toml) -> JSONL manifest
uv run raster-mosaic build-manifest --config mosaic.toml --manifest-out run/manifest.jsonl

# 1'. ...or from an explicit list of paths (one per line, '-' for stdin) instead of globs
uv run raster-mosaic build-manifest --config mosaic.toml --manifest-out run/manifest.jsonl \
    --inputs-from inputs.txt --summary-out run/summary.json

# 2. Reproject each input onto the common grid (local -- one line of the manifest at a time)
uv run raster-mosaic reproject-tile --manifest run/manifest.jsonl --index 1 --config mosaic.toml

# 2'. ...or submit the whole manifest as a PBS/SLURM job array
uv run raster-mosaic submit-reproject --scheduler pbs --manifest run/manifest.jsonl \
    --config mosaic.toml --container raster_mosaic.sif --account your_account

# 3. Mosaic all reprojected tiles into the final COG
uv run raster-mosaic mosaic --config mosaic.toml

# 3'. ...or submit as a single job, after the reproject array job finishes
uv run raster-mosaic submit-mosaic --scheduler pbs --config mosaic.toml \
    --container raster_mosaic.sif --account your_account \
    --depends-on <reproject_array_jobid>
```

`--scheduler` is `pbs` or `slurm` on both `submit-*` commands.

## Chunking many small/fast inputs

If your manifest has many inputs that each reproject quickly, submitting one PBS/SLURM array task
per input can be counter-productive: per-job scheduling overhead (apptainer startup, PBS/SLURM
dispatch) can dwarf the actual work. A real job was observed to spend most of its ~11s walltime on
overhead alone, and `cpupercent=100` in the PBS accounting despite requesting 4 CPUs -- a single
small warp doesn't scale with more threads, so those extra CPUs went unused.

`--chunk-size` on `submit-reproject` fixes both: each array task processes a *chunk* of several
manifest lines instead of just one (shrinking the array and amortizing overhead across more work
per job), and the tiles within a chunk run *concurrently* -- up to `--ncpus` at a time -- instead
of handing all those CPUs to one warp that can't use them:

```bash
uv run raster-mosaic submit-reproject --scheduler pbs --manifest run/manifest.jsonl \
    --config mosaic.toml --container raster_mosaic.sif --account your_account \
    --chunk-size 20 --ncpus 4
```

Defaults to 1 (one task per input, unchanged from before) since the right chunk size depends on
how many inputs you have and how fast each one is -- there's no universally-safe non-1 default.

## Driving it from a workflow orchestrator

The `submit-*` commands are one way to run the three stages on a cluster. A workflow orchestrator
(e.g. a Prefect flow fanning tasks out over Dask workers -- this package ships one,
`raster_mosaic.flow`, see [docs/running_on_hpc.md](docs/running_on_hpc.md)) can run them itself
instead, calling the same CLI inside this package's container: one `build-manifest`, then one `reproject-tile --index i`
per chunk, then one `mosaic`. A few options exist mainly for that:

- `build-manifest --inputs-from FILE` takes an explicit input list instead of the config's
  `input_globs` (exactly one of the two must be given) -- e.g. the outputs of an upstream workflow
  step. One path per line; blank lines and `#` comments are ignored; `-` reads stdin; relative
  paths are resolved against the current directory and stored absolute. Every listed path must
  exist.
- `build-manifest --summary-out FILE` writes a JSON summary next to the human-readable output:

  ```json
  {"manifest": "/abs/run/manifest.jsonl", "n_inputs": 3, "chunk_size": 2, "n_chunks": 2,
   "validity": {"mode": "mask", "nodata": null}}
  ```

  `n_chunks` is how many `reproject-tile --index` values to run (1..`n_chunks`). It assumes the
  config's own `chunk_size`, so don't also pass a different `--chunk-size` to `reproject-tile`.
- `reproject-tile --skip-existing` (or `skip_existing_tiles = true` in the config) skips inputs
  whose reprojected tile already exists, so a retried task only redoes what didn't finish. Tiles
  only get their final `.tif` name once every step on them is done, so an existing tile is always
  complete -- but it may have been built with a *different* config, which is why this is off by
  default.

`build-manifest` (both input modes) also rejects inputs that share a file stem, since reprojected
tiles are named by stem alone and same-stemmed inputs would otherwise overwrite each other's tile.

## HPC image

`submit-reproject` / `submit-mosaic` jobs run entirely inside the `--container` image via
`apptainer exec` -- no host-side `uv sync`, no live project checkout at job time. `.gitlab-ci.yml`
builds [`Dockerfile`](Dockerfile) (this package layered on the official GDAL CLI image) and pushes
it to this project's GitLab Container Registry on every push -- `:latest` tracks the default branch
and any git tag; every branch build is also pushed tagged with its ref slug and short commit SHA,
so a specific build stays pullable after the branch moves on.

Pull it with `apptainer` -- no docker needed on the HPC side, no local image build at all:

```bash
apptainer pull raster_mosaic.sif docker://registry.gitlab.com/jrsrp/sys/raster_mosaic:<tag>
```

The image is public, so no login is needed. Pick `<tag>` to match your installed package
(`raster-mosaic --version`): `X.Y.Z` for a release, the short commit SHA for a branch build, or
`latest` for the newest default-branch build. The Prefect flow refuses a mismatched image -- see
[docs/running_on_hpc.md](docs/running_on_hpc.md#the-container). Then point `--container` at the
pulled `raster_mosaic.sif`.

Pulling once and reusing the `.sif` for a whole job array (as in the examples above) is the
simplest. `apptainer exec docker://...` also works directly: the converted image is cached
(`~/.apptainer/cache`), so with that cache on a shared filesystem and warmed once beforehand, the
array tasks reuse it -- but each still checks the tag with the registry (compute nodes need network
access), and a mutable tag like `latest` could change mid-run, so use a fixed one.

## Running `submit-reproject` / `submit-mosaic` itself from inside a container

The image above has no PBS/SLURM client (`qsub`/`sbatch`) -- it only needs to run the actual
`reproject-tile`/`mosaic` work, which the job templates already invoke via `apptainer exec`. But if
you're also running `raster-mosaic submit-reproject`/`submit-mosaic` itself from inside a
container (rather than directly on a login node where `qsub`/`sbatch` are on `PATH`), pass
`--submit-via`:

```bash
uv run raster-mosaic submit-reproject --scheduler pbs ... --submit-via "ssh localhost"
```

This prefixes the `qsub`/`sbatch` invocation with the given (shell-split) command -- `ssh
localhost qsub ...` reaches the real PBS/SLURM client on the underlying host from inside a
container that doesn't have one installed. Every path this command itself hands to the rendered
job script (`WORKDIR`, `MANIFEST`, `CONFIG`, the job script's own temp-file path, and the job's own
`-o`/`--output` log location) is already absolute for exactly this reason: an SSH-wrapped
submission starts from the SSH session's home directory, not this command's own working directory,
so a relative path would resolve against the wrong place.

## Site-specific PBS/SLURM job configuration

The rendered job scripts (`hpc/templates/*.sh.tmpl`) are deliberately generic -- no site's
particular account, filesystem layout, or scheduler quirks are baked in. Where a site needs
something beyond the common flags (`--account`, `--walltime`, `--mem`, `--ncpus`, `--partition`,
`--qos`, `--module-load`), four escape hatches on `submit-reproject`/`submit-mosaic` cover it, all
opt-in and no-op by default:

- **`--extra-directive`** -- one extra raw scheduler directive line (`#PBS ...` / `#SBATCH ...` as
  appropriate for `--scheduler`), inserted verbatim into the job script. Repeatable for more than
  one line. For anything with no dedicated flag, e.g. a site that requires
  `#PBS -W sandbox=PRIVATE` before a job can use node-local scratch.
- **`--select-extra`** -- PBS only. Extra `:name=value` chunk resource(s) (leading `:` required)
  appended onto the auto-built `select=1:ncpus=...:mem=...` string, e.g.
  `:flashreq=500g:flashtmp=500g` on a site that provisions flash/NVMe scratch as its own `select`
  chunk resource rather than a top-level `-l` flag. Ignored on `--scheduler slurm` -- SLURM
  expresses extra resources via its own `--gres`/`--constraint` directives, which belong in
  `--extra-directive` instead.
- **`--apptainer-bind`** -- an extra `src[:dst[:opts]]` bind spec for the `apptainer exec`
  invocation itself (not a scheduler directive -- a flag to `apptainer`). See "Container visibility
  of site-specific paths" below; almost always needed alongside `--stage-dir` on a path apptainer
  doesn't auto-bind. Repeatable -- every occurrence is comma-joined into one `--bind` argument.
- **`--stage-dir`** (mosaic only) -- see "Staging the mosaic build through fast scratch" below.

Example: a PBS Pro site where `#PBS -W sandbox=PRIVATE` alone makes PBS create a per-job private
directory under `/flash/tmp` and cd the job into it, destroying it automatically when the job ends
(no `--select-extra` needed for capacity in this case -- only add it if your site *also* requires
an explicit resource reservation to get that scratch allocated, e.g.
`--select-extra ":flashreq=500g:flashtmp=500g"`):

```bash
uv run raster-mosaic submit-mosaic --scheduler pbs --config mosaic.toml \
    --container raster_mosaic.sif --account your_account \
    --extra-directive "#PBS -W sandbox=PRIVATE" \
    --apptainer-bind "/flash:/flash" \
    --stage-dir '$PBS_JOBDIR' \
    --depends-on <reproject_array_jobid>
```

`sandbox=PRIVATE`'s directory is exposed via an environment variable that PBS Pro populates once
the job starts -- `$PBS_JOBDIR` or `$TMPDIR` depending on your PBS version/site configuration, not
a fixed name across every site. Check yours with a trivial interactive job (`qsub -I -W
sandbox=PRIVATE` then `env | grep -i job`) before relying on either.

### Container visibility of site-specific paths

`apptainer exec` only auto-binds a small default set of paths inside the container -- your home
directory, `/tmp`, `/proc`, `/sys`, `/dev`, cwd, plus whatever a site added globally to
`apptainer.conf`'s `bind path`. A path outside that set -- most commonly wherever `--stage-dir`
points, e.g. a site's flash/NVMe mount -- exists on the host but is **invisible inside the
container** without an explicit bind. The symptom is a `FileNotFoundError` for a path that
demonstrably exists on the host (confirm with `ls` outside the container), which looks like a
`--stage-dir`/env-var resolution bug but usually isn't one -- it's a missing `--apptainer-bind`.
Quick way to confirm before chasing anything else: `apptainer exec $CONTAINER ls /your/mount` --
if that fails while `ls /your/mount` on the bare host succeeds, add
`--apptainer-bind "/your/mount:/your/mount"` (or whatever the mount's top-level path is; a bind of
the parent covers subdirectories created under it later, e.g. per-job `sandbox=PRIVATE`
directories, since apptainer resolves the bind live at container start, not as a snapshot).

### Staging the mosaic build through fast scratch

`mosaic`/`submit-mosaic --stage-dir` (equivalently, `MosaicConfig.stage_dir` in `mosaic.toml`)
copies the reprojected tiles into a temp subdirectory of `--stage-dir`, builds the VRT and COG
there, then copies only the finished COG back out to `output_dir` -- useful when
`reprojected_dir`/`output_dir` sit on a slower shared filesystem and that's the mosaic stage's
bottleneck compared to node-local fast scratch (flash/NVMe). Omit it (the default) to build
directly in `reprojected_dir`/`output_dir`, unchanged.

`--stage-dir`'s value on `submit-mosaic` is inserted **verbatim** into the rendered job script's
`raster-mosaic mosaic --stage-dir "..."` invocation, not passed through `qsub -v`/`sbatch
--export` the way `WORKDIR`/`CONTAINER`/`CONFIG` are -- this lets it be a scheduler-populated
environment variable (e.g. `$PBS_JOBDIR`, `$TMPDIR`, SLURM's `$SLURM_TMPDIR`) that's only resolved
once the job is actually running, in addition to a literal path. Quote it in single quotes on the
command line as in the example above so your local shell doesn't expand it early.

Provisioning the scratch directory itself (if your site requires an explicit resource request for
it, as PBS flash/NVMe chunk resources typically do) is exactly what `--select-extra`/
`--extra-directive` are for; making it visible inside the container is `--apptainer-bind`'s job
(see "Container visibility of site-specific paths" above) -- `--stage-dir` only tells `run_mosaic`
where to build once the scratch exists and is visible.

If your site's `sandbox=PRIVATE` (or equivalent) also changes the job's initial working directory,
that's independent of `--stage-dir` and nothing to worry about: the rendered job script does its
own `cd "$WORKDIR"` right after the directives (so logs/manifest/config still resolve against your
shared-filesystem run directory), and `--stage-dir`'s value is passed as its own absolute path/env
var reference, not read relative to whatever the cwd happens to be at that point.

## Config (`mosaic.toml`)

Unknown keys are an error, not ignored -- they're either a typo or an option the running version
doesn't have (e.g. a container image older than the install that wrote the config).

```toml
input_globs = ["egimages/*.tif"]   # glob pattern(s); ** for recursive matching. Omit if
                                    # using build-manifest --inputs-from instead
reprojected_dir = "run/reprojected"
output_dir = "run/out"
output_name = "mosaic.tif"

resolution = 10.0

dst_crs = "EPSG:3577"              # default; GDA94 Australian Albers
resampling = "near"                # "near" is the only safe choice for categorical/palette data
product_type = "categorical"       # default; "continuous" or "mask" also accepted
ncpus = 4

output_validity = "auto"           # default; "mask" or "nodata" also accepted -- see below
# output_nodata = 255.0            # only meaningful (and only ever required) with output_validity
                                    # = "nodata" -- see "Output validity" below

# optional -- omit entirely to skip land masking
land_mask = "land_mask/Land_any_epsg3577_100m.tif"

# optional -- omit entirely to build directly in reprojected_dir/output_dir; see "Staging the
# mosaic build through fast scratch" above. Usually site-specific, so more often passed per-run
# via --stage-dir instead of committed here.
# stage_dir = "/flash/tmp"

# optional -- extra -mo KEY=VALUE tags on the final mosaic, on top of cog_policy.md's own
# mandated PROCESSING_* tags. See "Extra metadata (extra_mo_tags)" below.
# [extra_mo_tags]
# PIXEL_SCALING = "GAIN:0.01 OFFSET:0.0 UNITS:fraction"

# optional -- GDAL's own native per-band scale/offset (-a_scale/-a_offset), applied uniformly to
# every band. Distinct from extra_mo_tags -- see "Extra metadata" below.
# output_scale = 0.01
# output_offset = 0.0

# optional -- extra metadata scoped to one band, keyed by band number as a string. See "Extra
# metadata" below.
# [extra_band_mo_tags.1]
# PIXEL_SCALING = "GAIN:0.01 OFFSET:0.0 UNITS:fraction"

# optional -- per-band scale/offset, for the (unusual, not recommended) case where bands genuinely
# need to differ. Mutually exclusive with output_scale/output_offset -- see "Extra metadata" below.
# [extra_band_scale]
# 1 = 0.01
# [extra_band_offset]
# 1 = 0.0

# optional -- compute and embed GDAL's standard per-band STATISTICS_* metadata. See "Band
# statistics" below.
# compute_stats = "exact"          # or "approx"; default "none" adds nothing

# optional -- embed a per-band histogram (STATISTICS_HISTO* metadata). See "Histograms" below.
# compute_histogram = true         # default false

# optional -- skip inputs whose reprojected tile already exists (for retries). See "Driving it
# from a workflow orchestrator" above.
# skip_existing_tiles = true       # default false

# optional -- build overviews on each reprojected tile so run_mosaic can reuse them instead of
# recomputing the mosaic's own overviews from scratch. See "Reusing tile overviews" below.
# build_tile_overviews = true      # default false
```

`product_type` drives the final mosaic's compression/predictor/overview-resampling and
`PROCESSING_*` provenance tags to match [`cog_policy.md`](cog_policy.md) -- see "COG policy
compliance" below. It's unrelated to `resampling`, which only affects the reprojection step.

`land_mask` is optional -- when set, every non-land pixel (per the mask raster, pixel value 1 =
land) is excluded from each tile's valid-data mask during `reproject-tile`, combined with (never
overriding back to valid) whatever the tile's own mask already marked invalid. Can be in any
CRS/resolution GDAL can read; no manual reprojection/prep step needed. Applied per-tile in the
parallel reproject stage, not once at the final serial mosaic step, so the cost is spread across
the already-parallel array job -- see `gdal_ops._apply_land_mask`.

## Extra metadata (`extra_mo_tags`)

`extra_mo_tags` is a plain escape hatch for embedding whatever extra `-mo KEY=VALUE` tags a caller
wants on the final mosaic, on top of (never instead of) cog_policy.md's own mandated `PROCESSING_*`
tags -- keys can't collide with those. This tool makes no attempt to compute the values itself:
`gdalbuildvrt` doesn't carry per-tile dataset metadata onto the VRT `run_mosaic` builds and reads
from (confirmed empirically -- a source tile's own tags simply aren't there once it's part of a
VRT), so there's nothing generic to summarise from the inputs even if this tried to. Working out
what should go here (e.g. an organisation-specific summary of the inputs' own provenance/scaling
metadata) is deliberately left to a site/product-specific wrapper around this package, same as
`land_mask`/`stage_dir`'s site-specific values.

`output_scale`/`output_offset` are the same passthrough idea but for GDAL's own native per-band
scale/offset (`gdalinfo`'s "Offset"/"Scale" fields, set via `gdal_translate -a_scale`/`-a_offset`)
rather than an arbitrary metadata string -- distinct from, and can be set alongside, an
`extra_mo_tags` entry describing the same gain/offset relationship in a custom format. Applied
uniformly to every band (confirmed empirically -- `-a_scale`/`-a_offset` take one value each, not
one per band). `None` (the default, for either) omits the flag entirely rather than writing GDAL's
own default (scale 1, offset 0) explicitly.

Both `-mo`/`-dmo` (`extra_mo_tags`) and `-a_scale`/`-a_offset` are **dataset-level only** -- neither
`gdal_translate` flag can target an individual band. (`output_scale`/`output_offset` are fine
regardless, since GDAL's Scale/Offset has no dataset-level form at all -- it's purely a
`GDALRasterBand` property, so even one uniform `-a_scale`/`-a_offset` value already lands on every
band correctly, which is what QGIS and other GDAL-based tools read for display.)

For genuine band-scoped arbitrary metadata, use `extra_band_mo_tags` instead -- keyed by band
number as a string (`"1"`, not `1`; TOML table keys are always strings), each value its own
`KEY=VALUE` mapping. This works differently from `extra_mo_tags`: since `gdal_translate` has no
band-scoped metadata flag at all, `run_mosaic` instead injects `<Metadata>`/`<MDI>` XML elements
directly into its own `_mosaic.vrt` (`gdal_ops._write_band_metadata`) before the final
`gdal_translate -of COG` reads it -- confirmed empirically that a VRT's per-band `<Metadata>`
survives onto the real output band, even though `gdalbuildvrt` never writes any there itself.
Editing that VRT isn't the "invalidates COG structure" case cog_policy.md's `gdal_edit.py` note
warns about -- it's `run_mosaic`'s own working intermediate, not the released COG product. A band
number with no matching band in the built mosaic raises rather than silently dropping the tag.

**Prefer GDAL's own Scale/Offset (`output_scale`/`output_offset`) over an `extra_mo_tags` string
like `PIXEL_SCALING` whenever you can.** The official mechanism is what QGIS and other GDAL-based
tools actually read for display -- a custom string tag isn't understood by anything unless it
parses that specific key itself. Use a custom tag (`PIXEL_SCALING` or similar) alongside it only
for backward compatibility with existing tooling/consumers that already expect that convention;
don't use it as a substitute for setting the real Scale/Offset.

**We recommend keeping scale/offset consistent across every band of a mosaic** -- i.e. reach for
`output_scale`/`output_offset`, not `extra_band_scale`/`extra_band_offset`, in the common case.
Different scale/offset per band *is* possible via `extra_band_scale`/`extra_band_offset` (same
VRT-injection mechanism as `extra_band_mo_tags`, via `gdal_ops._write_band_scale_offset`), but
isn't recommended -- it's easy for downstream consumers to silently misread pixel values without
per-band-aware code. The two mechanisms are mutually exclusive by construction: `MosaicConfig`
rejects setting `output_scale`/`output_offset` alongside `extra_band_scale`/`extra_band_offset`,
since `-a_scale`/`-a_offset` on the final `gdal_translate` step unconditionally overwrites whatever
per-band `<Scale>`/`<Offset>` the VRT already has (confirmed empirically) -- silently discarding
the per-band values rather than erroring would be worse than refusing to run.

## Band statistics

`compute_stats` ("none" default / "exact" / "approx") computes and embeds GDAL's standard per-band
`STATISTICS_MINIMUM`/`MAXIMUM`/`MEAN`/`STDDEV`/`VALID_PERCENT` metadata -- e.g. for consumers
expecting them the way HFA-produced imagery traditionally carries them. Unlike `extra_mo_tags`/
`extra_band_mo_tags`/`output_scale`/`extra_band_scale` above, this isn't a passthrough: GDAL
computes the values itself, directly from the mosaic's own output pixels, as part of the same final
`gdal_translate` call (`-stats`/`-approx_stats`) -- there's nothing for a caller to supply. Confirmed
empirically that this embeds the statistics in the output file itself at creation time (no
`.aux.xml` sidecar) and excludes invalid pixels from them.

**Mask-band mosaics compute their own statistics.** When the run resolves to a mask band (see
"Output validity" below), GDAL's statistics can't be trusted: on a VRT of masked tiles that don't
overlap each other (side-by-side tiles, or tiles with gaps between them), GDAL takes a shortcut
that ignores the masked pixels, reporting `STATISTICS_VALID_PERCENT=100` and a standard deviation
divided by the total rather than the valid pixel count (confirmed with GDAL 3.11; min/max/mean
are right, and there's no GDAL option to turn the shortcut off). So in mask mode, `run_mosaic`
computes the statistics itself, in one parallel pass over the mosaic, and doesn't use `-stats`/
`-approx_stats` at all. These are always exact, even with `compute_stats = "approx"`, so `"approx"`
gives no speed-up in mask mode. NoData-mode runs still use GDAL's own flags. Mosaics made before
this fix with `compute_stats` in mask mode may carry wrong `STATISTICS_STDDEV`/`VALID_PERCENT`
values.

`"exact"` reads every pixel (accurate, but a second full pass over a potentially very large mosaic);
`"approx"` samples via overviews instead (much cheaper, approximate values) -- pick based on how
much that cost matters for your mosaic size. `cog_policy.md` is silent on statistics entirely, so
this is purely a convenience, not a compliance requirement.

Note that `"approx"`'s speed depends on overviews actually being available to sample from -- with
no overviews anywhere in `run_mosaic`'s input VRT (the default, unless `build_tile_overviews` is
also set -- see below), `"approx"` still falls back to a slow, sparse full-resolution read before
it can produce anything, which can dominate a large mosaic's runtime. `"exact"` always reads every
pixel regardless of overviews, so `build_tile_overviews` doesn't speed it up at all -- if you need
`compute_stats`'s numbers to match a true full-resolution min/max exactly (`"approx"`'s sampling can
miss rare/sparse extreme values entirely, by design), use `"exact"`, not `"approx"` sped up by
overviews.

## Histograms

`compute_histogram = true` (default `false`) embeds a per-band histogram in the final mosaic, as
the `STATISTICS_HISTOBINVALUES` / `STATISTICS_HISTOMIN` / `STATISTICS_HISTOMAX` /
`STATISTICS_HISTONUMBINS` / `STATISTICS_HISTOBINFUNCTION` band metadata items that HFA-produced
imagery and RIOS's `calcstats` use. GDAL's own histogram support isn't used because it only
persists histograms to a `.aux.xml` sidecar, never into the GeoTIFF itself.

- It's computed over the mosaic (its VRT), not by adding up per-tile histograms: where tiles
  overlap, only the pixel the mosaic actually keeps is counted. Invalid pixels (mask band or NoData
  value) and non-finite floats are excluded.
- It's exact: every valid pixel is read twice, once to find each band's minimum/maximum (bin edges
  have to be fixed before counting starts) and once to count. Both passes run in parallel over
  windows, up to `ncpus` at a time.
- Binning: an integer band whose valid values span at most 256 distinct values gets one bin per
  value (`direct`). Any other band, including every float band, gets 256 equal-width bins between
  its minimum and maximum (`linear`). A band with no valid pixels gets no histogram.
- The tags are injected into the VRT the same way as `extra_band_mo_tags`, so they're written when
  the COG is created rather than edited in afterwards. `extra_band_mo_tags` can't also set
  `STATISTICS_HISTO*` keys while this is on.

**Interaction with `compute_stats`:** GDAL's `-stats`/`-approx_stats`, combined with the `-mo` tags
`run_mosaic` always writes, drop every `STATISTICS_*` item the VRT carries, including the
histogram (confirmed empirically). So with both options on, `run_mosaic` doesn't use those flags:
the histogram pass computes the `STATISTICS_*` values itself, from the same reads, at almost no
extra cost. That makes them exact even for `compute_stats = "approx"`.

## Reusing tile overviews (`build_tile_overviews`)

`run_mosaic`'s final `gdal_translate -of COG` step always used to pass `-co
OVERVIEWS=IGNORE_EXISTING`, forcing a full-resolution overview rebuild on every run. Setting
`build_tile_overviews = true` (default `false`) has `reproject-tile` build internal overviews on
each reprojected tile as it's produced -- using the same `OVERVIEW_RESAMPLING`
`cog_policy.md`/`product_type` mandate for the final mosaic -- so `run_mosaic` can potentially reuse
them instead (`-co OVERVIEWS=FORCE_USE_EXISTING`), skipping the recompute. This also makes
`compute_stats = "approx"` (above) much cheaper, since GDAL's approximate-stats sampling prefers
overviews over full-resolution reads when they're available.

`run_mosaic` never trusts this blindly: it checks that every reprojected tile actually has
overviews *and* is tagged with the exact resampling this run's `product_type` expects, falling back
to today's `IGNORE_EXISTING` behaviour (and logging why) if any tile disagrees or lacks
overviews -- e.g. `reprojected_dir` holding a mix of tiles from before and after this setting was
turned on, or from a run with a different `product_type`. This is purely a performance decision with
an always-safe fallback, so a mismatch never fails the run the way an `output_validity` conflict
does.

**The overviews are built during `reproject-tile`, not `mosaic`.** Turning this on and re-running
`mosaic` against tiles from an earlier `reproject-tile` run won't do anything -- those tiles have no
overviews on disk yet, so `run_mosaic` will fall back to `IGNORE_EXISTING` as before. You need
`build_tile_overviews = true` set before (re-)running `reproject-tile`/`submit-reproject` to
actually get the benefit on the following `mosaic` run. Building per-tile overviews costs real time
in the reproject stage itself (though it's already parallel across the array job), which is why this
stays opt-in rather than becoming the default -- worthwhile once `run_mosaic`'s own overview
build/`compute_stats="approx"` pass is the bottleneck, not before.

## COG policy compliance

`run_mosaic`'s final output follows [`cog_policy.md`](cog_policy.md) unconditionally: COG driver,
512x512 blocks, ZSTD level 9, `PREDICTOR`/`OVERVIEW_RESAMPLING` keyed off `cfg.product_type`
(`PRODUCT_TYPES` in `config.py`), and all six `PROCESSING_*` provenance tags written at creation
time. `reproject-tile`'s own per-input output is just a working intermediate for `gdalbuildvrt`,
not the released product the policy governs, so it isn't held to the same creation options.

Whether the output ends up with a mask band or a NoData value (see below) is itself following
`cog_policy.md`, not deviating from it: Section 4.7 makes an embedded NoData value the default for
"the standard large COG products", with an internal mask band as its narrower, explicit exception
("selected products where 0 is a valid data value... not used in the standard large COG products
covered by this specification"). This tool follows exactly that: a mask band only when the inputs
actually have one.

## Output validity: mask band vs NoData value

A raster can mark invalid pixels in two ways:

- **Mask band:** a separate valid/invalid layer (`gdalinfo`'s "Mask Flags: PER_DATASET"). Every
  pixel value stays usable as real data.
- **NoData value:** one reserved value means "invalid". A real pixel with that value is a problem.

Invalid pixels in the inputs are always respected while reprojecting and mosaicking, whichever way
each input marks them. `output_validity` (default `"auto"`) only decides which of the two the
reprojected tiles and the final mosaic carry.

### What `build-manifest` sees in each input

Each input is classified (from band 1) as one of:

- **mask:** has a mask band. This includes inputs with *both* a mask band and a NoData value: the
  mask takes precedence and the NoData value is ignored, as in GDAL itself (`gdalinfo`, QGIS and
  rasterio all read such a file that way). `gdalwarp` on its own does the opposite, so
  `reproject-tile` passes it `-srcnodata None` for these inputs to keep it consistent.
- **nodata=X:** has a NoData value and no mask band.
- **neither:** every pixel is valid.

### `"auto"` (the default)

| Inputs | Output | Notes |
|---|---|---|
| All masked | Mask band | No value is reserved, so nothing can collide. |
| All NoData, same value X | NoData = X | The standard product form under `cog_policy.md` §4.7. |
| Some NoData X, the rest neither | NoData = X | ⚠️ A real pixel equal to X in a "neither" input can't stay X: `gdalwarp` nudges it (e.g. 0 becomes 1) and only logs a warning. |
| All neither | Mask band | The output still needs validity: reprojection leaves fill around each tile, and the mosaic has gaps between tiles. A mask marks those invalid without reserving a value. |
| Masked mixed with unmasked | **Error** | Not guessed, because dropping a mask risks the collision this tool exists to avoid. |
| NoData values differ (e.g. 0 and 255) | **Error** | There's no single value to adopt. |

### Overriding it

- **`"mask"`** always carries a mask band, for every input; any NoData values are dropped from the
  output, though the pixels they marked stay invalid. It never alters data, so it's the safe fix
  for either error above.
- **`"nodata"`** always embeds a NoData value. Set `output_nodata` unless every input already has
  the same NoData value and none is masked. Every input's own invalid pixels (masked or NoData,
  whatever value) become `output_nodata` in the output. You're promising `output_nodata` never
  occurs as real data: a valid pixel equal to it gets nudged to a neighbouring value, with only a
  warning.

**Example: inputs with different NoData values** (A uses 0, B uses 255). Either set
`output_validity = "mask"` (both inputs' NoData pixels become masked, nothing reserved), or
`output_validity = "nodata"` with an `output_nodata` that can't be real data -- e.g. `255` for
0-100 percentage data -- in which case A's 0s and B's 255s both come out as 255.

This is resolved once for the whole run by `build-manifest` (over every input, not just one array
task's chunk) and persisted next to the manifest for `reproject-tile` to read back -- see
`src/raster_mosaic/validity.py`'s module docstring for why it can't be re-derived independently per
array task, and for the full decision logic.

## How the mask/NoData value is preserved

See the module docstrings in `src/raster_mosaic/gdal_ops.py` and `validity.py` for the full
mechanism. Briefly, for mask mode: `gdalwarp` already respects a source's internal mask band during
resampling by default (no `-srcnodata` needed), but doesn't carry it through to its output unless
asked. `reproject_one` asks via `-dstalpha`, then folds the resulting alpha band back into a real
internal mask band with `gdal_translate -mask`. `gdalbuildvrt` then respects each input's mask when
stacking overlapping tiles (a masked-invalid pixel from one tile does not overwrite valid data from
another tile underneath it), and the final `-of COG` write preserves the mask and excludes masked
pixels from overview generation. NoData mode works the same way in spirit but more simply -- a
single `gdalwarp -dstnodata` call, `gdalbuildvrt -srcnodata`/`-vrtnodata` for the same
overlapping-tile protection, and `-a_nodata` on the final `gdal_translate` -- since there's no mask
band to fold in the first place.

## Overlapping tiles

Where two or more reprojected tiles cover the same pixel, `gdalbuildvrt` resolves it per pixel as
it composites `_mosaic.vrt`'s sources: the **last-listed tile that's actually valid there wins** --
not a blend or average. If the later tile is invalid at that pixel (masked out, or NoData), its
invalidity doesn't overwrite a valid pixel from an earlier tile underneath (see "How the mask/NoData
value is preserved" above); if both are valid, the later one's value is simply used outright, with
no feathering or distance-weighting at the seam.

"Later" means later in the list `run_mosaic` hands to `gdalbuildvrt`, which is every reprojected
tile in `reprojected_dir` sorted alphabetically by filename (i.e. by each input's own stem) --
`input_globs`' own match order isn't preserved into this. There's currently no way to prefer one
overlapping tile over another for any other reason (e.g. a more recent acquisition, a cloud-free
scene, a quality flag) -- alphabetical order is all there is.

That's a meaningful distinction for anyone wanting different behaviour here:

- **A different priority order** (e.g. newest-wins, or an explicit per-input priority) only needs
  `run_mosaic` to hand `gdalbuildvrt` its tile list in a different order -- still no custom pixel
  processing, `gdalbuildvrt`'s own compositing already does the right thing once the order is right.
- **Anything other than a hard per-pixel cutover** -- feathering across the seam, distance-weighted
  blending, a quality-score-weighted average -- is not something `gdalbuildvrt`/`gdal_translate` can
  do at all. That's genuine custom pixel processing (reading tile arrays and computing the output
  directly, e.g. via `rasterio`/numpy, or a different tool entirely), a real departure from this
  package's "GDAL CLI only" architecture (see `gdal_ops.py`'s module docstring), not a small change.
