Metadata-Version: 2.4
Name: s2mosaic
Version: 2.0.0b6
Summary: Python library for making cloud-free Sentinel-2 mosaics
Author-email: Nick Wright <nicholas.wright@dpird.wa.gov.au>
License-Expression: MIT
Project-URL: Homepage, https://github.com/DPIRD-DMA/S2Mosaic
Project-URL: Repository, https://github.com/DPIRD-DMA/S2Mosaic
Project-URL: Issues, https://github.com/DPIRD-DMA/S2Mosaic/issues
Project-URL: Changelog, https://github.com/DPIRD-DMA/S2Mosaic/blob/main/CHANGELOG.md
Keywords: sentinel-2,satellite,mosaic,remote-sensing
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Developers
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: GIS
Classifier: Topic :: Scientific/Engineering :: Image Processing
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: geopandas>=1.0.3
Requires-Dist: multiclean>=0.4.0
Requires-Dist: numba>=0.61
Requires-Dist: numpy>=1.26
Requires-Dist: omnicloudmask>=1.7.1
Requires-Dist: opencv-python-headless>=4.10
Requires-Dist: pandas>=2.0
Requires-Dist: planetary-computer>=1.0.0
Requires-Dist: pyproj>=3.6
Requires-Dist: pystac>=1.10
Requires-Dist: pystac-client>=0.8.6
Requires-Dist: python-dateutil>=2.8
Requires-Dist: rasterio!=1.4.4,<2,>=1.3
Requires-Dist: shapely>=2.0
Requires-Dist: tqdm>=4.65
Requires-Dist: urllib3>=2.0
Dynamic: license-file

<p align="center">
  <img src="https://raw.githubusercontent.com/DPIRD-DMA/S2Mosaic/main/assets/s2mosaic-heading.svg" alt="S2Mosaic" width="600">
</p>

<div align="center">

