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.

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:

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:

uv sync --group dev --group test

Commands

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

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

  {"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 (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:

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

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

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

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

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.