satimg

A uniform, tile-friendly wrapper over satellite imagery. Point it at a product directory and you get the same interface regardless of sensor: two views on the pixels and one clean way to walk them in windows without ever losing track of where a window came from.

Supported products: Sentinel-1 GRD (IW / EW), Sentinel-2 L1C, Landsat C2 L1.

Concepts

.raw and .visual are plain lazy xarray.DataArrays, dims (band, y, x) — use them with the full xarray API directly.

Opening a product

satimg.open takes a product directory or a zip archive — detected by content, not by extension — and returns the matching Product either way.

A directory gives you full access, same as ever:

import satimg

product = satimg.open("/data/S2A_MSIL1C_20220114T103401_..._T33UUB_....SAFE")

product.raw                    # -> xarray.DataArray
for patch in product.patches(512):
    ...

A zip archive is read zip-native, straight off the archive — nothing is ever extracted to disk, and there's no temp dir to release, so a with block isn't needed just for a one-off read:

product = satimg.open("/data/S2A_MSIL1C_20220114T103401_..._T33UUB_....SAFE.zip")

product.timestamp
product.footprint
product.transformer
product.metadata               # all sensors except Landsat's per-pixel angles
product.thumbnail()

product.raw                    # raises ZipNativeUnsupportedError -- no
                                # pixels without extracting

For pixel access (raw / visual / patches / patches_at) from a zip, extract it first with satimg.open_zip:

with satimg.open_zip("/data/scene.zip") as product:
    for patch in product.patches(512):
        ...

This returns a fully capable, directory-backed product exactly like satimg.open on a directory would. It extracts to a temp dir that is removed — along with the product's open rasters — when the block exits; read what you need inside the block, since a lazy view can't be rebuilt afterwards. Pass dest= to control where that temp dir is created.

A product's open rasters (and, for a zip-native one, its open archive) are released when its with block exits, so use with satimg.open(path) as product: — zipped or not — if you loop over many products, to release handles promptly rather than waiting on garbage collection.

Two views on the pixels

import satimg

product = satimg.open("/data/S2A_MSIL1C_20220114T103401_..._T33UUB_....SAFE")

product.raw     # every native band, merged onto one grid, native dtype  -> xarray.DataArray
product.visual  # uint8 visualisation: 3-band true colour for optical    -> xarray.DataArray
                # sensors (Sentinel-2, Landsat), 1-band greyscale for SAR
                # (Sentinel-1)

Both are lazy DataArrays over the same grid, so a Window means the same thing in each. The band axis is labelled, so select by name:

product.raw.sel(band="red")        # Landsat / (SWIR etc.)  ; "B04" for Sentinel-2, "VV" for Sentinel-1
product.raw.band.values            # e.g. ['coastal', 'blue', 'green', 'red', ...]
product.raw.wavelength_nm          # Sentinel-2 only: central wavelength per band

Walking the image in patches

for patch in product.patches(512, overlap=64):
    tile = patch.raw.values        # np.ndarray, (band, 512, 512), native dtype
    x0, y0 = patch.col, patch.row  # exact offset in the full image
    ...

for patch in product is shorthand for product.patches() (raw bands, 512 px, no overlap). Every patch from product.patches() / patches_at() already carries .raw / .visual / .meta for that window (see Reading raw + visual + metadata together) -- product.patches() / patches_at() are the only iterators in the library.

product.patches() / patches_at() open raw and visual chunked to exactly the window size before iterating, so each patch read decodes just the one tile it covers.

A product releases its rasters when its with block exits, so looping over many products never leaks file handles:

with satimg.open(path) as product:
    ...

Edge handling

product.patches(size, edge=...) decides what happens at the right / bottom border:

Indices are always conserved: patch.window tells you where the tile belongs, patch.window.clip(product.width, product.height) gives its valid region.

Patches at specific points

When you already know where to look — detections, sample sites, a sparse set of tiles — patches_at visits only those points instead of sweeping the whole image. Each (col, row) is the centre of its patch, and every patch is exactly size (overhang zero-filled, so points near an edge still work):

points = [(4096, 2048), (10500, 512), (0, 0)]
for patch in product.patches_at(points, 512):
    tile = patch.raw.values         # (band, 512, 512), centred on the point

satimg.windows_at(points, size) gives the bare Windows if you don't need a patch bound.

Reading raw + visual + metadata together

Every patch from product.patches() / product.patches_at() already carries its own window's raw, visual and metadata -- one window, every view.

scene = product.metadata.attrs                 # scene-level scalars, constant per product
                                               # e.g. {"mean_sun_zenith": 77.5, ...}

for p in product.patches_at(points, 512):
    raw = p.raw                                # xarray.DataArray (band, 512, 512), band-labelled
    red = p.raw.sel(band="red").values         # np.ndarray (512, 512)   ("B04" on Sentinel-2)
    nir = p.raw.sel(band="nir08").values       #                          ("B08" on Sentinel-2)
    rgb = p.visual.values                       # np.ndarray (3, 512, 512) uint8, same window

    angles = p.meta.sample()                    # {"sun_zenith": 78.1, "view_zenith": 4.9, ...}
                                               #   -> per-pixel fields, bilinear at the patch centre
    sza_grid = p.meta["sun_zenith"].values      # np.ndarray (512, 512) — the field over the patch

Mapping results back