[![PyPI version](https://img.shields.io/pypi/v/s2mosaic.svg)](https://pypi.org/project/s2mosaic/) [![Python versions](https://img.shields.io/pypi/pyversions/s2mosaic.svg)](https://pypi.org/project/s2mosaic/) [![CI](https://github.com/DPIRD-DMA/S2Mosaic/actions/workflows/ci.yml/badge.svg)](https://github.com/DPIRD-DMA/S2Mosaic/actions/workflows/ci.yml) [![Downloads](https://static.pepy.tech/badge/s2mosaic)](https://pepy.tech/project/s2mosaic) [![License](https://img.shields.io/pypi/l/s2mosaic.svg)](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/LICENSE)

</div>

S2Mosaic is a Python package for creating cloud-free mosaics from Sentinel-2 satellite imagery. It can mosaic full MGRS tiles, rectangular bounds, or polygon AOIs over a chosen time window, with configurable scene ordering, masking, aggregation, and export options.

## Features

- Mosaic by MGRS grid tile (`grid_id`), rectangular bounds (`bounds`), or a single polygon (`aoi`). Bounds and AOIs can cross MGRS tile boundaries and are reprojected onto a common UTM grid in one step.
- Flexible scene ordering: by valid data percentage, oldest, newest, or a custom sort function.
- Multiple mosaic creation methods: mean, arbitrary percentile, median, medoid (per pixel, the scene closest to the per-band median, which preserves real observed spectra), or first valid pixel.
- Support for different spectral bands, including visual (RGB) composites.
- Cloud masking with OmniCloudMask, plus an SCL option that skips inference for cheaper bulk processing.
- STAC source selection: Microsoft Planetary Computer by default, Element 84 Earth Search on AWS Open Data, or Digital Earth Australia's Sentinel-2 NBART for Australia.
- Tile-streaming pipeline keeps peak memory low even for full-MGRS percentile mosaics over many scenes. The aggregation is parallelised across ~2048-pixel tiles, so only a handful of tile-sized buffers live in RAM at a time.
- Resilient to transient COG read failures. Per-scene fetches retry with exponential backoff, and a scene that still fails is logged and skipped so one bad asset doesn't abort the whole mosaic.
- Export mosaics as GeoTIFF files or return as NumPy arrays.

## Changelog

See [CHANGELOG.md](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/CHANGELOG.md) for the full release history.

## Choosing a cloud mask

S2Mosaic uses OmniCloudMask (OCM) by default for cloud and cloud-shadow masking. OCM runs much faster when an NVIDIA GPU or MPS accelerator is available. When compute is limited or throughput matters more than mask quality, pass `cloud_mask="SCL"` to skip the deep-learning model and use the Sentinel-2 L2A Scene Classification Layer instead.

## Choosing a source

`source` selects where scenes come from. All three return spectral bands on the same `reflectance * 10000` scale, so switching source doesn't change what the numbers mean.

| `source` | Provider | Coverage | Product | Notes |
|---|---|---|---|---|
| `"MPC"` (default) | Microsoft Planetary Computer | Global | ESA L2A (Sen2Cor) | SAS-signed URLs |
| `"AWS"` | Element 84 Earth Search | Global | ESA L2A (Sen2Cor) | Public COGs, no auth |
| `"DEA"` | Digital Earth Australia (Geoscience Australia) | Australia only | NBART: BRDF- and terrain-corrected surface reflectance | Public COGs, no auth, `final` datasets only |

DEA builds its product from ESA's L1C with its own atmospheric, BRDF and terrain correction, so its values differ from L2A for the same acquisition, and it has no Sen2Cor outputs. NBART's BRDF correction normalises to a fixed sun angle, so the gap is seasonal: over a Perth AOI the visible bands were within 7% of AWS in January but 25-45% brighter in June-July. Don't mix DEA and L2A mosaics in one analysis without accounting for this. On DEA, `cloud_mask="SCL"` reads DEA's fmask layer instead (same 20 m grid, classes translated to their SCL equivalents), and the `visual`, `SCL`, `AOT`, `WVP` and `B09` bands are unavailable. DEA's COGs only have overviews at 8x and coarser, so any request between native resolution and 80 m downloads full-resolution data: slow for whole tiles or large areas (S2Mosaic warns), barely noticeable for small bounds and AOIs. Use native resolution or 80 m and coarser.

```python
array, profile = mosaic(
    grid_id="50HMH",
    start_year=2023,
    duration_months=2,
    source="DEA",
    cloud_mask="SCL",  # fmask on DEA
)
```

See [Example use - DEA.ipynb](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/examples/Example%20use%20-%20DEA.ipynb) for a DEA walkthrough: a first mosaic, a side-by-side with AWS, fmask against OmniCloudMask, and a full tile.

## Try in Colab

[![Colab_Button]][Link]

[Link]: https://colab.research.google.com/github/DPIRD-DMA/S2Mosaic/blob/main/examples/Quick%20start.ipynb 'Try S2Mosaic In Colab'

[Colab_Button]: https://img.shields.io/badge/Try%20in%20Colab-grey?style=for-the-badge&logo=google-colab

The Colab badge opens [`examples/Quick start.ipynb`](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/examples/Quick%20start.ipynb), a minimal end-to-end demo.

## Installation

S2Mosaic 2.0 is currently in beta. Install with pip:
```
pip install --pre s2mosaic
```
Or with uv:
```
uv add --prerelease=allow s2mosaic
```
Drop `--pre` / `--prerelease=allow` once `2.0.0` ships stable.

## Quick start: MGRS grid tile

Mosaic an entire Sentinel-2 MGRS tile by ID, returning a NumPy array and rasterio profile. Find the tile ID for your area of interest with the [Sentinel-2 grid explorer](https://dpird-dma.github.io/Sentinel-2-grid-explorer/).

```python
from s2mosaic import mosaic

array, profile = mosaic(
    grid_id="50HMH",  # Sentinel-2 MGRS tile ID
    start_year=2022,
    start_month=1,
    duration_months=2,  # collect scenes over this window
    scene_order="valid_data",  # prioritise scenes with the most valid pixels
    mosaic_method="mean",  # combine valid pixels by mean
    bands=["B04", "B03", "B02", "B08"],
)

print(f"Mosaic array shape: {array.shape}")
print(f"CRS: {profile['crs']}")
```

To save a GeoTIFF instead of returning the array, pass `output_dir=Path("output")`
for an auto-generated filename, or `output_path=Path("output/custom.tif")` to
choose the exact filename. The function then returns the file path.

Auto-generated `output_dir` filenames include a readable summary of the target,
date range, bands, method, scene order, resolution, cloud-mask provider, source, and a short
deterministic hash of the output-affecting request fields, so similar requests do
not silently collide while the filename stays scan-friendly. Either way, auto-generated
or explicit, a matching `.json` sidecar is written next to the GeoTIFF with the
normalized request metadata.

## Quick start: arbitrary bounding box

Pass `bounds=(minx, miny, maxx, maxy)` instead of `grid_id` to mosaic any rectangular AOI, including ones that cross MGRS tile boundaries. Each intersecting scene is streamed through a rasterio `WarpedVRT` and aggregated onto a common UTM grid.

```python
from s2mosaic import mosaic

# ~7.6km x 3.4km AOI near Perth, WA, in UTM zone 50S (EPSG:32750)
bounds = (389410, 6462290, 397010, 6465700)

array, profile = mosaic(
    bounds=bounds,
    input_crs=32750,
    snap_to_source_grid=True,
    start_year=2023,
    start_month=6,
    duration_months=2,
    bands=["B04", "B03", "B02"],
    mosaic_method="median",
)

print(f"Shape: {array.shape}")
print(f"CRS:   {profile['crs']}")
print(f"Pixel: {profile['transform'].a}m")
```

`bounds=` always fills the requested rectangle: same-CRS uses it directly, cross-CRS (e.g. lon/lat input → UTM output) uses its reprojected axis-aligned envelope. There is no implicit polygon mask, so cross-CRS bounds don't produce nodata wedges at the corners; the envelope is just slightly larger than the original lon/lat region. The recommended pattern for single-zone AOIs is to pass bounds in the local UTM zone so `input_crs == output_crs` (no envelope inflation) with `snap_to_source_grid=True` (zero resampling against the source grid). Use lon/lat input for AOIs that genuinely span multiple UTM zones, or pass `aoi=shapely.geometry.box(*bounds)` instead if you want the lat/lon rectangle clipped after reprojection. `output_crs` defaults to the UTM zone containing the AOI centroid if omitted. Use `resolution` (in metres) and `resampling_method` (`nearest`, `bilinear`, ...) to control the output grid. See [Example use - bounds.ipynb](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/examples/Example%20use%20-%20bounds.ipynb) for cross-tile and lower-resolution examples, [Example use - aoi.ipynb](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/examples/Example%20use%20-%20aoi.ipynb) for polygon AOIs, and [Advanced - wide area visual export.ipynb](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/examples/Advanced%20-%20wide%20area%20visual%20export.ipynb) for a large, wide visual-only GeoTIFF export that streams tiles directly to disk.

## Mosaic method comparison

Different mosaic methods trade speed, smoothness, and spectral consistency. In this cloudy Singapore example, `medoid` keeps the sharper real-scene look of an observed pixel while avoiding the most scene-order-dependent artefacts from `first` and the blended look of `mean` / `median`.

![S2Mosaic method comparison showing first valid, mean, median, and medoid mosaics over a cloudy Singapore AOI](https://raw.githubusercontent.com/DPIRD-DMA/S2Mosaic/main/assets/mosaic-method-comparison.jpg)

The image was generated by [Mosaic method comparison.ipynb](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/examples/Mosaic%20method%20comparison.ipynb).

## Advanced usage

Every `mosaic()` parameter, with defaults shown in parentheses.

### Area (pass exactly one)

- `grid_id` (`None`): Sentinel-2 MGRS tile ID, e.g. `"50HMH"`. Mosaics the entire tile.
- `bounds` (`None`): `(minx, miny, maxx, maxy)` rectangle. Mosaics an arbitrary AOI, including ones that cross MGRS tile boundaries. See the bounds/AOI-mode-specific options below.
- `aoi` (`None`): single shapely `Polygon`. Mosaics the polygon bounds while skipping and masking pixels outside the polygon. Mutually exclusive with `grid_id` and `bounds`.

### Time window

- `start_year` (required), `start_month` (`1`), `start_day` (`1`): start of the scene-search window.
- `duration_years` (`0`), `duration_months` (`0`), `duration_days` (`0`): length of the search window. Inclusive of `start_*`, exclusive of the end.

### Output content

- `bands` (`["B04", "B03", "B02", "B08"]`): bands to include. Leave as `None` to use the default RGB+NIR set. Use `["visual"]` for the 3-band uint8 TCI RGB composite (mutually exclusive with other bands).
- `mosaic_method` (`"mean"`): how per-pixel scene stacks are reduced to one output value.
  - `"mean"`: per-band arithmetic mean of all valid scenes. Streams scenes incrementally so peak memory stays small.
  - `"first"`: first valid pixel in `scene_order`. The cheapest method. Reads only what's needed to fill each tile and stops as soon as it can.
  - `"percentile"`: per-band percentile across all valid scenes. Requires the `percentile` parameter (0-100).
  - `"median"`: shortcut for `"percentile"` with `percentile=50`.
  - `"medoid"`: picks the scene whose multi-band spectrum is closest (squared Euclidean) to the per-band median across all valid scenes for the pixel. Preserves real observed spectra, so band relationships stay coherent for indices and classifiers. This is the approximate-medoid formulation used by [LandTrendr on Google Earth Engine](https://github.com/eMapR/LT-GEE/blob/master/LandTrendr.js) ([Kennedy et al. 2018](https://doi.org/10.3390/rs10050691)) and the [Open-MRV](https://open-mrv.readthedocs.io/en/latest/image_composite_Web.html) tutorials, not the strict [Flood 2013](https://doi.org/10.3390/rs5126481) pairwise-distance medoid; the two often agree, but can differ.
- `percentile` (`None`): percentile to compute when `mosaic_method="percentile"` (0-100).
- `min_observations` (`None`): minimum valid observations to read per pixel for `"mean"`, `"percentile"`, and `"medoid"`. When set, tile aggregation stops reading later scenes once every coverable pixel has reached the target. This is not an output quality guarantee; pixels that cannot reach the target use whatever observations are available.
- `max_observations` (`None`): per-pixel cap on valid observations for `"mean"`, `"percentile"`, and `"medoid"`. Each pixel accepts at most this many scenes (in `scene_order`); later valid scenes are dropped for that pixel. Combined with `scene_order="oldest"` or `"newest"` this biases the mosaic toward early or late dates. Must be `>= min_observations` when both are set; ignored by `"first"` (effectively N=1).
- `include_observation_count` (`False`): append a final `Observation count` band with the number of valid source scenes that contributed to each output pixel. Works for both returned arrays and exported GeoTIFFs. For `bands=["visual"]`, enabling this writes/returns `uint16` output so the count band is not limited to 255; the RGB values remain in their usual 0-255 range.
- `tile_workers` (`8`): number of output tiles to aggregate concurrently. Tuned higher than CPU count because the work is I/O-bound on remote COG reads. Raise for faster networks; lower if memory or simultaneous-connection limits matter more than throughput.
- `adaptive_tiling` (`True`): split sparse output tiles based on the actual cloud-valid contribution masks. This reduces wasted reads for irregular AOIs, sparse coverage, and heavily masked scenes. Set to `False` to use fixed-size output tiles.

### Output destination

- `output_dir` (`None`): if set, writes a GeoTIFF to this directory using an auto-generated filename and returns the file path. The filename includes a readable request summary plus a short deterministic hash, and a matching `.json` sidecar records the normalized request metadata. Mutually exclusive with `output_path`.
- `output_path` (`None`): if set, writes a GeoTIFF to this exact `.tif`/`.tiff` path and returns it. The `.json` sidecar is written here too, alongside the GeoTIFF. Mutually exclusive with `output_dir`.
- `overwrite` (`True`): when exporting and the target path exists, controls whether to overwrite it.

### Output grid

- `output_crs` (`None`): EPSG of the output. Must be a projected CRS. Geographic CRSes (e.g. 4326) are rejected at validation, because `resolution` is metres in the target CRS and a geographic output would produce a degenerate grid. If you need a lat/lon raster, reproject the mosaic afterwards with `gdalwarp` / `rio warp`. In bounds/AOI mode, auto-picked as the UTM zone containing the AOI centroid if omitted. For AOIs wider than ~6° of longitude (one UTM zone), pass an explicit equal-area projection instead (e.g. `output_crs=3577` for Australia, `5070` for the contiguous US), because the auto-picked centroid UTM has growing scale distortion and a larger envelope overshoot far from its central meridian. Ignored in grid mode (the tile's native UTM zone is used).
- `resolution` (`10`): output pixel size in metres. Coarser resolutions read from the coarsest COG overview that is no coarser than the output pixel, in every mode, so much less data crosses the wire. Overviews are a fixed set of levels (MPC 2x-32x, AWS 2x-16x, DEA only 8x-32x); when a request falls well between them (at least 4x the pixels needed, e.g. 60 m from DEA's 10 m bands, which has no overview under 80 m), S2Mosaic logs a warning naming the nearest efficient resolutions. Each warning is shown once per session, so a loop of `mosaic()` calls logs it once.
- `resampling_method` (`"nearest"`): how the source is resampled to the output grid. Also accepts `"bilinear"`, `"cubic"`, `"average"`, `"lanczos"`.
- `snap_to_source_grid` (`False`): bounds/AOI mode only. When `True`, expand the output extent outward to whole multiples of `resolution` in the target CRS. This makes repeat runs over the same area produce identical grids, and at `resolution=10` aligns the output to the native Sentinel-2 pixel grid, so source COG reads become zero-cost copies rather than sub-pixel resamples. The output may grow by up to one pixel on each side; pixels outside an `aoi` polygon are still written as nodata.

### Scene selection

- `source` (`"MPC"`): STAC provider. `"MPC"` (default) uses Microsoft Planetary Computer with SAS-signed URLs. `"AWS"` uses Element 84's Earth Search on AWS Open Data: Sentinel-2 L2A scenes, public COGs, no auth, no SAS rotation. `"DEA"` uses Digital Earth Australia's Sentinel-2 NBART (Australia only; see [Choosing a source](#choosing-a-source)). All return spectral bands on one scale, `reflectance * 10000`: from processing baseline 04.00 (January 2022) ESA adds 1000 to every spectral DN, Element 84 removes it before publishing, and S2Mosaic removes it from MPC scenes at read time. A 04.00+ DN at or below 1000 (reflectance <= 0, e.g. dark water) reads as 1, not 0, so it stays a valid observation; 0 is always nodata.
- `additional_query` (`{"eo:cloud_cover": {"lt": 100}}`): extra STAC query filters in the Query extension's form, e.g. `{"eo:cloud_cover": {"lt": 80}}`. DEA's API only accepts CQL2, so for `source="DEA"` the filter is translated; the operators `eq`, `neq`, `lt`, `lte`, `gt`, `gte` and `in` are supported, and any other raises instead of being dropped. Property names are the provider's own (DEA's cloud cover is fmask's, as `eo:cloud_cover`).
- `min_coverage_fraction` (`None`): optional scene-edge trimming. When set, drops pixels covered by fewer than this fraction of the maximum scene-overlap count in the requested area. The default keeps the full requested coverage.
- `ignore_duplicate_items` (`True`): drop duplicate acquisitions, keeping the latest processing baseline.

### Scene ordering

- `scene_order` (`"valid_data"`): scene ordering. One of `"valid_data"`, `"oldest"`, or `"newest"`.
- `scene_sort_fn` (`None`): custom callable `fn(items: pd.DataFrame) -> pd.DataFrame`. Overrides `scene_order` when set.

### Cloud masking

- `cloud_mask` (`"OCM"`): mask provider. `"OCM"` runs the OmniCloudMask deep-learning model on R+G+NIR bands (most accurate); `"SCL"` reads the L2A Scene Classification Layer (much cheaper, lower accuracy). On `source="DEA"`, which has no SCL, `"SCL"` reads DEA's fmask layer.
- `ocm_batch_size` (`1`): OCM inference batch size. Only used with `cloud_mask="OCM"`.
- `ocm_inference_dtype` (`"fp32"`): OCM inference dtype. Defaults to `"fp32"`, which runs everywhere and is the fastest option on CPU. On GPU, use `"fp16"` for ~2× speedup and lower VRAM, or `"bf16"` on hardware that supports it. Only used with `cloud_mask="OCM"`.

### Diagnostics

- `show_progress` (`False`): show tqdm progress bars for the cloud-mask streaming and tile-aggregation phases. Useful in notebooks; leave off for headless/batch runs.

Example:

```python
array, profile = mosaic(
    bounds=(389410, 6462290, 397010, 6465700),
    input_crs=32750,
    snap_to_source_grid=True,
    start_year=2023,
    duration_months=2,
    bands=["visual"],
    mosaic_method="mean",
    include_observation_count=True,
)

rgb = array[:3]
observation_count = array[3]
print(array.shape)  # (4, height, width): Red, Green, Blue, Observation count
```

### Bounds/AOI-mode-specific options

- `input_crs` (`4326`): EPSG of `bounds` or `aoi`.

For the exact function signature and return types, see the `mosaic()` docstring in the source code.

## Logging

S2Mosaic emits progress logs at each pipeline stage (search, sort, fetch, cloud-mask, aggregate, export). Following standard Python logging convention, no output is produced unless logging is configured. The simplest way to enable it:

```python
import s2mosaic

s2mosaic.set_log_level("INFO")  # or "DEBUG" for more detail
```

If your application already configures the `logging` module, the package logger (`s2mosaic`) will respect that, and there is no need to call `set_log_level()`.

## Performance tips

- `cloud_mask`: Default `"OCM"` runs the OmniCloudMask deep-learning model, which is the most accurate but needs reasonable compute (GPU/MPS recommended). Switch to `"SCL"` on CPU-only machines or for bulk processing; it skips inference entirely and just reads the L2A Scene Classification Layer.
- `ocm_batch_size`: If using a GPU, setting this above the default value (1) will speed up cloud masking. In most cases, a value of 4 works well. If you encounter CUDA errors, try using a lower number.
- `ocm_inference_dtype`: defaults to `'fp32'`, which runs on every backend and is the fastest option on CPU (most CPUs don't have efficient fp16/bf16 paths). On GPU, switch to `'fp16'` for ~2× faster inference and lower VRAM use, or `'bf16'` on hardware that supports it (Ampere+ NVIDIA, Apple Silicon).
- `scene_order`: Using `"valid_data"` tends to work well with early stopping because clear scenes are considered first.
- `min_observations`: For large `"mean"`, `"percentile"`, or `"medoid"` jobs, set this to the number of observations per pixel you actually need to avoid reading later scenes for already-satisfied tiles.
- `max_observations`: Caps each pixel at N valid scenes. Combine with `scene_order="oldest"` (or `"newest"`) to bias the mosaic toward early/late dates over a long search window without paying for the extra reads.
- `resolution`: Data read scales with the source level each read comes from, not just the output size. Resolutions that land on a COG overview (native, or 2x, 4x, 8x... of it; DEA only has 8x and coarser) read the least; ones between levels read the finer level and downsample. S2Mosaic warns when that costs 4x or more.
- `mosaic_method`: Roughly ordered fastest to slowest, `"first"` < `"mean"` < `"percentile"`/`"median"` ≈ `"medoid"`. `"first"` only reads pixels needed to fill each tile and stops as soon as it can, so cloud-free scenes can finish a tile in one pass. `"mean"` streams every contributing scene but accumulates incrementally, so its memory stays small. `"percentile"`/`"median"` and `"medoid"` both hold each tile's per-scene stack in memory to compute the result, so they use more RAM and, without `min_observations`/`max_observations`, read every contributing scene. Both keep that stack as `uint16` plus a separate validity mask, so it is the same size for either method. `"medoid"` still peaks lower because its kernel returns `uint16` directly and stripe-blocks its scratch arrays, while percentile/median returns `float32` and is clipped back afterwards: on a 12-scene, 4-band, 2048px tile the stack is 403 MB and peak is roughly 690 MB for median against 540 MB for medoid. Peak scales with scene count (the same tile over 34 scenes peaks at about 1.5 GB) and every tile worker pays it, so `tile_workers` multiplies it. Set `min_observations` (and/or `max_observations`) to cap reads once every coverable pixel has enough samples.

## Known limitations

- **SCL is less accurate than OCM.** The L2A Scene Classification Layer is fast (one COG read per scene, no inference) but is consistently less accurate than OCM at identifying clouds and cloud shadow. Use SCL when compute is the bottleneck (CPU-only machines, bulk processing); use OCM when accuracy matters.

- **Microsoft Planetary Computer serves isolated single-band zeros over dark water.** At roughly 0.01-0.04% of valid pixels, one band reads `0` (NODATA) where its neighbours read near-zero reflectance. On the harmonised scale (see `source`) the value is off by only a few DN, but a zero in one band still distorts ratios and spectral indices at those pixels, where the denominators are already small. A pixel is only discarded when *every* requested band reads 0, so these survive: discarding them instead would throw away the good bands to avoid the bad one, and leave a hole no other scene could fill. Element 84's copy of the same acquisitions is unaffected, so `source="AWS"` avoids it entirely.

- **`bands=["visual"]` cannot distinguish black water from no data.** TCI is a quantised 8-bit render, so near-zero reflectance rounds to 0, and ESA also reserves 0 for NODATA. A pixel whose three channels all round to 0 is therefore indistinguishable from an unobserved one and is dropped - about 0.23% of valid pixels on a Perth AOI, all inside water that renders near black anyway. Request the spectral bands instead where dark water matters, since there DN 0 really is NODATA.

- **`source="DEA"` lags about two weeks and drops a few scenes.** Only `final` datasets are searched, and a final version usually appears 9-13 days after acquisition (up to ~50), so a window ending in the last couple of weeks has fewer scenes than on MPC or AWS. A small number of acquisitions (1-3 per tile over April-September 2026 on three sampled tiles) are only ever published as provisional `nrt`, and are never included. DEA also publishes no `visual` asset; request `["B04", "B03", "B02"]` and stretch for display.

## Contributing

Contributions are welcome. Open an issue or a pull request.

### Running the tests

Tests use `pytest`. The development environment uses the Python pinned in `.python-version` (3.14); CI runs the suite on every supported version, 3.10-3.14. The fast suite (unit tests + mocked pipelines) takes well under a minute and is what CI runs by default:

```bash
uv run pytest                       # full fast suite
uv run pytest tests/test_readers.py # one file
uv run pytest -k requeue            # match by name
```

End-to-end tests that hit the network and run a real mosaic are marked `slow` and excluded by default (see `addopts` in `pyproject.toml`). To run them explicitly:

```bash
uv run pytest -m slow               # only slow tests
uv run pytest -m ""                 # everything, including slow
```

The example notebooks run against live STAC sources and take several minutes, so they are not part of the per-push CI. A separate workflow executes them weekly, on pull requests that change `examples/`, and on demand from the Actions tab. To run them locally:

```bash
uv run pytest --nbmake examples/*.ipynb
```

Lint with ruff:

```bash
uv run ruff check s2mosaic/ tests/
```

For maintainers: the release flow (cut a tag, GitHub Actions builds + publishes to PyPI) is documented in [RELEASING.md](https://github.com/DPIRD-DMA/S2Mosaic/blob/main/RELEASING.md).

## License

This project is licensed under the MIT License.

## Acknowledgments

S2Mosaic is built on top of:

- [Sentinel-2](https://sentiwiki.copernicus.eu/web/s2-products): ESA's Copernicus Earth-observation mission, the imagery source.
- [Element 84 Earth Search](https://earth-search.aws.element84.com/): optional public AWS Open Data access to Sentinel-2 L2A COGs.
- [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/): the default STAC catalog and signed access to the Sentinel-2 L2A archive.
- [Digital Earth Australia](https://www.dea.ga.gov.au/): Geoscience Australia's Sentinel-2 NBART analysis-ready data and fmask, used by `source="DEA"`.
- [OmniCloudMask](https://github.com/DPIRD-DMA/OmniCloudMask): the deep-learning cloud and cloud-shadow mask used by the default `cloud_mask="OCM"` provider.
- L2A Scene Classification Layer (SCL): the published per-scene classification used by the optional `cloud_mask="SCL"` provider.
- [rasterio](https://rasterio.readthedocs.io/), [GeoPandas](https://geopandas.org/), [pystac-client](https://pystac-client.readthedocs.io/), [OpenCV](https://opencv.org/), [Numba](https://numba.pydata.org/), and [multiclean](https://github.com/DPIRD-DMA/multiclean): supporting libraries for I/O (including per-scene `WarpedVRT` reprojection), geometry, search, image ops, percentile aggregation, and mask post-processing.
