Visualizing multi-dimensional OPERA-DISP with titiler-multidim

Authors: Henry Rodman(Development Seed), Harshini Girish(UAH), Alex Mandel(Development Seed)

Date: September 25, 2025

Description: TiTiler-MultiDim is a TiTiler application designed to serve multidimensional rasters—like OPERA-DISP NetCDF—directly as web map tiles. It lets you pick a variable (e.g., water_mask or displacement layers) and slice along dimensions (time, burst, polarization, etc.), then renders those selections on-the-fly into XYZ/TileJSON that you can drop into Folium/Leaflet without downloading the file. The notebook centers on this workflow: open a multidim asset, choose variable + dims + styling (colormap, range), and stream dynamic tiles for quick, interactive visualization.

Run This Notebook

To access and run this tutorial within MAAP’s Algorithm Development Environment (ADE), please refer to the “Getting started with the MAAP” section of our documentation.

Disclaimer: it is highly recommended to run a tutorial within MAAP’s ADE, which already includes packages specific to MAAP, such as maap-py. Running the tutorial outside of the MAAP ADE may lead to errors.

Additional Resources

About the Dataset

The Level-3 OPERA Sentinel-1 Surface Displacement (DISP) product is generated through interferometric time-series analysis of Level-2 Coregistered Sentinel-1 Single Look Complex (CSLC) datasets. Using a hybrid Persistent Scatterer (PS) and Distributed Scatterer (DS) approach, this product quantifies Earth’s surface displacement in the radar line-of-sight. The DISP products enable the detection of anthropogenic and natural surface changes, including subsidence, tectonic deformation, and landslides.

The OPERA DISP suite comprises complementary datasets derived from Sentinel-1 and NISAR inputs, designated as DISP-S1 and DISP-NI, respectively. Each product, created per acquisition, adheres to a consistent structure, HDF5 file format, file-naming convention, and a 30 m spatial posting. This collection specifically includes DISP-S1 products, derived from Sentinel-1 data. For visualization and quick exploration, the Pangeo Image can be used for these datasets.

Source: OPERA Surface Displacement from Sentinel-1

Install/Import Packages

Make sure the following libraries are installed before running the notebook

[1]:

import os
import earthaccess
import httpx
from folium import GeoJson, LayerControl, Map, TileLayer
from pprint import pprint
from shapely.geometry import box, mapping
from pathlib import Path

Searching the Data

This performs a granule search using the earthaccess.granule_query() function on the OPERA Sentinel-1 displacement product collection.

[2]:

auth = earthaccess.login()
granule_query = (
    earthaccess.granule_query()
    .short_name("OPERA_L3_DISP-S1_V1")
    .bounding_box(-121, 38, -120, 39)
)

granules = granule_query.get(1)
granule = granules[0]

print(granule)

Collection: {'ShortName': 'OPERA_L3_DISP-S1_V1', 'Version': '1'}
Spatial coverage: {'HorizontalSpatialDomain': {'Geometry': {'GPolygons': [{'Boundary': {'Points': [{'Latitude': 39.16866, 'Longitude': -121.23333}, {'Latitude': 38.66165, 'Longitude': -124.06639}, {'Latitude': 37.8195, 'Longitude': -123.86328}, {'Latitude': 37.95784, 'Longitude': -122.92612}, {'Latitude': 37.35143, 'Longitude': -122.78615}, {'Latitude': 37.66388, 'Longitude': -120.92917}, {'Latitude': 39.16866, 'Longitude': -121.23333}]}}]}}}
Temporal coverage: {'RangeDateTime': {'BeginningDateTime': '2016-07-05T02:07:28Z', 'EndingDateTime': '2017-01-07T02:06:47Z'}}
Size(MB): 376.585862159729
Data: ['https://datapool.asf.alaska.edu/DISP/OPERA-S1/OPERA_L3_DISP-S1_IW_F09157_VV_20160705T020728Z_20170107T020647Z_v1.0_20250408T163918Z.nc', 'https://datapool.asf.alaska.edu/DISP/OPERA-S1/OPERA_L3_DISP-S1_IW_F09157_VV_20160705T020728Z_20170107T020647Z_v1.0_20250408T163918Z.zarr.json.gz', 'https://datapool.asf.alaska.edu/DISP/OPERA-S1/OPERA_L3_DISP-S1_IW_F09157_VV_short_wavelength_displacement.zarr.json.gz']