patch = next(iter(product.patches(1024)))
boxes = detect(patch.visual.values)     # boxes in patch pixels
boxes[:, [0, 2]] += patch.window.col    # -> full-image pixels
boxes[:, [1, 3]] += patch.window.row
latlon = product.transformer.rowcol_to_latlon(boxes[:, [1, 0]])

Working with a single window

from satimg import Window

win = Window(col=4096, row=2048, width=1024, height=1024)
chip = satimg.read_window(product.raw, win)   # lazy xarray.DataArray, (band, y, x)
chip.values                                   # numpy, zero-padded if the window overhangs

Per-pixel metadata

Every sensor ships scene geometry that varies across the image — Sentinel-1's incidence / elevation angles, Sentinel-2's and Landsat's sun / viewing angles — stored on a coarse grid. product.metadata exposes those as named fields that bilinearly interpolate to any pixel:

m = product.metadata
m.fields                          # ['sun_zenith', 'sun_azimuth', 'view_zenith', ...]
m.sun_zenith.at((row, col))       # -> float, bilinear on the coarse grid
m.sample((row, col))              # -> {field: value} for every field
m.sun_zenith.corners()            # -> {'top_left', 'top_right', 'bottom_left', 'bottom_right', 'center'}
m.corners()                       # -> {field: {'top_left': ..., ..., 'center': ...}}  (whole grid)
m.attrs                           # scalar scene-level metadata (mean angles, ...)

This is per-pixel geometry (varies across the scene). Per-band identity lives on product.raw instead — the labelled band coordinate, plus wavelength_nm for Sentinel-2. Scene-wide scalars are product.metadata.attrs.

Each field also materialises as a lazy full-grid (y, x) xarray.DataArray (m.sun_zenith.grid), and patches from product.patches() carry a matching lazy view over their window:

for patch in product.patches(512):
    patch.meta.sun_zenith          # lazy (512, 512) DataArray for this window
    patch.meta.at((y, x))          # -> {field: value} at a patch-local pixel
    patch.meta.sample()            # -> {field: value} at the patch centre
    patch.meta.corners()           # -> {field: {'top_left': ..., ..., 'center': ...}}

corners() samples five points (the four corners and the centre) without materialising the patch grid — a cheap way to see how a field varies across a tile. corners()[field]["center"] equals sample()[field].

Geometry without the product

Every product's geometry -- its transformer, image shape, coarse metadata grids and scene attrs -- is a small dict that survives strict JSON (e.g. a PostgreSQL jsonb column; NaN is stored as null). Store it while the product is open and keep converting coordinates and sampling metadata after the file is gone. For Sentinel-1, load it as a Sentinel1Geometry to also get the SAR geometry (Doppler shift, headings):

with satimg.open("S1A_IW_GRDH_....SAFE.zip") as product:
    stored = product.geometry.to_dict()      # ~30-80 KB

geometry = satimg.Sentinel1Geometry.from_dict(stored)   # SceneGeometry for other products
geometry.transformer.latlon_to_rowcol((lat, lon))
geometry.metadata.incidence_angle.at((row, col))
geometry.doppler_azimuth_shift((row, col), speed_mps, heading_deg)   # Sentinel-1 only

Landsat's angle fields are stored resampled to 32 x 32, and, like its metadata, need the product extracted (open_zip).

Downloading

search is anonymous; download needs a (free) account for the archive. Pass credentials to the connector, or set environment variables (CDSE_USERNAME / CDSE_PASSWORD for Copernicus, USGS_USERNAME / USGS_PASSWORD / optional USGS_TOKEN for USGS).

from satimg import get_connector

connector = get_connector("sentinel-2-l1c")                     # search only
items = connector.search(geojson=aoi, start_datetime=start, stop_datetime=stop)

connector = get_connector("sentinel-2-l1c",
                          username="me@example.com", password="...")
connector.download(items[0], "/data/scene.zip")

download always gives you a zip archive — see Opening a product for reading its metadata straight off the archive, or extracting it for pixel access.

Deep Zoom

from satimg.deepzoom import DeepZoom

DeepZoom(tile_size=512).build(product.visual, "/tmp/scene")  # -> /tmp/scene.dzi

Adding a product

Subclass Product, implement the abstract hooks (_open_raw(tile), _render_visual(raw, tile) — both return an xarray.DataArray dask-chunked to tile; normalise dims with satimg.as_band_yx(...) and name the band axis with satimg.label_bands(da, [...]) — transformer, timestamp, footprint, thumbnail), optionally override _read_metadata to expose per-pixel metadata, and register a filename pattern:

from satimg.registry import register

@register(r"^MY_SENSOR_\d{8}.*$")
class MySensorProduct(Product):
    ...

Importing satimg.products registers the built-ins.

Logging

satimg logs under the satimg.* logger namespace and never installs a handler of its own. Opt in from your application:

import logging

logging.basicConfig(level=logging.INFO)              # your app configures handlers
logging.getLogger("satimg").setLevel(logging.DEBUG)  # DEBUG for per-tile / per-request detail

INFO covers milestones (downloads, Deep Zoom builds, detection runs); DEBUG adds STAC queries, redirect hops, per-tile counts and product resolution. Credentials are never logged.

Development

uv sync
uv run pytest                      # fast unit tests + product integration tests
uv run pytest tests/test_geometry.py tests/test_tiling.py   # just the fast ones

Integration tests run against minified scenes under tests/data/ — real georeferencing and metadata, zeroed pixels.