---------------------------------------------------------------------- This is the API documentation for the s1bursts library. ---------------------------------------------------------------------- ## Searching and downloading Query the CDSE burst catalog and fetch single bursts over OData or S3. SearchQuery(time_range: 'TimeRange | None' = None, intersects: 'HasWkt | None' = None, burst_id: 'int | None' = None, abs_burst_id: 'int | None' = None, polarization: 'Polarization | None' = None, mode: 'Mode | None' = None, sub_swath: 'SubSwath | None' = None, pass_direction: 'PassDirection | None' = None, platform_serial_id: 'PlatformSerialID | None' = None, rel_orbit_number: 'int | None' = None) -> 'None' Query parameters for searching Sentinel-1 bursts via the CDSE API. All parameters are optional; omitting a parameter means no filtering on that field. Parameters ---------- time_range Acquisition time window. The search keeps every acquisition on any of its days, end included, and the same window can then name the composite built from the search. intersects Spatial filter; bursts must intersect this geometry, in WGS 84 longitude and latitude. A geometry that carries a ``crs``, such as an ``odc.geo`` ``Geometry``, must be in EPSG:4326. The catalog rejects a geometry with too many vertices, such as the exact outline of an MGRS cell; search with its bounding box and select with ``BurstCollection.intersecting()``. burst_id Sentinel-1 relative burst ID, unique within each sub-swath. Must be positive. abs_burst_id Sentinel-1 absolute burst ID, globally unique across all sub-swaths and orbits. Must be positive. polarization Polarization channel. mode Acquisition mode. sub_swath Sub-swath identifier. pass_direction Pass direction. platform_serial_id Sentinel-1 platform serial identifier. rel_orbit_number Relative orbit number (1-175). Raises ------ TypeError If ``time_range`` is not a ``TimeRange``. ValueError If any value falls outside its accepted range or set of identifiers, or ``intersects`` is not a geometry, is empty or is not in WGS 84. TimeRange(start: 'DateLike', end: 'DateLike') -> 'None' An inclusive window of calendar days, ``start`` through ``end``. The one object a pipeline hands to both the search and the composite: :class:`~s1bursts.query.SearchQuery` selects the acquisitions that fall on any of its days, and ``ResolutionWeighter.compose_*()`` names and tags the product with it. Because both take the same window, the product's dates describe exactly what was searched. Days have no time of day, so there is nothing to be exclusive about: a window ``2023-07-18`` through ``2023-10-15`` covers every acquisition on October 15. Filters that compare instants get the same span from :meth:`half_open`. Parameters ---------- start, end First and last day of the window. An instant names the UTC calendar day it falls on, the resolution every date tag is written at. Raises ------ ValueError If ``end`` is before ``start``. BurstCatalog(file_manager: 'FileManager', auth: 'CDSEAuth | None' = None, *, show_progress: 'bool' = True, continue_on_error: 'bool' = True, default_endpoint: 'Endpoint' = 'odata') -> 'None' Search and download Sentinel-1 bursts from the CDSE catalog. Provides search via the CDSE API and parallel download with progress tracking, skip-existing logic, and automatic endpoint fallback. Parameters ---------- file_manager File manager for determining SLC download paths. auth Authentication credentials for the CDSE API. If ``None``, a ``CDSEAuth`` instance is created lazily from the ``CDSE_USERNAME`` / ``CDSE_PASSWORD`` settings on the first download. Searching the catalog requires no authentication. show_progress Whether to show progress bars during downloads. continue_on_error Whether to continue downloading after individual failures. ``True``, the default, attempts every burst and records each failure in the summary, so one missing burst does not end a long download; ``False`` raises at the first burst that fails on both endpoints, once the downloads already running have finished, without starting the queued ones. A burst this catalog already found missing at the source is not requested again: it is recorded as failed without raising, whichever the setting, and a rerun on the same catalog reports it again rather than retrying it. Follow a run with ``OperationSummary.raise_on_failure()`` to fail once all bursts have been attempted. default_endpoint Default download endpoint. If a download fails, the other endpoint is tried automatically. Raises ------ ValueError If ``default_endpoint`` is neither ``"odata"`` nor ``"s3"``. ## Burst data model The acquisitions a search returns, and the interferometric pairs built from them. Burst(burst_id: 'int', sub_swath: 'SubSwath', polarization: 'Polarization', date: 'Timestamp', abs_burst_id: 'int | None' = None, cdse_id: 'str | None' = None, parent_product_name: 'str | None' = None, _rel_orbit_number: 'int | None' = None, _pass_direction: 'PassDirection | None' = None, _geometry: 'Polygon | MultiPolygon | None' = None) -> None Represents a single Sentinel-1 burst. A burst is a segment of SAR data acquisition characterized by its spatial location, temporal acquisition time, and imaging parameters. Two bursts are considered equal when their ``burst_id``, ``sub_swath``, ``polarization``, and calendar day (``date`` normalized to midnight) match, so a burst from a CDSE search compares equal to the same acquisition reconstructed from an SLC filename (which carries only the day). A given footprint is re-acquired at most once per repeat cycle, so the day is a unique key. BurstPair(primary: 'Burst', secondary: 'Burst') -> None Represents a pair of Sentinel-1 bursts for interferometric processing. The primary burst is typically the earlier acquisition, and the secondary is the later acquisition. The pair must have matching burst parameters. BurstCollection(df: 'gpd.GeoDataFrame | None' = None) -> 'None' Collection of Sentinel-1 bursts with geospatial and temporal data. Provides methods for working with multiple bursts, including finding interferometric pairs, filtering, and visualization. Parameters ---------- df DataFrame with burst data. If ``None``, an empty collection is created. Raises ------ TypeError If df is not a GeoDataFrame. ValueError If required columns are missing. BurstPairCollection(df: 'gpd.GeoDataFrame | None' = None) -> 'None' Collection of Sentinel-1 burst pairs for interferometric processing. Each entry contains a primary and secondary burst with matching spatial parameters but different acquisition times. ## Processing Batch SNAP GPT execution and the runtime budget it runs under. BurstProcessor(file_manager: 'FileManager', config: 'SNAPConfig | None' = None, *, show_progress: 'bool' = True, continue_on_error: 'bool' = True) -> 'None' Batch processor for Sentinel-1 bursts and burst pairs. Runs one SNAP GPT job per burst or pair with progress tracking, skips products that already exist with the requested settings, records failures per item, and returns an ``OperationSummary``. Takes bursts and burst pairs, as collections or lists, rather than product classes. Parameters ---------- file_manager File manager resolving input SLC and output product paths. config Runtime configuration for SNAP GPT processing; defaults to ``SNAPConfig()``. show_progress Whether to show progress bars. continue_on_error Whether to continue processing after an error. ``True``, the default, attempts every item and records each failure in the summary, so one bad burst does not end a long job; ``False`` raises at the first error, once the items already running have finished, without starting the queued ones. Follow a run with ``OperationSummary.raise_on_failure()`` to fail once all items have been attempted. SNAPConfig(bin: 'Path | None' = None, rasterio_profile: 'dict[str, Any] | None' = None, max_workers: 'int | None' = None, max_ram_gb: 'float | None' = None, threads_per_job: 'int | None' = None, cache_ratio: 'float | None' = None, auxdata_path: 'Path | None' = None, timeout: 'float | None' = ) -> None Runtime configuration for SNAP GPT processing. Parameters ---------- bin A SNAP installation's ``bin`` directory, used to locate the ``gpt`` executable. Defaults to ``SNAP_BIN``, then ``gpt`` discovered on ``PATH``. rasterio_profile Overrides for the rasterio profile used when rewriting GPT output GeoTIFFs. Merged over ``{"dtype": "float32", "nodata": "nan"}``. max_workers Number of parallel GPT subprocesses; at least 1. Defaults to ``SNAP_MAX_WORKERS`` (or 1). max_ram_gb Total RAM budget across all GPT processes in GB; must be positive and large enough to give every worker a 512 MB heap on top of 3 GB of JVM and OS overhead. Defaults to ``SNAP_MAX_RAM_GB`` (or 8). threads_per_job Threads per GPT process (``-q``); at least 1. Defaults to ``SNAP_THREADS_PER_JOB`` or ``cpu_count // max_workers``. cache_ratio Fraction of the per-process Java heap used for the GPT tile cache (``-c``); greater than 0 and at most 1. Defaults to ``SNAP_CACHE_RATIO`` (or 0.7). auxdata_path Shared SNAP auxdata cache (DEM tiles, orbit files). Defaults to ``SNAP_AUXDATA_PATH`` (or ``~/.snap/auxdata``). timeout Seconds one GPT run may take before it is stopped and recorded as failed; must be positive and finite, or ``None`` for no limit. Defaults to ``SNAP_TIMEOUT`` (or 3600), where ``none`` means no limit. Attributes ---------- resolved_gpt_path The ``gpt`` executable this configuration runs, from ``bin``, ``SNAP_BIN`` or ``PATH``. resolved_snap_version Version of the SNAP installation that ``gpt`` belongs to, such as ``(14, 0, 0)``, or ``None`` if it cannot be determined. Raises ------ ValueError If a resource argument or ``timeout`` is out of range, or if the RAM budget cannot cover the workers and ``max_ram_gb`` or ``max_workers`` was passed. ConfigurationError If a ``SNAP_*`` environment variable is invalid, if the RAM budget cannot cover the workers and both values come from the environment, or if GPT or the auxdata path cannot be resolved. SNAPClient(file_manager: 'FileManager', config: 'SNAPConfig | None' = None) -> 'None' Low-level client for SNAP GPT processing. Provides direct interface to SNAP GPT for processing Sentinel-1 bursts. For batch processing with progress tracking, use BurstProcessor instead. Parameters ---------- file_manager File manager resolving input SLC and output product paths. config Runtime configuration for SNAP GPT; defaults to ``SNAPConfig()``. ## Graph configuration The shape of the generated SNAP graph, one config per product type. Every config carries the shared Terrain-Correction tail parameters. RTCConfig(*, grid_spacing: Annotated[float, Gt(gt=0)] = 10.0, dem: str | pathlib.Path = 'Copernicus 30m Global DEM', dem_name: str | None = None, img_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', dem_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION', 'CUBIC_CONVOLUTION', 'BISINC_5_POINT_INTERPOLATION', 'BISINC_11_POINT_INTERPOLATION', 'BISINC_21_POINT_INTERPOLATION', 'BICUBIC_INTERPOLATION', 'DELAUNAY_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', align_to_standard_grid: bool = True, external_dem_nodata: float = 0.0, external_dem_apply_egm: bool = False, replace_zeros_with_nan: bool = True, multilook: tuple[Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])], Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])]] | None = (5, 1), dem_oversampling: Annotated[float, Gt(gt=0)] = 1.0, additional_overlap: Annotated[float, Ge(ge=0)] = 0.2) -> None Configuration for terrain-flattened backscatter (RTC). Gamma-0 backscatter, radiometrically terrain corrected: called Normalised Radar Backscatter (NRB) in CEOS-ARD and RTC-S1 by OPERA. The product also carries the local scattering area that terrain flattening divides by. Inherits every parameter of ``BackscatterConfig``. Parameters ---------- dem_oversampling DEM oversampling factor used by ``Terrain-Flattening``. Must be positive. additional_overlap Additional overlap percentage used by ``Terrain-Flattening`` when searching for terrain contributions outside the tile. GTCConfig(*, grid_spacing: Annotated[float, Gt(gt=0)] = 10.0, dem: str | pathlib.Path = 'Copernicus 30m Global DEM', dem_name: str | None = None, img_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', dem_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION', 'CUBIC_CONVOLUTION', 'BISINC_5_POINT_INTERPOLATION', 'BISINC_11_POINT_INTERPOLATION', 'BISINC_21_POINT_INTERPOLATION', 'BICUBIC_INTERPOLATION', 'DELAUNAY_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', align_to_standard_grid: bool = True, external_dem_nodata: float = 0.0, external_dem_apply_egm: bool = False, replace_zeros_with_nan: bool = True, multilook: tuple[Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])], Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])]] | None = (5, 1)) -> None Configuration for geocoded backscatter (GTC). Beta-0 backscatter, geometrically terrain corrected only: geocoded, but not terrain-flattened. Takes exactly the parameters of ``BackscatterConfig``. CoherenceConfig(*, grid_spacing: Annotated[float, Gt(gt=0)] = 10.0, dem: str | pathlib.Path = 'Copernicus 30m Global DEM', dem_name: str | None = None, img_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', dem_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION', 'CUBIC_CONVOLUTION', 'BISINC_5_POINT_INTERPOLATION', 'BISINC_11_POINT_INTERPOLATION', 'BISINC_21_POINT_INTERPOLATION', 'BICUBIC_INTERPOLATION', 'DELAUNAY_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', align_to_standard_grid: bool = True, external_dem_nodata: float = 0.0, external_dem_apply_egm: bool = False, replace_zeros_with_nan: bool = True, window: tuple[typing.Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])], typing.Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])]] = (25, 5), subtract_flat_earth_phase: bool = True, subtract_topographic_phase: bool = True, correct_coherence_bias: bool = False, remove_local_phase_ramp: bool = False) -> None Configuration for the coherence pipeline. Inherits every parameter of ``GeocodingConfig``. No ``Multilook`` step is applied, unlike ``BackscatterConfig``. It would save no runtime, since it can only sit after the coregistration and phase-removal stages that dominate this graph, and placing it after the ``Coherence`` operator averages coherence magnitudes and so cannot reduce the small-window estimator bias. Use ``window`` to set the looks. Parameters ---------- window Coherence estimation window ``(n_rg, n_az)`` in single-look samples. This is the only control over the number of looks entering the estimator. SNAP's ``squarePixel`` option, which derives the azimuth window from the range window and the scene geometry, is pinned off, so the window given is the window applied. subtract_flat_earth_phase Remove the flat-earth phase before estimating coherence. subtract_topographic_phase Remove the topographic phase using ``dem``. correct_coherence_bias Remove the positive bias a finite estimation window leaves in the coherence magnitude. Requires SNAP 14; building the graph for an older or undetermined version raises ``ConfigurationError``. Off by default, unlike SNAP's own default, so that the estimator stays the one the window describes and stays comparable across SNAP versions. remove_local_phase_ramp Let SNAP estimate and remove a local phase ramp within each estimation window before averaging, which recovers coherence lost to residual fringes. Requires SNAP 14, as above. BackscatterConfig(*, grid_spacing: Annotated[float, Gt(gt=0)] = 10.0, dem: str | pathlib.Path = 'Copernicus 30m Global DEM', dem_name: str | None = None, img_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', dem_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION', 'CUBIC_CONVOLUTION', 'BISINC_5_POINT_INTERPOLATION', 'BISINC_11_POINT_INTERPOLATION', 'BISINC_21_POINT_INTERPOLATION', 'BICUBIC_INTERPOLATION', 'DELAUNAY_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', align_to_standard_grid: bool = True, external_dem_nodata: float = 0.0, external_dem_apply_egm: bool = False, replace_zeros_with_nan: bool = True, multilook: tuple[Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])], Annotated[int, FieldInfo(annotation=NoneType, required=True, metadata=[Gt(gt=0)])]] | None = (5, 1)) -> None Settings shared by the backscatter pipelines, RTC and GTC. Abstract: pass ``RTCConfig`` for terrain-flattened gamma-0 backscatter or ``GTCConfig`` for beta-0 backscatter that is geocoded only. Inherits every parameter of ``GeocodingConfig``. Parameters ---------- multilook ``(n_rg_looks, n_az_looks)`` applied after ``TOPSAR-Deburst``, before terrain flattening; ``None`` skips the ``Multilook`` step. Multilooking here cuts runtime sharply, because the reduced raster passes through both terrain flattening and terrain correction. SNAP's ``grSquarePixel`` option, which derives the azimuth looks from the range looks and the local geometry and discards the requested ``n_az_looks`` is pinned off. Therefore, the looks given are the looks applied. GeocodingConfig(*, grid_spacing: Annotated[float, Gt(gt=0)] = 10.0, dem: str | pathlib.Path = 'Copernicus 30m Global DEM', dem_name: str | None = None, img_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', dem_resampling: Literal['NEAREST_NEIGHBOUR', 'BILINEAR_INTERPOLATION', 'CUBIC_CONVOLUTION', 'BISINC_5_POINT_INTERPOLATION', 'BISINC_11_POINT_INTERPOLATION', 'BISINC_21_POINT_INTERPOLATION', 'BICUBIC_INTERPOLATION', 'DELAUNAY_INTERPOLATION'] = 'BILINEAR_INTERPOLATION', align_to_standard_grid: bool = True, external_dem_nodata: float = 0.0, external_dem_apply_egm: bool = False, replace_zeros_with_nan: bool = True) -> None Terrain-Correction tail parameters shared by every product type. Abstract: inherited by each per-product config, which supplies the graph that precedes the shared Terrain-Correction tail via ``to_xml``. Parameters ---------- grid_spacing Output pixel spacing in meters of the target CRS. Must be positive. dem DEM used for geocoding. A string selects one of SNAP's auto-download DEMs; a path (or a string with a file suffix) is passed as an external DEM. dem_name Name of an external DEM, written to product metadata and compared before an existing product is skipped. Required when ``dem`` is a path, and must not be blank; ignored for SNAP's built-in DEMs, whose ``dem`` value is already their name. Name each DEM differently, such as by its source and version. img_resampling Resampling of the image during terrain correction. SNAP's cubic and sinc kernels are not supported, because they corrupt the pixels next to missing data. dem_resampling Resampling of the DEM during terrain correction. align_to_standard_grid Snap the output to a fixed global grid, so products from different bursts and dates share pixel boundaries and can be stacked without resampling. external_dem_nodata No-data value of an external DEM. external_dem_apply_egm Convert an external DEM from EGM96 geoid heights to WGS84 ellipsoidal heights, which is the reference terrain correction requires. SNAP does this automatically for its built-in ACE / Copernicus / SRTM DEMs, all of which are geoid-referenced; this flag is the equivalent for an external one and is ignored unless ``dem`` is a path. Set it ``True`` for an orthometric DEM (heights above mean sea level, the usual case for national DEMs) and leave it ``False`` only when the DEM already carries ellipsoidal heights. Note this default differs from SNAP's own, which is ``True``. Leave it ``False`` for a DEM returned by ``s1bursts.prepare_dem_for_snap``, which has already been converted to ellipsoidal heights. SNAP only applies EGM96, so preprocess DEMs in other vertical datums such as EGM2008 instead of enabling this option. replace_zeros_with_nan Convert SNAP's zero fill value to NaN, so no-data does not enter later averages as a valid low value. ## Storage Deterministic product paths and discovery, on local disk or S3. FileManager(base_dir: 'str | Path | UPath', slc_dir: 'str | Path | UPath | None' = None, prod_dir: 'str | Path | UPath | None' = None, comp_dir: 'str | Path | UPath | None' = None) -> 'None' Manages file system operations for Sentinel-1 products. All filename generation logic is centralized here. Storage directories (``base_dir``, ``slc_dir``, …) may point to either local paths or S3 locations via ``upath.UPath``. A directory given as a string is read as a local path, with ``~`` expanded. Parameters ---------- base_dir Root directory for all products. May be a local ``Path`` or a ``UPath`` pointing to an S3 bucket. slc_dir Directory for SLC burst files. Defaults to ``base_dir/s1-slc``. prod_dir Directory for processed products. Defaults to ``base_dir/s1-prod``. comp_dir Directory for composites. Defaults to ``base_dir/s1-comp``. Raises ------ ValueError If a directory is given as an empty string, or as a URL string rather than a ``UPath``. ## Compositing Locally resolution weighted and per-pass composites on the 10 km squares of the MGRS grid, re-exported from freezebase. ResolutionWeighter(file_manager: 'FileManager', *, rasterio_profile: 'dict[str, Any] | None' = None, max_workers: 'int' = 8, checksum: 'bool' = False) -> 'None' Compute LRW and per-pass composites for Sentinel-1 products. Each compose call receives its target ``geobox``, so one instance can be reused across MGRS grid cells. Parameters ---------- file_manager Used to resolve input file paths and output product paths. max_workers Number of threads for parallel file I/O during product loading. rasterio_profile GeoTIFF profile applied to every product, or profiles keyed by product type with an optional ``"default"`` entry. Product-specific settings override the default. If omitted, uses the library defaults with float32 data and NaN no-data. checksum Compute a SHA-256 multihash of every written file's exact bytes and record it in ``checksums``, keyed by output path. Costs holding the encoded COG in memory a second time per write. Attributes ---------- checksums : dict[str, str] SHA-256 multihash of every file written so far, keyed by output path as a string. Populated only when ``checksum=True``. References ---------- .. [1] Small, D., Rohner, C., Miranda, N., Ruetschi, M., & Schaepman, M. E. (2022). Wide-Area Analysis-Ready Radar Backscatter Composites. IEEE Transactions on Geoscience and Remote Sensing, 60, 5201814. https://doi.org/10.1109/TGRS.2021.3055562 .. [2] CEOS-ARD Product Family Specification: Synthetic Aperture Radar, version 1.3.1 (2026). Committee on Earth Observation Satellites. https://ceos.org/ard/files/PFS/SAR/v1.3.1/CEOS-ARD_PFS_SAR_v1.3.1.pdf MGRSGrid(filter_geometry: 'GeoSeries | GeoArray | BaseGeometry', resolution: 'float' = 10.0) -> 'None' The 10 km MGRS squares intersecting an area of interest. Adjacent UTM zones may overlap because each zone is gridded independently. Parameters ---------- filter_geometry : GeoSeries | GeoArray | BaseGeometry Area of interest in WGS84. Intersecting squares are included whole. resolution : float, default 10.0 Pixel size in metres. Must divide 10 km evenly. Raises ------ ValueError If ``filter_geometry`` carries a CRS other than WGS84. MGRSGeoBox(mgrs_code: 'str', shape: 'SomeShape', affine: 'Affine', crs: 'CRS') -> 'None' A 10 x 10 km MGRS square backed by ``odc.geo.GeoBox``. Parameters ---------- mgrs_code : str MGRS reference. shape : SomeShape Raster dimensions. affine : Affine Pixel-to-CRS transform. crs : CRS Cell CRS. ## Single-burst products Metadata records returned by the processing services. They describe a product; the pixels live in the GeoTIFF at its `FileManager` path. RTCProduct(burst_id: 'int', sub_swath: 'SubSwath', crs: 'CRSLike', date: 'Timestamp', polarization: 'Polarization', *, pass_direction: 'PassDirection', _rel_orbit_number: 'int | None' = None) -> None Radiometrically Terrain Corrected (RTC) backscatter product. Gamma nought backscatter with terrain flattening applied. Parameters ---------- burst_id Burst identifier. sub_swath Sub-swath identifier. crs The coordinate reference system. date Acquisition date. polarization The polarization. pass_direction Pass direction, ``"Ascending"`` or ``"Descending"``. Keyword-only. GTCProduct(burst_id: 'int', sub_swath: 'SubSwath', crs: 'CRSLike', date: 'Timestamp', polarization: 'Polarization', *, pass_direction: 'PassDirection', _rel_orbit_number: 'int | None' = None) -> None Geometrically Terrain Corrected (GTC) backscatter product. Beta nought backscatter with geometric terrain correction only. Parameters ---------- burst_id Burst identifier. sub_swath Sub-swath identifier. crs The coordinate reference system. date Acquisition date. polarization The polarization. pass_direction Pass direction, ``"Ascending"`` or ``"Descending"``. Keyword-only. CoherenceProduct(burst_id: 'int', sub_swath: 'SubSwath', crs: 'CRSLike', date_primary: 'Timestamp', date_secondary: 'Timestamp', polarization: 'Polarization', *, pass_direction: 'PassDirection', _rel_orbit_number: 'int | None' = None) -> None Coherence product from burst pair. Parameters ---------- burst_id Burst identifier. sub_swath Sub-swath identifier. crs The coordinate reference system. date_primary Primary acquisition date. date_secondary Secondary acquisition date. polarization The polarization. pass_direction Pass direction, ``"Ascending"`` or ``"Descending"``. Keyword-only. ## Composite products MGRS grid-square composites. Each holds its measurement and, optionally, its quality-assurance layers. RTCLRWComposite(mgrs_code: 'str', start_date: 'Timestamp', end_date: 'Timestamp', polarization: 'Polarization', *, layers: 'frozenset[Layer]' = frozenset({})) -> None Locally resolution weighted RTC intensity composite. Only RTC (radiometrically terrain corrected) products can be combined using LRW, as they account for terrain effects. The stored value is intensity (power), i.e. the squared amplitude. Parameters ---------- mgrs_code MGRS grid cell identifier. start_date Start of the compositing period (inclusive). end_date End of the compositing period (inclusive). polarization The polarization. layers The layers the product holds: the measurement, ``QA-NUM`` and ``QA-CQM``. Keyword-only; defaults to the measurement alone. CohLRWComposite(mgrs_code: 'str', start_date: 'Timestamp', end_date: 'Timestamp', polarization: 'Polarization', baseline_days: 'int', *, layers: 'frozenset[Layer]' = frozenset({})) -> None Locally resolution weighted coherence composite. Parameters ---------- mgrs_code MGRS grid cell identifier. start_date Start of the compositing period (inclusive). end_date End of the compositing period (inclusive). polarization The polarization. baseline_days Temporal baseline in days. layers The layers the product holds: the measurement, ``QA-NUM`` and ``QA-CQM``. Keyword-only; defaults to the measurement alone. RTCPassComposite(mgrs_code: 'str', start_date: 'Timestamp', end_date: 'Timestamp', polarization: 'Polarization', pass_direction: 'PassDirection', *, layers: 'frozenset[Layer]' = frozenset({})) -> None Unweighted temporal mean composite for a single pass direction. Parameters ---------- mgrs_code MGRS grid cell identifier. start_date Start of the compositing period (inclusive). end_date End of the compositing period (inclusive). polarization The polarization. pass_direction Pass direction. layers The layers the product holds: the measurement and ``QA-NUM``. Keyword-only; defaults to the measurement alone. CohPassComposite(mgrs_code: 'str', start_date: 'Timestamp', end_date: 'Timestamp', polarization: 'Polarization', pass_direction: 'PassDirection', baseline_days: 'int', *, layers: 'frozenset[Layer]' = frozenset({})) -> None Unweighted temporal mean coherence composite for a single pass direction. Parameters ---------- mgrs_code MGRS grid cell identifier. start_date Start of the compositing period (inclusive). end_date End of the compositing period (inclusive). polarization The polarization. pass_direction Pass direction. baseline_days Temporal baseline in days. layers The layers the product holds: the measurement and ``QA-NUM``. Keyword-only; defaults to the measurement alone. GTCPassComposite(mgrs_code: 'str', start_date: 'Timestamp', end_date: 'Timestamp', polarization: 'Polarization', pass_direction: 'PassDirection', *, layers: 'frozenset[Layer]' = frozenset({})) -> None Unweighted temporal composite of GTC beta0 for a single pass direction. Parameters ---------- mgrs_code MGRS grid cell identifier. start_date Start of the compositing period (inclusive). end_date End of the compositing period (inclusive). polarization The polarization. pass_direction Pass direction. layers The layers the product holds: the measurement and ``QA-NUM``. Keyword-only; defaults to the measurement alone. Layer(*values) One file of a composite product: its measurement or one of its quality layers. A quality layer adds its value to the composite's file name as a token, such as ``S1_RTC_LRW_QA-NUM_…``; the measurement file carries no token. ## Results and errors What a batch operation reports, and the exception hierarchy that separates per-item failures from environment-wide ones. OperationSummary(n_total: 'int', n_processed: 'int', n_skipped: 'int', n_failed: 'int', products: 'list[T]', failed_items: 'list[ErrorResult]', processing_time: 'float' = 0.0) -> None Summary of batch operation results (processing, downloading, etc.). Parameters ---------- n_total Total number of operations attempted. n_processed Number of successfully processed items. n_skipped Number of skipped operations (already existed). n_failed Number of failed operations. products The products successfully processed or already existed. failed_items List of operation results for failed items. processing_time Total operation time in seconds. OkResult(source: 'Any', result: 'T', skipped: 'bool' = False) -> None Successful or skipped operation result. Parameters ---------- source The source burst or burst pair. result The resulting product. skipped Whether the operation was skipped (product already exists). ErrorResult(source: 'Any', error: 'Exception') -> None Failed operation result. Parameters ---------- source The source burst or burst pair. error The exception that caused the failure. S1BurstsError Base class for all s1bursts errors. OperationFailedError(summary: 'OperationSummary[Any]', ignore: 'tuple[type[BaseException], ...]' = ()) -> 'None' Items of a batch download or processing run failed. Raised by ``OperationSummary.raise_on_failure()`` once the batch has finished, so it is never a per-item failure itself. The message lists the failed items and their errors. The error carries the whole summary, so the products of the items that succeeded stay within reach. Parameters ---------- summary Summary of the batch, with at least one failed item. ignore Error types whose failures did not raise this error. They are counted in the message but not listed. Attributes ---------- summary Summary of the batch: ``summary.products`` holds the products created or found, ``summary.failed_items`` every failure, ignored ones included. ignore Error types whose failures did not raise this error. ConfigurationError Invalid or missing configuration or environment. Raised for problems that affect every item — missing credentials, an unusable PROJ/GDAL install, a missing SNAP GPT executable. These propagate and fail the run fast; batch handlers must not catch them. TransformationGridError A PROJ grid needed for a DEM's vertical conversion is missing or could not be downloaded. BurstMapUnavailableError The Sentinel-1 burst-ID map is needed but not downloaded or unreadable. CatalogError(*args: 'object', status_code: 'int | None' = None) -> 'None' A CDSE catalog request failed, or its response cannot be used. Raised when the catalog cannot be reached, answers with an HTTP error, for example HTTP 414 for a search geometry with too many vertices, or returns a malformed response. The underlying ``requests`` error is kept as the cause. Parameters ---------- *args The message, as for any exception. status_code HTTP status of the catalog's response. Attributes ---------- status_code HTTP status of the catalog's error response, such as ``414`` or ``503``. ``None`` if the catalog could not be reached or answered successfully with an unusable response. TooManyResultsError More bursts match a search than one catalog request returns. The catalog serves at most 1000 bursts per request. A query that matches more is refused rather than paged through, so that an oversized AOI or time window fails after one request. Narrow it, or split it into several searches and combine their results with ``+``. BurstNotInMapError A burst's ``(burst_id, sub_swath)`` is not in the Sentinel-1 burst-ID map. Every burst the mission acquires is in the map, so this points to a malformed burst rather than one missing at the source, which is ``BurstNotFoundError``. ProductSettingsMismatchError(*args: 'object', path: 'Path | UPath | None' = None, differences: 'Mapping[str, tuple[str | None, str | None]] | None' = None) -> 'None' An existing product was made with other settings than requested. Also raised for a per-burst product made by another SNAP version than the running one, unless either version is undetermined, and for a product that does not record every requested setting, such as one written by an earlier s1bursts. Raised instead of skipping the product when ``overwrite=False``: by ``BurstProcessor.create_*`` and ``SNAPClient.create_*`` for a per-burst product, and by ``ResolutionWeighter.compose_*`` for a composite. Pass ``overwrite=True`` to replace the product, or give the other settings their own ``FileManager`` directories to keep both. Parameters ---------- *args The message, as for any exception. path The existing file whose settings differ. differences The differing settings, as for the attribute. Attributes ---------- path The existing file whose settings differ. differences Each differing setting, keyed by the tag or grid property compared, such as ``"GRID_SPACING"`` or ``"TRANSFORM"``, mapped to its recorded and its requested value; for ``"PROCESSING_SOFTWARE_VERSION"``, the requested value is the running SNAP version. ``None`` stands for a tag that is absent. BurstError Base class for recoverable, per-item burst failures. Batch loops catch this type (not bare ``ValueError``/``OSError``) to record an ``ErrorResult`` and continue with the remaining items. BurstNotFoundError A requested burst or its parent product is not available at the source. BurstDownloadError Downloading or extracting a single burst failed. NoUTMZoneError A burst's footprint lies outside the UTM zones, so ``crs="auto"`` has none for it. Raised by ``Burst.utm_zone`` and ``BurstPair.utm_zone`` for a footprint center north of 84 degrees or south of 80 degrees. With ``crs="auto"`` it belongs to that burst alone, so ``BurstProcessor.create_*`` and ``SNAPClient.create_*`` record it for the item rather than stopping. SnapProcessingError(*args: 'object', returncode: 'int | None' = None, output_tail: 'str | None' = None) -> 'None' Processing a single burst or burst pair with SNAP failed. Raised when the GPT graph fails, and for an error staging the input SLCs or writing the product, which is kept as the cause. Parameters ---------- *args The message, as for any exception. returncode Exit code of ``gpt``, negative for the signal that ended it. output_tail The last lines ``gpt`` printed. Attributes ---------- returncode Exit code of ``gpt``, negative for the signal that ended it, such as ``-9`` for ``SIGKILL``. ``None`` if ``gpt`` timed out or did not run, or the failure came from staging or writing. output_tail The last lines ``gpt`` printed to stdout and stderr, or ``None`` if it printed nothing or did not run. CompositeError Compositing one grid cell failed; another cell may still succeed. Raised by ``ResolutionWeighter.compose_*`` for problems with the data of the cell, never for an invalid argument, which raises ``ValueError``, or for inputs that cannot be composited together, which raise ``InconsistentInputsError``. EmptyCompositeError No input observation has a valid pixel in the target cell. A burst whose footprint touches a cell need not reach any of its pixels once geocoded. A composite of such inputs would hold no value, so none is written, and the cell is left without one, as if no burst had matched. InconsistentInputsError The input products passed to ``compose_*`` cannot be composited together. They are on another CRS, pixel size or pixel grid than the target cell, were made with different DEMs, processing settings or SNAP versions, or differ in polarization or temporal baseline. An undetermined SNAP version matches any. The inputs are chosen by the call, so every other cell given such inputs fails alike: this is not a ``CompositeError``, and a loop over cells stops at it. CompositeInputError An input product of a cell is missing or cannot be read. Missing includes the RTC product that supplies an observation's local scattering area in a coherence composite: an RTC that failed to process leaves only the cells its burst reaches without one. ## Reference data and low-level access The optional SAR-MPC burst-ID map, external-DEM preparation, and the CDSE clients the high-level services are built on. download_burst_map(cache_dir: 'str | Path | None' = None, *, mode: 'Mode' = 'IW') -> 'Path' Download and prepare the Sentinel-1 burst-ID map for a given acquisition mode. The map is a large (~500 MB per mode) static dataset from SAR-MPC. It is only needed to look up the relative orbit number, pass direction, or footprint geometry of a burst that was **not** obtained from a CDSE search or an SLC filename (e.g. a manually constructed burst). It is not downloaded automatically; call this once to opt in. The raw download is verified against a pinned SHA-256. Parameters ---------- cache_dir Root directory for cached dataset files. Defaults to the ``FREEZEBASE_CACHE`` environment variable or the system cache directory. mode Sentinel-1 acquisition mode. Returns ------- Path Path to the prepared GeoParquet file. Raises ------ BurstMapUnavailableError If the download fails, the downloaded archive is corrupt or holds no burst-ID map for ``mode``, or a map prepared earlier is unreadable. A failed download keeps the ``requests`` or freezebase error as the cause; a corrupt archive or unreadable prepared map keeps the read error and names the file to delete. Sentinel1Bursts(cache_dir: 'str | Path | None' = None, *, mode: 'Mode' = 'IW') -> 'None' Sentinel-1 Burst ID map. Sentinel-1 performs systematic acquisition of bursts in both IW and EW modes. The bursts overlap almost perfectly between different passes and are always located in the same place. With the deployment of the SAR processor S1-IPF 3.4, a new element has been added to the products annotations: the Burst ID, which should help the end-user to identify a burst area of interest and facilitate searches. **Important**: The Burst ID is unique in each sub-swath. For this reason to identify a burst univocally it is necessary to specify both the "*subswath_name*" (see below) and the "*burst_id*". Two maps are provided: one for EW and one for IW acquisition modes. The maps have a validity that covers the entire time span of the mission and they are global, i.e., they include as well information where no SAR data is acquired. Parameters ---------- cache_dir The root directory where the dataset is stored. Defaults to the 'FREEZEBASE_CACHE' environment variable or the system cache directory (in that order). mode The Sentinel-1 acquisition mode. prepare_dem_for_snap(dem: 'str | Path', *, source_crs: 'str | None' = None, download_missing_grids: 'bool' = False, cache_dir: 'Path | None' = None, snap_bin: 'Path | None' = None) -> 'Path' Return an external DEM prepared for SNAP's sampling and height conventions. The output geotransform is shifted by half a pixel because SNAP interprets an external DEM's origin as the center of its first pixel, and so samples the file half a pixel south-east of where its georeferencing places it [1]. SNAP fixed that in the 14.0.2 update of its ``snap-dem`` module, so the shift is applied only for an installation older than that, and for one whose version cannot be determined. When ``source_crs`` is supplied, elevations are also transformed from that compound horizontal and vertical CRS to WGS 84 ellipsoidal heights (EPSG:4979). SNAP cannot do that conversion itself: the only geoid it knows is EGM96, which it applies even to the Copernicus DEM, whose heights are referenced to EGM2008 [2]. The result is correct **only** when read by the SNAP installation it was prepared for, and must only be passed as an ``externalDEMFile``. Set SNAP's ``externalDEMApplyEGM`` parameter to ``False``. Prepared copies are cached by source content and conversion settings. Parameters ---------- dem Local DEM to prepare. source_crs Compound CRS describing the DEM's horizontal and orthometric references, for example ``"EPSG:4326+3855"`` for Copernicus GLO-30. If omitted, elevations pass through unchanged. download_missing_grids Download required open-license PROJ grids when they are not installed. cache_dir Directory for prepared copies. Defaults to ``SNAP_DEM_CACHE`` (``~/.s1bursts/dem``). snap_bin A SNAP installation's ``bin`` directory, which decides whether the half-pixel shift is needed. Defaults to ``SNAP_BIN``, then ``gpt`` discovered on ``PATH``. Returns ------- Path Path to the prepared copy. Raises ------ FileNotFoundError If ``dem`` does not exist. ValueError If the source CRS is invalid or inconsistent with the DEM. TransformationGridError If no exact vertical transformation is available, because a PROJ grid it needs is missing or could not be downloaded. pyproj.exceptions.ProjError If any DEM pixel falls outside the selected transformation's valid area. ConfigurationError If ``snap_bin`` is given but contains no ``gpt`` executable. References ---------- .. [1] External DEM is sampled half a pixel off. ESA STEP forum. https://forum.step.esa.int/t/external-dem-is-sampled-half-a-pixel-off/46042 .. [2] Copernicus DEM built-in path applies EGM96 while the data is EGM2008. ESA STEP forum. https://forum.step.esa.int/t/copernicus-dem-built-in-path-applies-egm96-while-the-data-is-egm2008/46041 CDSEAuth(username: 'str | None' = None, password: 'str | None' = None) -> 'None' CDSE authentication handler using OpenIDConnect bearer tokens. Automatically refreshes access tokens when expired. Parameters ---------- username CDSE username used to request the initial token. If not provided, falls back to the ``CDSE_USERNAME`` environment variable. password CDSE password used to request the initial token. If not provided, falls back to the ``CDSE_PASSWORD`` environment variable. Raises ------ ConfigurationError If credentials are not provided as arguments and not found in the environment, or if CDSE rejects them (HTTP 400, 401 or 403). requests.RequestException If the token request fails otherwise. The error keeps neither the request nor the response, so that the password does not travel with it. CDSEBurstClient(auth: 'CDSEAuth | None' = None) -> 'None' Client for searching and downloading Sentinel-1 bursts from CDSE using the OData API. Implemented according to the CDSE API documentation [1]. Parameters ---------- auth Authentication instance. Required for downloading bursts. References ---------- .. [1] Sentinel-1 SLC Burst API. Copernicus Data Space Ecosystem documentation. https://documentation.dataspace.copernicus.eu/APIs/Sentinel-1%20SLC%20Burst.html CDSEResponse(context: 'str', value: 'list[dict[str, Any]]', count: 'int | None' = None) -> None Response object for Sentinel-1 burst search operations. Parameters ---------- context OData context metadata string (usually includes entity type). value List of result items returned by the search. count Total number of matches on the server (if requested). extract_burst_from_s3(save_path: 'str | Path | UPath', *, parent_product_name: 'str', polarization: 'str', sub_swath: 'str', burst_id: 'int') -> 'Path | UPath' Extract a single Sentinel-1 SLC burst from CDSE S3. All intermediate files (edited XMLs, extracted TIFF) are written to a local temporary directory. Only the final ``.zip`` is moved or uploaded to ``save_path``, which may be a local path or a remote UPath (e.g. S3). Parameters ---------- save_path Full destination path (including filename) for the ``.zip``. parent_product_name Name of the parent ``*.SAFE`` product on CDSE. polarization Polarization channel, e.g. ``"VV"``, ``"VH"``, ``"HH"``, ``"HV"``. sub_swath Subswath identifier, e.g. ``"IW1"``, ``"IW2"``, ``"IW3"``. burst_id Relative burst ID within the subswath. Returns ------- Path or UPath Path to the produced ``.zip`` file. Raises ------ ValueError If ``sub_swath`` is not a recognized subswath identifier. BurstNotFoundError If the parent product or the requested burst cannot be found. BurstDownloadError If the product's annotation, calibration, noise, RFI or manifest file is malformed. The parse error is kept as the cause. ConfigurationError If the CDSE S3 credentials are not configured. xml_to_mermaid(graph_xml: 'str') -> 'str' Convert a SNAP GPT XML processing graph (string) to a Mermaid script. Parameters ---------- graph_xml SNAP GPT graph XML document. Returns ------- str ``flowchart LR`` Mermaid diagram source. build_backscatter_graph(cfg: 'BackscatterConfig') -> 'str' Emit the SNAP GPT XML for the backscatter (RTC or GTC) pipeline. build_coherence_graph(cfg: 'CoherenceConfig', *, snap_version: 'tuple[int, ...] | None' = None) -> 'str' Single-shot SLC→Coherence→Terrain-Correction graph. ---------------------------------------------------------------------- This is the User Guide documentation for the package. ---------------------------------------------------------------------- ## Getting Started # Introduction s1bursts makes Sentinel-1 synthetic aperture radar (SAR) data easy to use from Python. It searches and downloads single-look complex (SLC) data from the [Copernicus Data Space Ecosystem](https://dataspace.copernicus.eu/) (CDSE), processes it with [ESA SNAP](https://step.esa.int/main/toolboxes/snap/) into backscatter and coherence, and combines the results into composites on a fixed grid. SLC is the Sentinel-1 product that keeps the phase of the radar signal, which coherence needs. Other products, such as ground range detected (GRD) scenes, are not supported. ## Products Sentinel-1's radar sends microwave pulses to the ground and records their echo, by day and night and through clouds. s1bursts creates two products from it: - **Backscatter** is how much of the pulse returns to the satellite. It depends on the roughness, moisture and structure of the surface. Wet snow, for example, looks much darker than dry snow. Use it to follow the state of a surface through time. - **Coherence** compares two acquisitions of the same place, from 0 (unrelated) to 1 (unchanged). It drops where the surface melted, moved, grew or was disturbed between the two dates. Use it to detect change. Backscatter and coherence can then be combined into composites of many acquisitions. [Products and bands](products.qmd) describes every product in detail. ## What is a burst? s1bursts handles Sentinel-1 data one burst at a time. A burst is the smallest part of a Sentinel-1 scene that can be processed on its own. Each scene is made of three strips side by side (five in EW mode), called sub-swaths, and each sub-swath is a sequence of overlapping bursts about 20 km long. A burst download is roughly 100 MB, compared with about 4 GB for the full scene, so a study area of a few kilometers no longer costs a whole scene. The same footprint is imaged again on every orbit repeat, which makes bursts a natural unit for time series. s1bursts identifies a burst by four fields: | Field | Meaning | |---|---| | `burst_id` | Number of the burst footprint, unique within a sub-swath | | `sub_swath` | Sub-swath of the footprint: `IW1` to `IW3` (or `EW1` to `EW5`) | | `polarization` | `VV`, `VH`, `HH` or `HV` | | `date` | Acquisition time | Together, `burst_id` and `sub_swath` name a footprint: the same patch of ground on every repeat. ## From search to composite A typical workflow has four steps, each handled by one class: | Step | Class | Result | |---|---|---| | Search | `BurstCatalog.search()` | A `BurstCollection` of catalog records | | Download | `BurstCatalog.download()` | One SLC file per burst | | Process | `BurstProcessor` | One GeoTIFF per burst or burst pair | | Composite | `ResolutionWeighter` | GeoTIFFs on 10 km MGRS grid squares | A `FileManager` decides where every file is stored, on local disk or on S3, so each step finds the output of the previous one. ## Next steps Start with [Installation](installation.qmd), then run the [Quickstart](quickstart.qmd). The [tutorials](glacier-backscatter.qmd) show complete analyses on real data. # Installation s1bursts needs a Python package, a free CDSE account, and ESA SNAP. Searching the catalog works with the package alone. Downloading needs the account, and processing needs SNAP. ## Install the package s1bursts requires Python 3.12 or newer and installs from PyPI: ::: {.panel-tabset group="installer"} ### pip ```{.bash} pip install s1bursts ``` ### uv ```{.bash} uv add s1bursts ``` ### pixi ```{.bash} pixi add --pypi s1bursts ``` ### conda s1bursts is not packaged for conda. Install it from PyPI into your conda environment: ```{.bash} conda activate my-env pip install s1bursts ``` ::: ## Install SNAP Processing runs `gpt`, the command-line tool of SNAP version 13 or 14. Install it with ESA's installer, or as a conda package: ::: {.panel-tabset group="snap"} ### ESA installer Download [SNAP](https://step.esa.int/main/download/snap-download/) from ESA and include the Microwave Toolbox, which holds the Sentinel-1 tools. The installer does not always put `gpt` on your `PATH`. If it does not, set `SNAP_BIN` to the installation's `bin` directory (see [Environment variables](environment-variables.qmd)). ### conda or pixi The [`esa-snap-s1tbx-gpt`](https://anaconda.org/channels/sarforge/packages/esa-snap-s1tbx-gpt/overview) package installs `gpt` with the radar toolboxes and puts it on your `PATH` whenever the environment is active. It is a community package, not published by ESA, available for Linux, macOS and Windows. ```{.bash} pixi workspace channel add sarforge pixi add "esa-snap-s1tbx-gpt 14.*" ``` ```{.bash} conda install -c sarforge -c conda-forge "esa-snap-s1tbx-gpt=14" ``` ::: Use one SNAP version for one study. The versions give slightly different results, so a time series that mixes them contains a small artificial step (see [SNAP processing](snap-processing.qmd#snap-versions)). ## Check that SNAP is found s1bursts uses the `gpt` in `SNAP_BIN` if that is set, and otherwise the one on your `PATH`. Check which one it finds by running this in the environment where you installed s1bursts: ```{python} from s1bursts import SNAPConfig config = SNAPConfig() print("gpt:", config.resolved_gpt_path) print("SNAP version:", config.resolved_snap_version) ``` If s1bursts cannot find `gpt`, `SNAPConfig()` raises an error that says what to set. ## Set up credentials Downloads need a free [CDSE account](https://dataspace.copernicus.eu/). s1bursts reads the username and password from the environment variables `CDSE_USERNAME` and `CDSE_PASSWORD`. [Environment variables](environment-variables.qmd) shows how to set them and lists the optional settings, including keys for a second download route. # Environment variables s1bursts reads credentials and SNAP settings from environment variables. Set them in your shell, or keep them in a `.env` file: ```{.bash filename=".env"} CDSE_USERNAME=your-account-name CDSE_PASSWORD=your-password SNAP_MAX_RAM_GB=16 ``` s1bursts does not read `.env` by itself. Run your script with `uv run --env-file .env script.py`, or call `dotenv.load_dotenv()` from [python-dotenv](https://pypi.org/project/python-dotenv/) at its start. ::: {.callout-warning} Keep `.env` out of version control, since it holds your password. ::: ## CDSE credentials Downloads need a [CDSE account](https://dataspace.copernicus.eu/), but searching needs none. The S3 keys are optional and enable a second download route, which [takes over](finding-and-downloading.qmd#when-the-burst-service-fails) when CDSE's burst service fails: | Variable | Purpose | |---|---| | `CDSE_USERNAME`, `CDSE_PASSWORD` | CDSE account, used by the burst service | | `CDSE_S3_ACCESS_KEY_ID`, `CDSE_S3_SECRET_ACCESS_KEY` | CDSE S3 keys, from the [S3 key manager](https://eodata-s3keysmanager.dataspace.copernicus.eu/) | | `CDSE_S3_PROFILE` | Instead of the two keys: a profile in `~/.aws/credentials` that holds them | ## SNAP SNAP needs plenty of memory. `SNAP_MAX_RAM_GB` is the total budget shared by all parallel SNAP processes. s1bursts reserves 3 GB of it for overhead and splits the rest between the processes, so the default of 8 GB gives one process 5 GB. That is enough for backscatter and coherence. Raise the budget when you run several processes in parallel: | Variable | Purpose | Default | |---|---|---| | `SNAP_BIN` | The SNAP installation's `bin` directory | `gpt` on `PATH` | | `SNAP_MAX_WORKERS` | SNAP processes run in parallel | 1 | | `SNAP_MAX_RAM_GB` | Memory budget across those processes, in GB | 8 | | `SNAP_THREADS_PER_JOB` | Threads per SNAP process | CPU count divided by workers | | `SNAP_CACHE_RATIO` | Share of each process's memory used as tile cache | 0.7 | | `SNAP_TIMEOUT` | Seconds one SNAP process may run before it is stopped and its burst recorded as failed, or `none` for no limit | 3600 | | `SNAP_AUXDATA_PATH` | Cache for orbit files and DEM tiles | `~/.snap/auxdata` | | `SNAP_DEM_CACHE` | Cache for [prepared external DEMs](external-dems.qmd) | `~/.s1bursts/dem` | `SNAPConfig` sets the same values in code and takes precedence over the environment. See [Processing with SNAP](processing.qmd#run-several-jobs-at-once). ## Temporary files SNAP writes large uncompressed files while it works, in the folder that `TMPDIR` names. Point `TMPDIR` at a local disk with several GB free, since on many systems the default `/tmp` is held in memory. # Quickstart This example processes two acquisitions of one burst over the Aletsch glacier in the Swiss Alps. It produces terrain-corrected backscatter for each date and the coherence between them. Complete [Installation](installation.qmd) first, and run the script with your [credentials in the environment](environment-variables.qmd). The example downloads about 200 MB and runs SNAP. ```{python} #| include: false import os from pathlib import Path from dotenv import load_dotenv root = next(p for p in (Path.cwd(), *Path.cwd().parents) if (p / "pyproject.toml").exists()) load_dotenv(root / ".env") (root / "cache/docs").mkdir(parents=True, exist_ok=True) os.chdir(root / "cache/docs") ``` ## Find two acquisitions A search returns the catalog records that match a query. This one covers two weeks of one burst footprint, which Sentinel-1 imaged twice, 12 days apart: ```{python} from s1bursts import BurstCatalog, BurstProcessor, FileManager, SearchQuery, TimeRange fm = FileManager("data/quickstart") catalog = BurstCatalog(fm) bursts = catalog.search(SearchQuery( time_range=TimeRange("2022-08-01", "2022-08-13"), burst_id=30346, sub_swath="IW1", polarization="VV", )) ``` ## Download and process The download stores both acquisitions under `data/quickstart/s1-slc`. Processing then writes one GeoTIFF per product under `data/quickstart/s1-prod`. With `crs="auto"`, each burst is projected in the UTM zone that holds it, here zone 32 North (EPSG:32632). Backscatter is processed per burst. Coherence needs a pair of acquisitions, which `find_pairs()` builds from the same bursts: ```{python} #| output: false catalog.download(bursts).raise_on_failure() processor = BurstProcessor(fm) rtc = processor.create_backscatter(bursts, crs="auto").raise_on_failure() pairs = bursts.find_pairs(baseline=12) coherence = processor.create_coherences(pairs, crs="auto").raise_on_failure() ``` A burst that fails does not stop the others. `raise_on_failure()` raises an error naming every failed burst, and otherwise returns the summary of the run. Running the script again skips every file that already exists, after checking that it was made with the same settings. ## Look at the result Each product is a GeoTIFF that any GIS or raster library can open, and `fm.open()` opens one with rasterio. This plots the backscatter of the first date in decibels, with the coherence below it: ```{python} import matplotlib.pyplot as plt import numpy as np with fm.open(rtc.products[0]) as src: gamma0 = src.read(1, masked=True) with fm.open(coherence.products[0]) as src: coh = src.read(1, masked=True) fig, (top, bottom) = plt.subplots(2, 1, figsize=(7.8, 4.7), layout="constrained") fig.colorbar(top.imshow(10 * np.log10(gamma0), cmap="viridis", vmin=-20, vmax=0), ax=top, label="Backscatter (dB)") fig.colorbar(bottom.imshow(coh, cmap="magma", vmin=0, vmax=1), ax=bottom, label="Coherence") for ax in (top, bottom): ax.set_axis_off() plt.show() ``` The burst covers a strip of the Alps about 85 km long. Coherence is high (light) where the ground changed little in the 12 days, such as bare rock and meadows, and low (dark) where it changed, for example on the melting glaciers and in forests. ## What you get The processing calls return two kinds of GeoTIFF: - **Backscatter** (one file per date): how strongly the ground reflects the radar signal, corrected for terrain. - **Coherence** (one file): how similar the two acquisitions are, from 0 (no similarity) to 1. [Products and bands](products.qmd) explains each band and its units. The [glacier tutorials](glacier-backscatter.qmd) take the same burst further. ```{python} #| echo: false #| output: asis from datetime import date from importlib.metadata import version from s1bursts import SNAPConfig snap_version = ".".join(map(str, SNAPConfig().resolved_snap_version)) print( '::: {.callout-note collapse="true"}\n## Run details\n\n' f"Run on {date.today():%d %B %Y} with s1bursts {version('s1bursts')} and " f"ESA SNAP {snap_version}. Contains modified Copernicus Sentinel data (2022).\n:::" ) ``` ## How-to Guides # Finding and downloading bursts A search queries the CDSE catalog for the bursts of a place and a time window, and a download fetches them. Searching needs no account and downloads nothing, so you can refine a query freely before committing to a download. ## Search an area Pass an area of interest as a shapely geometry in longitude and latitude, together with a time window. The search returns every burst that overlaps the area and shows their footprints on a map: ```{python} from shapely.geometry import box from s1bursts import BurstCatalog, FileManager, SearchQuery, TimeRange catalog = BurstCatalog(FileManager("data")) bursts = catalog.search(SearchQuery( intersects=box(8.00, 46.40, 8.08, 46.48), time_range=TimeRange("2022-01-01", "2022-01-31"), polarization="VV", )) bursts ``` Even a small area usually returns several footprints. Bursts from different orbits and sub-swaths overlap the same few kilometers, each seeing the ground from a different direction. The area selects whole bursts and does not crop them. A `TimeRange` covers whole days, and both end dates are included. ## Follow one footprint To follow one footprint through time, search by `burst_id` and `sub_swath` instead of by area. This is the usual query once an area search has shown which footprint covers your site best: ```{.python} bursts = catalog.search(SearchQuery( time_range=TimeRange("2022-01-01", "2022-12-31"), burst_id=30346, sub_swath="IW1", polarization="VV", )) ``` `SearchQuery` also filters by orbit, pass direction and acquisition mode. Its [reference](../reference/query.SearchQuery.qmd) lists every field. ::: {.callout-note} A search returns every matching burst, up to 1000, the size of one catalog page. If a query matches more, it raises `TooManyResultsError` rather than return an incomplete list. Narrow the time window or the area, or split the query and [combine the results](#combine-searches). ::: ## Combine searches Search results add up with `+`, which drops any duplicates. Use this to collect several searches into one download, for example the same footprint in winter and in summer: ```{python} footprint = dict(burst_id=30346, sub_swath="IW1", polarization="VV") winter = catalog.search(SearchQuery(time_range=TimeRange("2022-01-01", "2022-01-31"), **footprint)) summer = catalog.search(SearchQuery(time_range=TimeRange("2022-07-01", "2022-07-31"), **footprint)) both = winter + summer print(len(winter), len(summer), len(both)) ``` ## Select part of a search A search result can be narrowed without a new query. Each selection returns a new collection, so they can be chained: | Method | Keeps the bursts | |---|---| | `within(time_range)` | Acquired on a day of a `TimeRange` | | `intersecting(geometry)` | Whose footprint touches a geometry, such as an area or an [MGRS grid square](compositing.qmd#get-the-grid-squares) | | `filter(**fields)` | Whose fields equal a value, or any value of a list | ```{.python} january = bursts.within(TimeRange("2022-01-01", "2022-01-31")) descending = bursts.filter(pass_direction="Descending", rel_orbit_number=[66, 139]) ``` ## Pair acquisitions Coherence compares two acquisitions of the same footprint. `find_pairs()` turns a search result into pairs a fixed number of days apart, which can be downloaded and [processed](processing.qmd#coherence) like single bursts: ```{.python} pairs = bursts.find_pairs(baseline=12) ``` The baseline is a whole number of days, and both arguments are passed by keyword. `max_delta` sets how far the actual spacing may deviate from it, six hours by default. A pair needs both of its acquisitions, so search a little beyond the end of the period you want and keep the pairs that start within it: ```{.python} pairs = bursts.find_pairs(baseline=12).starting_within(january) ``` ## Save a search Save a search result if you want to rerun an analysis on exactly the same acquisitions later. The records live in a GeoDataFrame, which can be written to GeoParquet and loaded again without querying the catalog: ```{.python} import geopandas as gpd from s1bursts.burst import BurstCollection bursts.df.to_parquet("inventory.parquet") bursts = BurstCollection(gpd.read_parquet("inventory.parquet")) ``` ## Download the bursts `BurstCatalog.download()` fetches the bursts of a search result, or both acquisitions of every pair. Bursts already on disk are skipped, so an interrupted download can be run again: ```{.python} summary = catalog.download(bursts, max_workers=4) summary.print_summary() ``` `max_workers` sets how many bursts download in parallel. Each file is saved in the `s1-slc` folder below the `FileManager`'s base directory. ## When the burst service fails CDSE's burst service occasionally fails for some bursts. If [CDSE S3 keys](environment-variables.qmd#cdse-credentials) are set, s1bursts then extracts the burst from the full scene on CDSE's S3 storage instead. This route is slower but usually succeeds. Without the keys, the burst is recorded as failed. Pass `endpoint="s3"` to try the S3 route first. ## Failed downloads A failed burst does not stop the download by default. It is recorded in the returned `OperationSummary` and the remaining bursts continue. Check the summary before processing: `raise_on_failure()` raises an `OperationFailedError` that lists every failed burst, and otherwise returns the summary, so it can follow the call directly: ```{.python} catalog.download(bursts).raise_on_failure() ``` To inspect the failures yourself, read them from the summary: ```{.python} for failure in summary.failed_items: print(failure.source.date, failure.error) ``` Most failures are temporary server errors and succeed on a second run, which downloads only the missing bursts. To stop at the first failure instead, create the catalog with `BurstCatalog(fm, continue_on_error=False)`. A burst that CDSE reports as missing is not requested again by the same `BurstCatalog`. Create a new one to retry it later. To accept such bursts and continue, ignore their error type: ```{.python} from s1bursts.errors import BurstNotFoundError catalog.download(bursts).raise_on_failure(ignore=BurstNotFoundError) ``` ## Reuse downloaded bursts `FileManager` lists what is already stored, without a new search. Use this to process bursts downloaded in an earlier session: ```{.python} from s1bursts import FileManager from s1bursts.burst import Burst fm = FileManager("data") bursts = fm.find(Burst, burst_id=30346, time_range=TimeRange("2022-01-01", "2022-06-30")) ``` The result is a list of bursts, which the processing methods accept as it is. To [form pairs](#pair-acquisitions), turn it into a collection first with `BurstCollection.from_bursts(bursts)`. `find()` lists every product and composite type the same way: pass the type, such as `RTCProduct` or `RTCLRWComposite` from `s1bursts.product`, and filter on the fields its filename records. These listings are rebuilt from filenames, so they contain less information than the original search result. [Save the search](#save-a-search) if you need the full catalog record. ### The burst-ID map A burst from a search includes its relative orbit, pass direction and footprint, and product and SLC filenames record the orbit and pass direction, so the usual search, download, process and composite workflow needs no further data. Only a burst that did not come from a search needs to look these up, in the [burst-ID map](https://sar-mpc.eu/sar-mpc-products/) published by SAR-MPC. That applies to a burst you construct by hand, and to the footprint of a burst found with `find(Burst)`, which its filename does not record. Download the map once, about 500 MB per acquisition mode, when you need it: ```{.python} from s1bursts import download_burst_map download_burst_map(mode="IW") ``` # Processing with SNAP `BurstProcessor` turns downloaded bursts into geocoded GeoTIFFs. It runs SNAP once for every burst, or once for every pair of bursts, and skips products that already exist. ## How a processing call works Every `create_*` method takes the bursts, the UTM zone of the output, and optionally a settings object. From the settings, s1bursts builds the SNAP processing graph, the chain of SNAP steps that turns a burst into a product. The settings describe how to process and suit every burst alike. The UTM zone depends on where the bursts are, so it is passed separately as `crs`, the coordinate reference system (CRS): ```{.python} from s1bursts import BurstProcessor, RTCConfig processor = BurstProcessor(fm) rtc = processor.create_backscatter(bursts, crs="auto", config=RTCConfig(grid_spacing=20)) ``` Products are always projected in a WGS 84 / UTM zone. `crs="auto"` takes each burst's own zone, its `utm_zone`, which suits bursts used on their own, such as a time series: the zone depends only on the footprint, so every date of a burst lands in the same one. To fix the zone instead, pass it as a string, an EPSG code or a `pyproj.CRS`, such as `crs="EPSG:32632"` for zone 32 North. If you plan to [composite](compositing.qmd), choose the zone before processing: a composite accepts only products in the UTM zone and pixel size of its grid square, so process them with `crs=square.crs`. Without `config`, the defaults apply, and they suit most uses. Each product has its own method and settings object: | Method | Input | Settings | Product | |---|---|---|---| | `create_backscatter()` | Bursts | `RTCConfig` or `GTCConfig` | Backscatter of each acquisition | | `create_coherences()` | Pairs | `CoherenceConfig` | Coherence between two acquisitions | The settings objects share most fields. The one you are most likely to change is the pixel size, `grid_spacing`, which is 10 m by default. A smaller pixel size does not add detail beyond what the radar resolves. The [API reference](../reference/snap_graph.RTCConfig.qmd) lists every setting and its default. [SNAP processing](snap-processing.qmd) explains the steps each graph runs and how to inspect the graph your settings produce. ## Backscatter Backscatter describes how strongly the ground reflects the radar signal on one date. In mountains, it also depends on the angle of each slope to the satellite: slopes facing the satellite look bright, and slopes facing away look dark. By default, s1bursts removes this effect with radiometric terrain correction (RTC), following the terrain-flattening method of Small (2011)[^small]. RTC divides each pixel by the area of ground that contributed to it, so the remaining differences in brightness come from the surface itself. The result is called gamma-0. Use RTC unless you have a reason not to. Only RTC products can be composited with [resolution weighting](resolution-weighting.qmd). The settings object selects the product: `RTCConfig`, the default, gives RTC, and `GTCConfig` gives geometrically terrain-corrected (GTC) backscatter, called beta-0, which places pixels correctly on the map but keeps the slope brightness: ```{.python} from s1bursts import GTCConfig gtc = processor.create_backscatter(bursts, crs="auto", config=GTCConfig()) ``` ## Coherence Coherence measures how similar two acquisitions of the same footprint are, from 0 to 1. It drops wherever the surface changed between the two dates, for example through snowfall, melt, vegetation growth or movement, such as the flow of a glacier. The pairs come from [`find_pairs()`](finding-and-downloading.qmd#pair-acquisitions): ```{.python} pairs = bursts.find_pairs(baseline=12) coherence = processor.create_coherences(pairs, crs="auto") ``` The product contains the coherence magnitude: the absolute value of the complex correlation between the two acquisitions, estimated over a small window of neighboring pixels set by `CoherenceConfig(window=...)`. A larger window gives smoother values at a coarser resolution. The phase of the correlation is not kept, and s1bursts makes no interferograms. ASF's [HyP3 burst InSAR](https://hyp3-docs.asf.alaska.edu/guides/burst_insar_product_guide) products cover the same bursts with wrapped and unwrapped phase. ## Run several jobs at once `SNAPConfig` sets how many SNAP processes run in parallel and how much memory they share. Pass it to the processor to override the [environment variables](environment-variables.qmd#snap): ```{.python} from s1bursts import SNAPConfig processor = BurstProcessor(fm, config=SNAPConfig(max_workers=2, max_ram_gb=16)) ``` More processes finish a batch sooner but split the memory budget between them. If SNAP runs out of memory, use fewer processes or a larger budget. ## Reprocessing and failures A product whose file already exists is not processed again, so a script can be rerun after an interruption. ::: {.callout-warning} ## Changed settings raise an error Filenames describe the burst and date, not the settings, so an existing file is checked before it is skipped. If it was made with other settings, such as another `grid_spacing` or DEM, the call raises a `ProductSettingsMismatchError` naming each difference, before SNAP runs. The same error is raised for a file that does not record every setting, such as one written by an earlier s1bursts version, and for a file made by another SNAP version, since the versions give slightly different results (see [SNAP versions](snap-processing.qmd#snap-versions)). A SNAP version that could not be determined is not compared. Write each experiment to its own `FileManager` directory, or pass `overwrite=True` to replace the old products. ::: A failed burst does not stop the batch by default. Each call returns an `OperationSummary`. Check it before using the results, since a time series with silent gaps can mislead. `raise_on_failure()` raises an error listing every failed burst, and otherwise returns the summary, so its products are one step away: ```{.python} rtc = processor.create_backscatter(bursts, crs="auto").raise_on_failure().products ``` Pass `continue_on_error=False` to `BurstProcessor` to stop at the first failure instead. [^small]: Small, D. (2011). Flattening gamma: Radiometric terrain correction for SAR imagery. *IEEE Transactions on Geoscience and Remote Sensing*, 49(8), 3081–3093. # Compositing A composite combines many burst products into one image per grid square. `ResolutionWeighter` writes composites on the 10 km squares of the Military Grid Reference System (MGRS), so results from different bursts, orbits and dates line up pixel for pixel. ## Get the grid squares `MGRSGrid` lists the grid squares that overlap an area. Its `resolution` is the pixel size of the composite and must equal the `grid_spacing` the products were processed with: ```{.python} from shapely.geometry import box from s1bursts import MGRSGrid aoi = box(7.03, 46.06, 7.12, 46.14) squares = list(MGRSGrid(aoi, resolution=20)) ``` Each square is an `MGRSGeoBox`, which knows its MGRS code, its pixel grid and its UTM zone, `square.crs`. Both classes come from [freezebase](https://lqgentner.github.io/freezebase/), whose guide [The MGRS grid](https://lqgentner.github.io/freezebase/user-guide/mgrs-grid.html) explains the grid and why it suits composites. The products must be in the zone of the square: process them with `crs=square.crs` and keep the default `align_to_standard_grid=True`. A product on a different grid is rejected. ::: {.callout-warning} ## Do not process products for a composite with `crs="auto"` `"auto"` takes each burst's own zone, and a burst whose center lies in the neighboring zone would be projected there, not in the zone of the square. ::: [From an area to composites](area-to-composites.qmd) shows the whole loop from an area to the composites of its squares. ## Compose with local resolution weighting `compose_rtc_lrw()` for backscatter and `compose_coherence_lrw()` for coherence combine all products that cover a square, whatever their orbit or direction. Each pixel is weighted toward the acquisitions with the highest spatial resolution at that pixel. This is called local resolution weighting (LRW). [Local resolution weighting](resolution-weighting.qmd) explains how it works. ```{.python} from s1bursts import ResolutionWeighter, TimeRange weighter = ResolutionWeighter(fm) window = TimeRange("2025-08-01", "2025-08-31") for square in squares: composite = weighter.compose_rtc_lrw( rtc.products, square, write_qa=True, time_range=window, ) ``` The method returns the composite, a record whose file `fm.get_product_path(composite)` holds the backscatter. Pass the `TimeRange` of your search as `time_range`, so that every square records the same compositing window. The weights come from the RTC products, so `compose_coherence_lrw()` takes the RTC products of the same dates as well: ```{.python} coherence_composite = weighter.compose_coherence_lrw( coherence.products, rtc.products, square, write_qa=True, time_range=window, ) ``` ## Quality layers With `write_qa=True`, two quality layers are written beside an LRW composite: - The number of observations (`QA-NUM`): how many acquisitions contributed to each pixel. - The composite quality map (`QA-CQM`): how fine the combined detail of each pixel is. The composite lists them in its `layers`, and `fm.open()` reads each layer's file: ```{.python} from s1bursts.product import Layer with fm.open(composite, Layer.QA_CQM) as src: cqm = src.read(1) ``` ## Composites per pass direction The pass composites keep ascending and descending orbits apart: | Method | Combines | Use it for | |---|---|---| | `compose_rtc_by_pass()` | RTC | Backscatter by viewing direction | | `compose_gtc_by_pass()` | GTC | Backscatter without terrain flattening | | `compose_coherence_by_pass()` | Coherence | Coherence by viewing direction | They return a dictionary of composites keyed by pass direction, `"Ascending"` and `"Descending"`. They take a plain average without resolution weights, and their only quality layer is the number of observations. ## Choose how to combine Every method uses `aggregation="mean"` by default. Pass `aggregation="median"`, or a number between 0 and 1 for a quantile, to reduce the influence of outliers. Composites are written in linear units, like the products they combine. Convert backscatter to decibels after reading it, as `10 * np.log10(gamma0)`, so that the averaging has already happened on linear values. All inputs of one call must share a polarization, the processing settings and the SNAP version, the inputs of both pass directions included, and coherence inputs must share a temporal baseline. A pass composite that already exists must match the settings and SNAP version of the direction being added. A SNAP version that could not be determined matches any. The RTC products that weight a coherence composite must share their settings among themselves, and the DEM and SNAP version with the coherence products. The composite records the settings of the coherence products only. Like processing, compositing skips a square whose files already exist, after checking that they were made on the same grid and with the same settings. Pass `overwrite=True` to replace composites made with other settings. # From an area to composites This page joins the steps of the previous pages into one script that takes an area of interest (AOI) to composites on the MGRS grid squares that cover it. s1bursts provides each step, and the script that joins them is yours to adapt. The script below assumes that all squares lie in one UTM zone, as they do for most areas. [The last section](#areas-across-a-utm-zone-boundary) adapts it to an area across a zone boundary. ## Outline Compositing is the only step that works one square at a time. Searching, downloading and processing serve all squares at once, so they run once for the whole area, and SNAP gets one batch of bursts that it can process in parallel: ```{.python} squares = list(MGRSGrid(aoi, resolution=20)) # Once for the whole area bursts = catalog.search(...) # every burst that reaches a square catalog.download(bursts) rtc = processor.create_backscatter(bursts, ...) # one batch, in the zone of the squares # Once per square for square in squares: square_rtc = ... # the products whose bursts reach it weighter.compose_rtc_lrw(square_rtc, square, ...) ``` The next three sections fill in these steps. Together, their code blocks form the complete script. The [Martigny tutorial](martigny-composite.qmd) runs the same script on real data. ## Find the bursts of every square `MGRSGrid` lists the squares that the AOI touches. Its `resolution` is the pixel size of the composites, so it is taken from the processing settings: ```{.python} import shapely from shapely.geometry import box from s1bursts import ( BurstCatalog, BurstProcessor, FileManager, MGRSGrid, ResolutionWeighter, RTCConfig, SearchQuery, TimeRange, ) fm = FileManager("data") catalog = BurstCatalog(fm) processor = BurstProcessor(fm) weighter = ResolutionWeighter(fm) aoi = box(7.03, 46.06, 7.12, 46.14) window = TimeRange("2025-08-01", "2025-08-31") config = RTCConfig(grid_spacing=20) squares = list(MGRSGrid(aoi, resolution=config.grid_spacing)) search_area = shapely.union_all([box(*square.geographic_extent.boundingbox) for square in squares]) bursts = catalog.search(SearchQuery(time_range=window, intersects=search_area, polarization="VV")) ``` The search covers the squares, not the AOI. `MGRSGrid` includes every square that the AOI touches, whole, and each composite covers its whole square. A search of the AOI would miss a burst that reaches only the part of an edge square outside the AOI, and that square would be composed from fewer acquisitions without a warning. A later run would then skip the incomplete composite as finished. `search_area` joins the bounding boxes of the squares in longitude and latitude, since the exact outline of a square has too many vertices for a catalog query. A bounding box is slightly larger than its square, so no burst is missed. ## Process in the zone of the squares A composite accepts only products in the UTM zone of its square, so the bursts are processed with that zone as `crs`, not with `crs="auto"`. [Compositing](compositing.qmd#get-the-grid-squares) explains why. One batch has one `crs`, so check first that all squares share a zone: ```{.python} if len({square.crs for square in squares}) > 1: raise ValueError("The squares span several UTM zones; process each zone on its own.") crs = squares[0].crs catalog.download(bursts).raise_on_failure() rtc = processor.create_backscatter(bursts, crs=crs, config=config).raise_on_failure().products ``` To set how many bursts SNAP processes at once, see [Run several jobs at once](processing.qmd#run-several-jobs-at-once). ## Compose each square Each square is composed from the products whose bursts cover it. `intersecting()` selects those bursts from the search result, and a product belongs to a burst footprint by its `burst_id` and `sub_swath`: ```{.python} from s1bursts.errors import CompositeError def reaching(products, bursts, square): """Return the products whose burst footprint reaches the square.""" footprints = {(b.burst_id, b.sub_swath) for b in bursts.intersecting(square)} return [p for p in products if (p.burst_id, p.sub_swath) in footprints] failed = {} for square in squares: square_rtc = reaching(rtc, bursts, square) if not square_rtc: continue # no burst reaches the square try: weighter.compose_rtc_lrw(square_rtc, square, write_qa=True, time_range=window) except CompositeError as error: failed[square.mgrs_code] = error ``` Passing the full product list to each square would give the same composites, since a product outside a square contributes nothing to it, but each product would be read once per square. Some squares cannot be composed, because of their own data. A square that no burst reached in the time window has no products, so the loop skips it. A burst footprint can also overlap a square without a valid pixel in it, and if that holds for all products of a square, the call raises an `EmptyCompositeError`. Such an error concerns one square and is a `CompositeError`, so the loop records it in `failed` and continues with the next square. An error that every square would meet is not a `CompositeError` and stops the loop: an `InconsistentInputsError` for products in another UTM zone or at another pixel size, or a `ProductSettingsMismatchError` for composites made with other settings. A composite contains one polarization. For VV and VH, run the script once with each. Rerunning the script is cheap. Downloads, products and composites that exist are skipped, so after a failure only the missing pieces are made. ## Add coherence Coherence composites follow the same outline, with pairs in place of bursts, and take the RTC products of the same dates as weights. A pair that starts at the end of the window ends after it, so the search reaches beyond the window by the baseline, and the bursts and pairs that start within the window are kept, as in [Pair acquisitions](finding-and-downloading.qmd#pair-acquisitions). These lines take the place of the search and of the download and processing calls: ```{.python} from datetime import timedelta from s1bursts import CoherenceConfig found = catalog.search(SearchQuery( time_range=TimeRange(window.start, window.end + timedelta(days=12)), intersects=search_area, polarization="VV", )) bursts = found.within(window) pairs = found.find_pairs(baseline=12).starting_within(window) catalog.download([*bursts, *pairs]).raise_on_failure() rtc = processor.create_backscatter(bursts, crs=crs, config=config).raise_on_failure().products coh = processor.create_coherences( pairs, crs=crs, config=CoherenceConfig(grid_spacing=config.grid_spacing), ).raise_on_failure().products ``` `download()` fetches both acquisitions of each pair and every burst only once. Each square's coherence is then composed with that square's RTC products: ```{.python} for square in squares: square_coh = reaching(coh, bursts, square) if not square_coh: continue # no pair reaches the square try: weighter.compose_coherence_lrw( square_coh, reaching(rtc, bursts, square), square, write_qa=True, time_range=window, ) except CompositeError as error: failed[square.mgrs_code] = error ``` ## Areas across a UTM zone boundary Each UTM zone is gridded on its own, so for an AOI that crosses a zone boundary, `MGRSGrid` returns squares of both zones, and they overlap along the boundary. `box(5.9, 46.0, 6.1, 46.1)`, for example, crosses 6° E and gives four squares in zone 31 and four in zone 32, among them 32TKR69, which covers half of 31TGL29. Composing all of them gives that ground one composite in each zone. The squares of each zone need their own processing batch, in their zone. Group the squares by zone, and run the script for each group, with `zone_squares` in place of `squares` and without the zone check: ```{.python} from collections import defaultdict squares_by_zone = defaultdict(list) for square in squares: squares_by_zone[square.crs].append(square) for crs, zone_squares in squares_by_zone.items(): ... # search, download, process and compose zone_squares as above ``` A burst near the boundary is found by the searches of both zones. It is downloaded once and processed once per zone. Another way, for any number of zones, is to run the whole script for each square on its own: search the bounding box of the square, keep the bursts `intersecting()` it, then download them, process them with `crs=square.crs` and compose. Each search is small, which suits an area whose zones contain more bursts than one search returns. Neighboring squares share most of their bursts, but each burst is still downloaded and processed only once per zone, since the later squares find its files and skip them. The cost is parallelism: each batch hands SNAP only the few bursts that are not processed yet. # Storing products on S3 `FileManager` can keep bursts, products and composites in an S3 bucket instead of on local disk. Give it an S3 path from [freezebase](https://lqgentner.github.io/freezebase/) instead of a local `Path`. Everything else in your workflow stays the same. ## Point the FileManager at a bucket `make_s3_upath()` builds the S3 path. Its `profile` names an entry in your AWS configuration files, which holds the access keys and, for S3 services other than AWS, the endpoint URL: ```{.python} from pathlib import Path from freezebase.s3 import make_s3_upath from s1bursts import FileManager fm = FileManager( make_s3_upath("s3://my-bucket/s1bursts", profile="my-profile"), slc_dir=Path("cache/slc"), ) ``` freezebase's [Accessing S3](https://lqgentner.github.io/freezebase/user-guide/accessing-s3.html) guide shows how to set up a profile and the other ways to give credentials. This example writes products and composites to the bucket but keeps the downloaded bursts on local disk, which saves storage costs for files you can download again. Omit `slc_dir` to store everything in the bucket. Read products with `fm.open()`, which passes the bucket's endpoint and credentials to GDAL. `rasterio.open()` on the path from `get_product_path()` does not. Your bucket is unrelated to CDSE's S3 storage. The `CDSE_S3_*` settings are only for [downloading bursts](finding-and-downloading.qmd#when-the-burst-service-fails). ::: {.callout-note} SNAP reads and writes only local files. s1bursts therefore copies each input to a temporary local folder, runs SNAP there, and uploads the result, so some local disk space is still needed. ::: # Using an external DEM Terrain correction needs a digital elevation model (DEM). By default SNAP downloads the Copernicus 30 m DEM by itself, so most users never need this page. Use an external DEM when a better local one exists, such as a national elevation model, or to avoid SNAP's [height error](snap-processing.qmd#corrections-s1bursts-applies) with the built-in Copernicus DEM. ## Prepare the DEM An external DEM must be prepared before SNAP reads it, because SNAP has two reported problems with external DEMs: - **It samples them half a pixel off**, to the south-east ([forum report](https://forum.step.esa.int/t/external-dem-is-sampled-half-a-pixel-off/46042)). SNAP fixed this in the 14.0.2 update of its DEM module. - **It knows only one height reference, EGM96**, which most modern DEMs do not use ([forum report](https://forum.step.esa.int/t/copernicus-dem-built-in-path-applies-egm96-while-the-data-is-egm2008/46041)). `prepare_dem_for_snap()` corrects both and writes a copy, leaving your original untouched: ```{.python} from pathlib import Path from s1bursts import prepare_dem_for_snap prepared = prepare_dem_for_snap( Path("swissalti3d_30m.tif"), source_crs="EPSG:2056+5728", download_missing_grids=True, ) ``` `source_crs` states both the map projection and the height reference of the DEM, joined by `+`. Here, `EPSG:2056` is the Swiss LV95 grid and `5728` the Swiss LN02 heights. Look these up in your DEM's documentation. Converting the heights can need correction grids from PROJ, which `download_missing_grids=True` fetches when they are freely available. The half-pixel shift is applied only for SNAP versions that still need it. `prepare_dem_for_snap()` checks the SNAP installation it finds the same way processing does, so the prepared copy is correct only for that installation and only for SNAP. Prepared copies are cached, so preparing the same file again costs nothing. ## Use the prepared DEM Pass the prepared file to the processing settings, together with a short name, which is required. The name, not the path, is written into every product's metadata, so you can tell later which DEM was used, and a rerun checks it before skipping an existing product: ```{.python} from s1bursts import RTCConfig config = RTCConfig( dem=prepared, dem_name="swissALTI3D 30 m", external_dem_apply_egm=False, ) ``` Keep `external_dem_apply_egm=False`: the preparation has already converted the heights, and SNAP would otherwise convert them a second time. ## Downsample a fine DEM first SNAP's terrain flattening and terrain correction (geocoding) take much longer on a fine DEM grid. So a DEM much finer than your products should be downsampled before preparing it, by averaging all fine pixels within each coarse one. freezebase's `vrt.warp()` averages by default. This call reduces the 2 m swissALTI3D to 30 m: ```{.python} from freezebase import vrt from freezebase.raster import write_cog write_cog(vrt.warp("swissalti3d_2m.tif", resolution=30), "swissalti3d_30m.tif") ``` ## Tutorials # Glacier backscatter time series This tutorial follows the Aletsch glacier in the Swiss Alps through 2022, one acquisition per month. It shows how backscatter changes with the seasons and how to turn a stack of products into a time series for one site. ```{python} #| include: false import os from pathlib import Path from dotenv import load_dotenv root = next(p for p in (Path.cwd(), *Path.cwd().parents) if (p / "pyproject.toml").exists()) load_dotenv(root / ".env") (root / "cache/docs").mkdir(parents=True, exist_ok=True) os.chdir(root / "cache/docs") ``` ## Process one acquisition per month The search covers a year of one burst footprint over the glacier. Keeping the first acquisition of each month leaves twelve, which are downloaded and processed into RTC backscatter at 20 m pixel size, in the burst's own UTM zone: ```{python} from s1bursts import BurstCatalog, BurstProcessor, FileManager, RTCConfig, SearchQuery, TimeRange from s1bursts.burst import BurstCollection fm = FileManager("data/glacier") catalog = BurstCatalog(fm) bursts = catalog.search(SearchQuery( time_range=TimeRange("2022-01-01", "2022-12-31"), burst_id=30346, sub_swath="IW1", polarization="VV", )) df = bursts.df.sort_values("date") monthly = BurstCollection(df.groupby(df.date.dt.month).head(1)) print(f"{len(bursts)} acquisitions in 2022, {len(monthly)} selected") ``` ```{python} #| output: false catalog.download(monthly).raise_on_failure() processor = BurstProcessor(fm) rtc = processor.create_backscatter( monthly, crs="auto", config=RTCConfig(grid_spacing=20), ).raise_on_failure() ``` ## Compare January and May The maps show the glacier in January and in May. The tributaries of the Great Aletsch Glacier meet at Konkordiaplatz, near the top center, and the glacier flows south-east, then south-west toward its tongue. In May the glacier and much of the terrain around it are far darker than in January. ```{python} #| fig-cap: "RTC backscatter around the Aletsch glacier in January and May 2022. The white rectangle is the sample area of the time series." import matplotlib.pyplot as plt import numpy as np from rasterio.warp import transform_bounds from rasterio.windows import from_bounds from shapely.geometry import box view = (7.95, 46.38, 8.15, 46.55) # lon/lat bounds of the map sample = box(8.065, 46.472, 8.072, 46.477) # on the glacier below Konkordiaplatz products = sorted(rtc.products, key=lambda p: p.date) fig, axes = plt.subplots(1, 2, figsize=(7.8, 3.7), sharey=True, layout="constrained") for ax, product in zip(axes, [products[0], products[4]]): with fm.open(product) as src: left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *view) gamma0 = src.read(1, window=from_bounds(left, bottom, right, top, src.transform)) x0, y0, x1, y1 = transform_bounds("EPSG:4326", src.crs, *sample.bounds) extent = [v / 1000 for v in (left, right, bottom, top)] image = ax.imshow(10 * np.log10(gamma0), extent=extent, cmap="viridis", vmin=-25, vmax=0) ax.add_patch(plt.Rectangle((x0 / 1000, y0 / 1000), (x1 - x0) / 1000, (y1 - y0) / 1000, fill=False, edgecolor="white")) ax.set(title=f"{product.date:%d %B %Y}", xlabel="Easting (km)") axes[0].set_ylabel("Northing (km)") fig.colorbar(image, ax=axes, shrink=0.8, label="Backscatter (dB)") plt.show() ``` ## Extract a time series The time series is the median backscatter within the sample rectangle on the glacier. The median is taken in linear units and converted to decibels afterwards, because averaging decibels gives a biased result. ```{python} import pandas as pd from rasterio.mask import mask from rasterio.warp import transform_geom from shapely.geometry import mapping records = [] for product in products: with fm.open(product) as src: region = transform_geom("EPSG:4326", src.crs, mapping(sample)) values, _ = mask(src, [region], indexes=1, crop=True, filled=False) valid = values.compressed() valid = valid[np.isfinite(valid)] records.append({"date": product.date, "gamma0_db": 10 * np.log10(np.median(valid))}) series = pd.DataFrame(records) ``` ```{python} #| fig-cap: "Median RTC backscatter in the sample rectangle, one acquisition per month in 2022." fig, ax = plt.subplots(figsize=(7.8, 3.3)) ax.plot(series.date, series.gamma0_db, "o-") ax.set(xlabel="Acquisition date", ylabel="Median backscatter (dB)") ax.grid(alpha=0.25) fig.autofmt_xdate() plt.show() ``` ## Interpret the curve Backscatter on the glacier stays near −10.5 dB from January to April, drops to −18.5 dB on 9 May, and recovers to about −9 dB in August and September. Wet snow is a likely cause of the drop: liquid water in the snowpack absorbs much of the radar signal, so wet snow appears dark (Rott, 1984)[^rott]. From October, values settle around −12 to −13 dB. The [next tutorial](glacier-coherence.qmd) computes coherence for the same burst. [^rott]: Rott, H. (1984). Synthetic aperture radar capabilities for snow and glacier monitoring. *Advances in Space Research*, 4(11), 241–246. ```{python} #| echo: false #| output: asis from datetime import date from importlib.metadata import version from s1bursts import SNAPConfig snap_version = ".".join(map(str, SNAPConfig().resolved_snap_version)) print( '::: {.callout-note collapse="true"}\n## Run details\n\n' f"Run on {date.today():%d %B %Y} with s1bursts {version('s1bursts')} and " f"ESA SNAP {snap_version}. Contains modified Copernicus Sentinel data (2022).\n:::" ) ``` # Glacier coherence through the seasons This tutorial computes coherence for one 12-day pair per month over the Aletsch glacier in 2022. Coherence shows where the surface stayed the same between two acquisitions, so it separates moving ice and changing snow from stable ground. ```{python} #| include: false import os from pathlib import Path from dotenv import load_dotenv root = next(p for p in (Path.cwd(), *Path.cwd().parents) if (p / "pyproject.toml").exists()) load_dotenv(root / ".env") (root / "cache/docs").mkdir(parents=True, exist_ok=True) os.chdir(root / "cache/docs") ``` ## Process one pair per month The search is the same as in the [backscatter tutorial](glacier-backscatter.qmd). `find_pairs()` forms every pair of acquisitions 12 days apart. Keeping the first pair of each month leaves twelve, which are processed into coherence: ```{python} from s1bursts import BurstCatalog, BurstProcessor, CoherenceConfig, FileManager, SearchQuery, TimeRange from s1bursts.burst import BurstPairCollection fm = FileManager("data/glacier") catalog = BurstCatalog(fm) bursts = catalog.search(SearchQuery( time_range=TimeRange("2022-01-01", "2022-12-31"), burst_id=30346, sub_swath="IW1", polarization="VV", )) df = bursts.find_pairs(baseline=12).df.sort_values("date_1") monthly = BurstPairCollection(df.groupby(df.date_1.dt.month).head(1)) print(f"{len(monthly)} pairs selected") ``` ```{python} #| output: false catalog.download(monthly).raise_on_failure() processor = BurstProcessor(fm) coherence = processor.create_coherences( monthly, crs="auto", config=CoherenceConfig(grid_spacing=20), ).raise_on_failure() ``` ## Coherence in winter and summer The maps compare a winter pair with a summer pair. In January almost the whole area is snow-covered and incoherent. Snow decorrelates quickly: new snow and melt change it between the two dates, and wet snow can lose its coherence within hours (Strozzi et al., 1999)[^strozzi]. In August the slopes around the glacier are highly coherent, while the glacier is dark, which makes its outline stand out. ```{python} #| fig-cap: "Coherence of a 12-day pair in winter and in summer 2022. The rectangles are the sample areas on the glacier (white) and off it (cyan)." import matplotlib.pyplot as plt from rasterio.warp import transform_bounds from rasterio.windows import from_bounds from shapely.geometry import box view = (7.95, 46.38, 8.15, 46.55) # lon/lat bounds of the map areas = { "glacier": (box(8.065, 46.472, 8.072, 46.477), "white"), "no glacier": (box(8.104, 46.448, 8.112, 46.453), "tab:cyan"), } products = sorted(coherence.products, key=lambda p: p.date_primary) fig, axes = plt.subplots(1, 2, figsize=(7.8, 3.7), sharey=True, layout="constrained") for ax, product in zip(axes, [products[0], products[7]]): with fm.open(product) as src: left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *view) coh = src.read(1, window=from_bounds(left, bottom, right, top, src.transform)) for area, color in areas.values(): x0, y0, x1, y1 = transform_bounds("EPSG:4326", src.crs, *area.bounds) ax.add_patch(plt.Rectangle((x0 / 1000, y0 / 1000), (x1 - x0) / 1000, (y1 - y0) / 1000, fill=False, edgecolor=color)) image = ax.imshow(coh, extent=[v / 1000 for v in (left, right, bottom, top)], cmap="magma", vmin=0, vmax=1) ax.set(title=f"{product.date_primary:%d %b} to {product.date_secondary:%d %b %Y}", xlabel="Easting (km)") axes[0].set_ylabel("Northing (km)") fig.colorbar(image, ax=axes, shrink=0.8, label="Coherence") plt.show() ``` ## Follow coherence through the year The time series takes the median coherence of each pair within the two rectangles, one on the glacier and one on a slope east of it without glacier. ```{python} #| fig-cap: "Median coherence of one 12-day pair per month in 2022, on the glacier and off it." import numpy as np import pandas as pd from rasterio.mask import mask from rasterio.warp import transform_geom from shapely.geometry import mapping records = [] for product in products: with fm.open(product) as src: for name, (area, _) in areas.items(): region = transform_geom("EPSG:4326", src.crs, mapping(area)) values, _ = mask(src, [region], indexes=1, crop=True, filled=False) valid = values.compressed() records.append({"date": product.date_primary, "area": name, "coherence": np.median(valid[np.isfinite(valid)])}) series = pd.DataFrame(records).pivot(index="date", columns="area", values="coherence") fig, ax = plt.subplots(figsize=(7.8, 3.3)) for name, (_, color) in areas.items(): ax.plot(series.index, series[name], "o-", label=name, color="tab:blue" if color == "white" else color) ax.set(xlabel="Date of the first acquisition", ylabel="Median coherence", ylim=(0, 1)) ax.legend() ax.grid(alpha=0.25) fig.autofmt_xdate() plt.show() ``` ## Interpret the result On the glacier, coherence stays near 0.12 all year. Values this low mean that the two acquisitions are essentially unrelated. They do not reach 0, because the coherence magnitude estimated from a finite window of pixels is biased toward higher values, most strongly where coherence is low (Touzi et al., 1999)[^touzi]. The glacier flows, and its surface melts or gains new snow, so 12 days are too long for it to stay coherent. Off the glacier, coherence follows the snow cover. It is low from November to May, rises in June, and stays between 0.57 and 0.74 until October, when the snow returns. Summer coherence therefore separates the glacier from the terrain around it, even where both look similar in backscatter. Shorter baselines keep more coherence. In 2022 only Sentinel-1A was operating, so 12 days is the shortest baseline available for these dates. [^strozzi]: Strozzi, T., Wegmüller, U., & Mätzler, C. (1999). Mapping wet snowcovers with SAR interferometry. *International Journal of Remote Sensing*, 20(12), 2395–2403. [^touzi]: Touzi, R., Lopes, A., Bruniquel, J., & Vachon, P. W. (1999). Coherence estimation for SAR imagery. *IEEE Transactions on Geoscience and Remote Sensing*, 37(1), 135–149. ```{python} #| echo: false #| output: asis from datetime import date from importlib.metadata import version from s1bursts import SNAPConfig snap_version = ".".join(map(str, SNAPConfig().resolved_snap_version)) print( '::: {.callout-note collapse="true"}\n## Run details\n\n' f"Run on {date.today():%d %B %Y} with s1bursts {version('s1bursts')} and " f"ESA SNAP {snap_version}. Contains modified Copernicus Sentinel data (2022).\n:::" ) ``` # Local resolution weighted composite of Martigny This tutorial runs the whole pipeline, from search to composites, for an area around Martigny in the Swiss canton of Valais. The town is surrounded by mountains, whose steep slopes are challenging terrain for SAR. The tutorial compares the single-pass composites, one of the ascending and one of the descending acquisitions, with a composite that combines both by local resolution weighting (LRW). ```{python} #| include: false import os from pathlib import Path from dotenv import load_dotenv root = next(p for p in (Path.cwd(), *Path.cwd().parents) if (p / "pyproject.toml").exists()) load_dotenv(root / ".env") (root / "cache/docs").mkdir(parents=True, exist_ok=True) os.chdir(root / "cache/docs") ``` ## Find the bursts in the area of interest The search finds the VV bursts that overlap the area's grid square in a 12-day window, 2025-08-13 to 2025-08-24. The area lies in one 10 km square of the MGRS grid, and the search covers the entire square, as in [From an area to composites](area-to-composites.qmd#find-the-bursts-of-every-square): ```{python} from shapely.geometry import box from s1bursts import ( BurstCatalog, BurstProcessor, FileManager, MGRSGrid, ResolutionWeighter, RTCConfig, SearchQuery, TimeRange, ) fm = FileManager("data/martigny") catalog = BurstCatalog(fm) aoi = box(7.01, 46.09, 7.055, 46.125) window = TimeRange("2025-08-13", "2025-08-24") (square,) = MGRSGrid(aoi, resolution=20) # the area lies in one square print(square.mgrs_code, square.crs) search_area = box(*square.geographic_extent.boundingbox) found = catalog.search(SearchQuery(time_range=window, intersects=search_area, polarization="VV")) bursts = found.intersecting(square) bursts.df.groupby(["pass_direction", "rel_orbit_number"]).agg( acquisitions=("date", lambda dates: dates.dt.date.nunique()), bursts=("burst_id", "size"), sub_swaths=("sub_swath", lambda names: ", ".join(sorted(names.unique()))), satellites=("parent_product_name", lambda names: ", ".join(sorted(names.str[:3].unique()))), ) ``` One ascending and two descending orbits cover the square, each imaged twice in these 12 days, with 14 bursts in total. ## Process backscatter The bursts are processed in the UTM zone of the square, at 20 m grid spacing: ```{python} #| output: false catalog.download(bursts).raise_on_failure() processor = BurstProcessor(fm) rtc = processor.create_backscatter( bursts, crs=square.crs, config=RTCConfig(grid_spacing=20), ).raise_on_failure().products ``` ## Compose the composites To compare local resolution weighting with single-pass composites, the square is composed both ways. `compose_rtc_lrw()` combines both pass directions and also writes the quality layers, and `compose_rtc_by_pass()` averages each direction on its own, without weights. For an area of several squares, [From an area to composites](area-to-composites.qmd#compose-each-square) loops over them. ```{python} weighter = ResolutionWeighter(fm) lrw = weighter.compose_rtc_lrw(rtc, square, write_qa=True, time_range=window) by_pass = weighter.compose_rtc_by_pass(rtc, square, time_range=window) ``` ## Compare the composites How finely a slope is resolved depends on how it faces the sensor. Sentinel-1 looks sideways, to the right of its track: east on ascending passes, west on descending ones. Slopes facing the sensor are compressed by foreshortening or layover and poorly resolved, while slopes facing away are well resolved or, if steep enough, fall into radar shadow, where that pass direction has no value ([Why mountains need it](resolution-weighting.qmd#why-mountains-need-it)). ```{python} #| fig-cap: "Backscatter composites of 2025-08-13 to 2025-08-24 in a 4 by 4 km window west of Martigny: the two ascending acquisitions (left), the four descending acquisitions (middle), and all six combined by local resolution weighting (right). White marks radar shadow, where a composite has no value." import matplotlib.pyplot as plt import numpy as np from s1bursts.product import Layer def read(composite, layer=Layer.MEASUREMENT): """Read one layer of a composite, with its extent in kilometers for imshow.""" with fm.open(composite, layer) as src: left, bottom, right, top = (v / 1000 for v in src.bounds) return src.read(1).astype("float32"), [left, right, bottom, top] view = (346, 350, 5106, 5110) # easting and northing in km, UTM zone 32N composites = { "Ascending": read(by_pass["Ascending"]), "Descending": read(by_pass["Descending"]), "LRW": read(lrw), } viridis = plt.get_cmap("viridis").copy() viridis.set_bad("white") # pixels without a value: radar shadow fig, axes = plt.subplots(1, 3, figsize=(7.6, 3.2), sharey=True, layout="constrained") for ax, (name, (data, extent)) in zip(axes, composites.items()): image = ax.imshow(10 * np.log10(data), extent=extent, cmap=viridis, vmin=-15, vmax=0) ax.set(title=name, xlabel="Easting (km)", xlim=view[:2], ylim=view[2:]) axes[0].set_ylabel("Northing (km)") fig.colorbar(image, ax=axes, shrink=0.5, label="Backscatter (dB)") plt.show() ``` Each single-pass composite has gaps where the steepest slopes facing away from the sensor fall into radar shadow, and compresses the slopes facing it into coarse streaks. The resolution-weighted composite averages all six acquisitions but gives the better-resolved ones more weight, so it has a value everywhere in this window, and its resolution varies less from slope to slope. ## Check quality layers The LRW composite comes with two quality layers, shown here for the same window: ```{python} #| fig-cap: "Quality layers of the resolution-weighted composite (left): the number of acquisitions per pixel (middle) and the local resolution (right), darkest where the ground is most coarsely resolved." # title: (layer, colorbar label, colormap, colorbar ticks, the first and last being the limits) layers = { "LRW composite": (Layer.MEASUREMENT, "Backscatter (dB)", "viridis", [-15, -10, -5, 0]), "Observations": (Layer.QA_NUM, "Number of observations", "cividis", [0, 2, 4, 6]), "Quality": (Layer.QA_CQM, "Local resolution (dB)", "cividis", [-12, -8, -4, 0, 4, 8]), } fig, axes = plt.subplots(1, 3, figsize=(7.6, 3.9), sharey=True, layout="constrained") for ax, (title, (layer, label, cmap, ticks)) in zip(axes, layers.items()): data, extent = read(lrw, layer) if layer != Layer.QA_NUM: data = 10 * np.log10(data) image = ax.imshow(data, extent=extent, cmap=cmap, vmin=ticks[0], vmax=ticks[-1]) fig.colorbar(image, ax=ax, orientation="horizontal", label=label, ticks=ticks) ax.set(title=title, xlabel="Easting (km)", xlim=view[:2], ylim=view[2:]) axes[0].set_ylabel("Northing (km)") plt.show() ``` The number of observations (`QA-NUM`) counts the acquisitions outside radar shadow that contribute to each pixel: all six almost everywhere, fewer where one pass direction is in shadow. The composite quality map (`QA-CQM`) is the local resolution of the composite itself. It is highest on slopes that face away from the sensor in some orbits, and lowest on the steepest slopes, which every orbit either compresses or leaves in shadow. A pixel with a low quality describes the ground more coarsely than its 20 m grid spacing suggests. [Local resolution weighting](resolution-weighting.qmd) explains both layers. ```{python} #| echo: false #| output: asis from datetime import date from importlib.metadata import version from s1bursts import SNAPConfig snap_version = ".".join(map(str, SNAPConfig().resolved_snap_version)) print( '::: {.callout-note collapse="true"}\n## Run details\n\n' f"Run on {date.today():%d %B %Y} with s1bursts {version('s1bursts')} and " f"ESA SNAP {snap_version}. Contains modified Copernicus Sentinel data (2025).\n:::" ) ``` ## Background # Products and bands Every product is a GeoTIFF, described in Python by a small metadata record. The record holds the product's identity, such as burst, date and type, and `FileManager.get_product_path()` turns it into the file's location. This page describes what the files contain. ## Burst products Each burst product is one GeoTIFF, with its bands named in the band descriptions: | Product | Bands | Values | |---|---|---| | RTC backscatter | `gamma0`, `scattering_area` | `gamma0` in linear power, `scattering_area` without units | | GTC backscatter | `beta0` | Linear power | | Coherence | `coherence` | Magnitude, 0 to 1, without units | Backscatter is stored as linear power, not decibels, in products and composites alike. Convert with `10 * log10(value)` only after averaging, since the mean of decibels is not the decibel of the mean. The `scattering_area` band describes how much ground each pixel covers. [Local resolution weighting](resolution-weighting.qmd) explains how composites use it. ## Composites A composite has one band, the quantity of its inputs: gamma-0, beta-0 or coherence, in the same units. Its band description is the polarization, such as `VV`. The bands of the quality layers are named `num_obs` and `cqm`. With `write_qa=True`, each quality layer of a composite is written to a file of its own: - The number of observations (`QA-NUM`) counts the acquisitions that contributed to each pixel. - The composite quality map (`QA-CQM`) rates how fine the combined detail of each pixel is, without units. Only LRW composites have it. [Local resolution weighting](resolution-weighting.qmd#the-quality-layers) explains both. ## Missing values Pixels without a value are `NaN`, for example outside the burst. Keep them missing in statistics and plots rather than replacing them with zero. The one exception is the number of observations, an integer layer that marks pixels without data as 0. Both are declared as the file's no-data value, so `src.read(1, masked=True)` masks them. ## Metadata tags Each GeoTIFF also has tags that say how to read it and how it was made. They include the units of the stored values, the acquisition dates, every processing setting, such as the DEM, the resampling or whether the topographic phase was removed, and the SNAP and s1bursts versions: ```{.python} with fm.open(product) as src: print(src.descriptions, src.units) print(src.tags()) ``` These tags make a file understandable on its own, even when separated from the code that wrote it. A composite records the processing settings of its inputs, which must all match. A resolution-weighted coherence composite records those of its coherence products, not of the RTC products that weight them. Where the [CEOS-ARD SAR specification](https://ceos.org/ard/) has a term, the tag names use it, so it can be looked up there. s1bursts borrows the vocabulary only and does not claim conformance with the specification. ## Filenames A filename contains only what distinguishes one product from another, in fixed fields separated by `_`. Each product type has its own pattern: | Product | Filename pattern | |---|---| | SLC burst | `S1_SLC_____.zip` | | RTC backscatter | `S1_RTC______.tif` | | GTC backscatter | `S1_GTC______.tif` | | Coherence | `S1_COH_______.tif` | | RTC LRW composite | `S1_RTC_LRW[_]____.tif` | | Coherence LRW composite | `S1_COH_LRW[_]____.tif` | | RTC or GTC pass composite | `S1_RTC_COMP_[_QA-NUM]____.tif`, with `GTC` for GTC | | Coherence pass composite | `S1_COH_COMP_[_QA-NUM]____.tif` | The part in square brackets appears only in the name of a quality layer. The fields are: | Field | Meaning | Example | |---|---|---| | `` | Relative orbit number, three digits | `015` | | `` | Pass direction: `A` ascending, `D` descending | `A` | | `` | ESA burst ID, six digits | `030346` | | `` | Sub-swath | `IW1` | | `` | Polarization | `VV` | | `` | Acquisition date, `YYYYMMDD` | `20220801` | | ``, `` | Dates of the primary and secondary acquisitions of a pair | `20220801`, `20220813` | | `` | Days between the acquisitions of a pair, two digits | `12` | | `` | UTM zone of the product's CRS | `32N` | | `` | MGRS grid square | `32TLS40` | | ``, `` | First and last day of the `time_range`, or of the acquisitions without one | `20250801`, `20250831` | | `` | Quality layer: `QA-NUM` or `QA-CQM` | `QA-CQM` | The RTC product of the [Quickstart](quickstart.qmd), burst 30346 in sub-swath IW1 on ascending orbit 15, is named `S1_RTC_015A_030346_IW1_VV_20220801_32N.tif`, and the number of observations of an LRW composite is named `S1_RTC_LRW_QA-NUM_32TLS40_VV_20250801_20250831.tif`. `FileManager` keeps SLC bursts in `s1-slc`, burst products in `s1-prod` and composites in `s1-comp`, below the directory it is given. Processing settings are not part of the name. They are recorded in the tags instead, and a rerun compares them before it skips an existing file. See [Reprocessing and failures](processing.qmd#reprocessing-and-failures). # Local resolution weighting Local resolution weighting (LRW) combines acquisitions pixel by pixel and gives more weight to the acquisitions that resolve the terrain more finely at that pixel. s1bursts follows Small et al. (2022)[^small] and applies the method to backscatter and coherence alike. ## Why mountains need it A radar looks sideways, so the ground resolution of a mountain slope depends on how the slope faces the sensor. A slope facing the sensor is compressed into few pixels and coarsely resolved. A slope facing away is stretched over many pixels and finely resolved, unless it is so steep that it lies in radar shadow. Ascending and descending orbits look from opposite sides, so each resolves slopes that the other resolves poorly. ## The weights The `scattering_area` band of an RTC product contains the normalized ground area $A_i$ that acquisition $i$ observes in each pixel. Its inverse is the local resolution, $r_i = 1/A_i$. At each pixel, the weights are the local resolutions, normalized to sum to one over the $M$ acquisitions with a valid observation there: $$ W_i = \frac{r_i}{\sum_{j=1}^{M} r_j} $$ The composite of a measured value $x$ is the weighted mean $$ x_c = \sum_{i=1}^{M} W_i \, x_i , $$ or the weighted median or quantile with `aggregation="median"` or a quantile level. The weights depend only on the viewing geometry, so $x$ can be backscatter $\gamma^0$ or coherence. A coherence product has no scattering area of its own, so each pair is weighted by the RTC product of its primary acquisition, which is why `compose_coherence_lrw()` takes RTC products as well. An acquisition in radar shadow has no valid area there and does not count toward $M$. Bursts of the same orbit and date are mosaicked first, so the overlap of neighboring bursts counts once. The edge of every burst is dropped, because multilooking or the coherence window darkens it. `burst_edge_trim` turns this off, and its width follows the looks or window and the grid spacing. Each pixel is taken from the burst whose data edge is farther away, which splits every overlap at its midline, as SNAP's TOPSAR-Deburst and TOPSAR-Merge do. A pixel's scattering area comes from the same burst as its value, or for coherence from the RTC product of that burst. ## Capping the local resolution Close to radar shadow, $A_i$ approaches zero and $r_i$ grows without bound, although the simulated area is least reliable there. `clip_max_resolution` sets a cap $r_\text{max}$ and uses $\min(r_i, r_\text{max})$ in place of $r_i$ in the weights, so that one acquisition cannot dominate such pixels. No cap is applied by default. For scale, flat terrain has a local resolution of about 0.6 to 1.0. ## The quality layers The number of observations (`QA-NUM`) is $M$, the number of acquisitions that contributed to each pixel. The composite quality map (`QA-CQM`) is the local resolution of the composite: the inverse of the mean scattering area, weighted as the composite is: $$ \text{CQM} = \left( \sum_{i=1}^{M} W_i \, A_i \right)^{-1} $$ A higher value means a finer resolution. The quality map describes the viewing geometry only, not the accuracy of the measured value. With a median or quantile, it is still the weighted mean from [The weights](#the-weights), so it describes the weighting rather than the acquisitions that supplied each value. [^small]: Small, D., Rohner, C., Miranda, N., Ruetschi, M., & Schaepman, M. E. (2022). Wide-area analysis-ready radar backscatter composites. *IEEE Transactions on Geoscience and Remote Sensing*, 60, 5201814. # SNAP processing s1bursts relies on ESA SNAP for the SAR and InSAR processing. This page describes the steps SNAP applies to a burst, how to inspect them for your own settings, which SNAP problems s1bursts corrects, and how the supported SNAP versions differ. ## The processing steps Each product runs through a fixed chain of SNAP operators, called a processing graph. s1bursts builds the graph from the product's settings object, so the settings decide which steps run and with which parameters. The RTC backscatter graph has eight steps with the default settings: | Step | What it does | Settings | |---|---|---| | `Apply-Orbit-File` | Replaces the predicted satellite orbit with the precise one | | | `ThermalNoiseRemoval` | Removes the instrument's own noise | | | `Calibration` | Converts the raw values to beta-0 radar brightness | | | `TOPSAR-Deburst` | Keeps only the burst's valid lines and samples, dropping the empty margin around it | | | `Multilook` | Averages neighboring pixels to reduce noise and runtime | `multilook` (`None` skips it) | | `Terrain-Flattening` | Corrects brightness for the slope of the terrain (RTC only) | `dem_oversampling`, `additional_overlap`, DEM settings | | `ReplaceZerosWithNoData` | Marks the empty border of the burst as missing ([why](#corrections-s1bursts-applies)) | `replace_zeros_with_nan` (`False` skips it) | | `Terrain-Correction` | Places every pixel on the map, using the DEM | `grid_spacing`, `img_resampling`, `align_to_standard_grid`, DEM settings | The DEM settings are `dem`, `dem_resampling` and, for an [external DEM](external-dems.qmd), `external_dem_nodata` and `external_dem_apply_egm`. The GTC graph of `GTCConfig` is the same without `Terrain-Flattening`. Coherence compares two acquisitions, so its graph differs. Both acquisitions get precise orbits, then `Back-Geocoding` aligns the second onto the first pixel by pixel, using the DEM. `Coherence` estimates the similarity of the two within a moving window, set by `window`. `TOPSAR-Deburst`, `ReplaceZerosWithNoData` and `Terrain-Correction` then run as in the backscatter graph. The coherence graph has no noise removal, calibration or multilooking. ## Inspect a graph Every settings object can show the graph it produces, as a diagram or as the XML file that SNAP runs, to check what your settings change. ### As a diagram `to_mermaid()` returns the graph as text in the [Mermaid](https://mermaid.js.org/) diagram language. In a notebook, the [mermaid-py](https://github.com/ouhammmourachid/mermaid-py) package draws it. Install it separately, since s1bursts does not depend on it: ```{.python} from mermaid import Mermaid from s1bursts import RTCConfig config = RTCConfig(grid_spacing=20) Mermaid(config.to_mermaid()) ``` mermaid-py draws the diagram through the online mermaid.ink service. Without internet access, paste the text into a `mermaid` code block on GitHub or in Quarto instead. For `RTCConfig(grid_spacing=20)`, the graph looks like this: ```{python} #| echo: false #| output: asis from s1bursts import RTCConfig config = RTCConfig(grid_spacing=20) print(f"```{{mermaid}}\n{config.to_mermaid().replace('flowchart LR', 'flowchart TD')}\n```") ``` ### As XML `to_xml()` returns the graph in the XML format of SNAP's Graph Processing Tool (GPT), with every parameter SNAP receives: ```{python} xml = config.to_xml() ``` ```{python} #| echo: false #| output: asis print(f':::{{.callout-note collapse="true"}}\n## The full graph\n\n```xml\n{xml}\n```\n:::') ``` The graph leaves two placeholders, `${source}` for the burst file and `${crs}` for the UTM zone. s1bursts fills them in for every burst it processes, which is why the UTM zone is passed to the `create_*` methods rather than stored in the settings. `CoherenceConfig` returns its graph the same way, with `${primary}` and `${secondary}` in place of `${source}`. ### For your SNAP version The coherence graph differs slightly between SNAP 13 and 14. Without a version, s1bursts builds the graph for SNAP 13. Pass the version of your installation to see the exact graph that processing runs: ```{.python} from s1bursts import CoherenceConfig, SNAPConfig version = SNAPConfig().resolved_snap_version xml = CoherenceConfig().to_xml(snap_version=version) ``` ## Corrections s1bursts applies s1bursts works around four SNAP problems, so you do not need to: - **External DEMs are read half a pixel off** ([forum report](https://forum.step.esa.int/t/external-dem-is-sampled-half-a-pixel-off/46042)). `prepare_dem_for_snap()` shifts the DEM to compensate, but only before the 14.0.2 update of SNAP's DEM module, which fixed it. - **DEM heights use the wrong sea-level model** ([forum report](https://forum.step.esa.int/t/copernicus-dem-built-in-path-applies-egm96-while-the-data-is-egm2008/46041)). SNAP knows only the EGM96 model, even for DEMs based on another. `prepare_dem_for_snap()` converts external DEM heights itself (see [Using an external DEM](external-dems.qmd)). The built-in Copernicus DEM is affected too. To avoid it, download the Copernicus DEM and prepare it as an external DEM with `source_crs="EPSG:4326+3855"`. - **Burst edges come out too dark** ([forum report](https://forum.step.esa.int/t/sentinel-1-slc-terrain-correction-operator-uses-zero-values-in-interpolation-instead-of-treating-them-as-no-data/45336)). SNAP mixes the empty border of the input into the edge pixels during terrain correction, which leaves dark seams along the burst edges. s1bursts marks that border as missing before terrain correction. - **Areas without DEM coverage become zeros.** Where the DEM has no data, for example over the sea at the edge of a burst, SNAP can write zeros instead of missing values. s1bursts turns these into `NaN`. ## SNAP versions s1bursts works with SNAP 13 and 14. It reads the installed version and adapts the graphs to it, so the same code runs on either. The version is recorded in every product's `PROCESSING_SOFTWARE_VERSION` tag. An existing product of another version is not skipped, and products of two versions are not combined into one composite. ::: {.callout-note} ## Keep one SNAP version per study SNAP versions give slightly different results for the same burst, so a time series that switches version contains an artificial step. Between SNAP 13 and 14, backscatter differs by a median of about 0.05 dB and coherence by about 0.004 on average. ::: SNAP 14 also adds two coherence options, exposed as `CoherenceConfig.correct_coherence_bias` and `remove_local_phase_ramp`. Both are off by default and need SNAP 14 or newer. Products record them in the `COHERENCE_BIAS_CORRECTED` and `LOCAL_PHASE_RAMP_REMOVED` tags: `false` under SNAP 13, and `N/A` when the SNAP version could not be determined, since SNAP 14 then applies its own default, which corrects the bias. ## Validation s1bursts products are checked against products of the same bursts from independent processors: backscatter against [OPERA RTC-S1](https://hyp3-docs.asf.alaska.edu/guides/opera_rtc_product_guide/) from NASA/JPL, and coherence against [ASF HyP3](https://hyp3-docs.asf.alaska.edu/guides/burst_insar_product_guide) burst InSAR products. Both use the same Copernicus DEM, so the comparison tests the processing itself. The comparisons are tests in the repository, [`test_opera_validation.py`](https://github.com/lqgentner/s1bursts/blob/main/tests/test_opera_validation.py) and [`test_insar_validation.py`](https://github.com/lqgentner/s1bursts/blob/main/tests/test_insar_validation.py). [`REFERENCE.md`](https://github.com/lqgentner/s1bursts/blob/main/tests/fixtures/REFERENCE.md) describes the reference products and what each test checks. ## Development # Contributing Contributions are welcome as issues and pull requests on [GitHub](https://github.com/lqgentner/s1bursts). This page covers setting up a development environment and building this documentation. [`AGENTS.md`](https://github.com/lqgentner/s1bursts/blob/main/AGENTS.md) holds the full development guidelines and the reasoning behind the design. ## Set up and test The project uses [uv](https://docs.astral.sh/uv/). These commands install the environment and run the checks that every pull request must pass: ```{.bash} uv sync uv run pytest uv run ruff format uv run ruff check --fix uv run pyrefly check ``` The default tests run offline. Tests that download data or run SNAP are marked `integration` and run with `uv run pytest -m integration` if CDSE credentials and a SNAP installation are available. ## Check the prose [Vale](https://vale.sh/docs/install) checks the Markdown and Quarto files. It is installed separately, not through uv. `.vale.ini` sets the rules and gives the reason for each rule that is off. Download the rule packages once, then check the tracked files: ```{.bash} vale sync vale $(git ls-files '*.md' '*.qmd') ``` When the spelling check reports a correct term, add it to `.vale/styles/config/vocabularies/s1bursts/accept.txt`. ## Build the documentation The site is built with [Great Docs](https://posit-dev.github.io/great-docs/), which also needs [Quarto](https://quarto.org/docs/get-started/). Guide pages live in `user_guide/`, the API reference comes from the docstrings, and `great-docs.yml` configures the site: ```{.bash} uv sync --group docs uv run great-docs build uv run great-docs preview ``` The build writes the site to `great-docs/_site/`, which is generated and never edited by hand. ## Frozen tutorials The tutorials download real data and run SNAP, which takes too long for every build. Their results are therefore frozen: `_freeze/` stores each page's rendered output and figures, and builds reuse it without running the code. CI builds the site from `_freeze/` alone, without credentials or SNAP. Refresh a frozen page after changing its code, its text or the API it uses. The frozen output includes the text, so the site shows an edit only after the page is refreshed. This reruns the page, so it needs `.env` and SNAP. Downloads and products are cached in `cache/docs/` and reused. Refresh one page per command, then rebuild: ```{.bash} uv run great-docs freeze user_guide/10-glacier-backscatter.qmd uv run great-docs build ``` Review the new figures and the `_freeze/` diff before committing. Each frozen page ends with its run details: the date, the s1bursts version and the SNAP version it was run with.