Downloading OPERA-DISP granules

Creates the local folder to fetch the selected OPERA-DISP files.Then the progress indicators for the download tasks.

[ ]:

!mkdir -p /projects/my-public-bucket/opera-disp

[ ]:

download_dir = "/projects/my-public-bucket/opera-disp/"
downloaded = earthaccess.download(
    granules,
    download_dir,
)

Building nc_path and the S3 nc_key

Finds the first downloaded NetCDF by extension then constructs a shared-workspace S3 URI so the file can be referenced.

[ ]:

nc_path = next(
    (
        str(Path(fn).relative_to(download_dir))
        for fn in downloaded
        if Path(fn).suffix == ".nc"
    ),
    None,
)

assert nc_path is not None, f"No NetCDF (.nc) file downloaded to {download_dir}"

nc_key = (
    f"s3://maap-ops-workspace/shared/{os.getenv('CHE_WORKSPACE_NAMESPACE')}"
    f"/opera-disp/{nc_path}"
)

Requesting TileJSON from TiTiler-MultiDim

Send a request to /WebMercatorQuad/tilejson.json with url=nc_key, the target variable, and styling. TiTiler reads/slices the NetCDF on the fly and returns TileJSON with bounds, center, zoom limits, and a tiles URL template. Then use the returned tiles template to create a Leaflet tileLayer and render the layer in the map without downloading the file.

[6]:

TITILER_MULTIDIM_ENDPOINT = "https://staging.openveda.cloud/api/titiler-multidim"

params = {
    "url": nc_key,
    "rescale": [[-0.1, 0.05]],
    "colormap_name": "rdbu",
    "variable": "displacement",
    "minzoom": 7,
}

tilejson = httpx.get(
    f"{TITILER_MULTIDIM_ENDPOINT}/WebMercatorQuad/tilejson.json",
    params=params,
    timeout=None,
).json()

pprint(tilejson)

{'bounds': [-124.14140929951725,
            37.08503078714323,
            -120.82130131653709,
            39.24070040130846],
 'center': [-122.48135530802716, 38.16286559422585, 7],
 'maxzoom': 12,
 'minzoom': 7,
 'scheme': 'xyz',
 'tilejson': '2.2.0',
 'tiles': ['https://staging.openveda.cloud/api/titiler-multidim/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?url=s3%3A%2F%2Fmaap-ops-workspace%2Fshared%2Fharshinigirish%2Fopera-disp%2FOPERA_L3_DISP-S1_IW_F09157_VV_20160705T020728Z_20170107T020647Z_v1.0_20250408T163918Z.nc&rescale=%5B-0.1%2C+0.05%5D&colormap_name=rdbu&variable=displacement'],
 'version': '1.0.0'}

Visualising tile layer and granule bounds in Leaflet

This builds a Leaflet TileLayer from the TileJSON (tiles[0], minzoom) so the OPERA-DISP raster streams as XYZ tiles. It also constructs a GeoJSON polygon from tilejson["bounds"] and styles it (orange outline) to show the granule footprint. Finally, it centers a map on the bounds midpoint, adds both layers, and enables a LayerControl for toggling.

[12]:

# create an XYZ tile layer for the leaflet map
tiles = TileLayer(
    name="OPERA-DISP",
    tiles=tilejson["tiles"][0],
    min_zoom=tilejson["minzoom"],
    opacity=1,
    attr="NASA",
    overlay=True,
    control=True,
)

# load the granule bounding box so we can plot it on the map
geojson = {"type": "Feature", "geometry": mapping(box(*tilejson["bounds"])), "properties": {}}

geojson_layer = GeoJson(
    data=geojson,
    style_function=lambda x: {
        "opacity": 1,
        "dashArray": "1",
        "fillOpacity": 0,
        "weight": 2,
        "color": "orange",
    },
    name="granule bounds",
    overlay=True,
    control=True,
    show=True,
)

m = Map(
    tiles="OpenStreetMap",
    location=(
        (tilejson["bounds"][1] + tilejson["bounds"][3]) / 2,
        (tilejson["bounds"][0] + tilejson["bounds"][2]) / 2
    ),
    zoom_start=8,
)
tiles.add_to(m)
geojson_layer.add_to(m)

LayerControl(collapsed=False).add_to(m)

m

[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Note: The plot may take 20–30 seconds to load.