Visualize NASA OPERA layers¶
This notebook renders OPERA granules on a GeoLibre map through titiler-cmr. titiler-cmr reads the Cloud-Optimized GeoTIFFs and handles Earthdata authentication on the server, so the notebook needs no NASA login. This is the same path the GeoLibre OPERA plugin's Display button uses.
It covers four products over Valencia, Spain, around the October 2024 DANA flood:
- DSWx-HLS and DSWx-S1 water classification, with categorical colormaps
- RTC-S1 SAR backscatter, converted to decibels with a band-math expression
- DIST-ALERT vegetation disturbance
The default endpoint
https://titiler-cmr.opengeos.orgis a shared demo service. For heavy use, pointTITILER_CMRat a titiler-cmr deployment you control.
import json
import re
import geopandas as gpd
import requests
from geolibre import Map
from shapely.geometry import Polygon, box
Helpers¶
search_granules is the CMR client from notebook 01. opera_tiles turns one granule and band into an XYZ tile URL.
CMR_GRANULES_URL = "https://cmr.earthdata.nasa.gov/search/granules.umm_json"
BAND_PATTERN = re.compile(r"_(B\d{2}_[A-Za-z0-9-]+|[VH]{2}|[A-Z][A-Za-z0-9-]+)\.tiff?$")
def granule_geometry(umm):
"""Build a footprint polygon from a UMM-G granule record.
Args:
umm: The ``umm`` block of a CMR UMM-G granule item.
Returns:
A shapely Polygon, or None when the record has no usable geometry.
"""
geometry = umm.get("SpatialExtent", {}).get("HorizontalSpatialDomain", {})
geometry = geometry.get("Geometry", {})
rects = geometry.get("BoundingRectangles") or []
polygons = geometry.get("GPolygons") or []
if rects:
r = rects[0]
return box(
r["WestBoundingCoordinate"],
r["SouthBoundingCoordinate"],
r["EastBoundingCoordinate"],
r["NorthBoundingCoordinate"],
)
if polygons:
points = polygons[0]["Boundary"]["Points"]
return Polygon([(p["Longitude"], p["Latitude"]) for p in points])
return None
def search_granules(short_name, bbox, start_date, end_date, count=100):
"""Search NASA CMR for OPERA granules (no Earthdata login required).
Args:
short_name: OPERA collection short name, e.g. ``"OPERA_L3_DSWX-HLS_V1"``.
bbox: Search area as ``[west, south, east, north]`` in degrees.
start_date: First day of the search window, ``"YYYY-MM-DD"``.
end_date: Last day of the search window, ``"YYYY-MM-DD"``.
count: Maximum number of granules to return.
Returns:
A GeoDataFrame with one footprint per granule, sorted by date.
"""
params = {
"short_name": short_name,
"bounding_box": ",".join(str(v) for v in bbox),
"temporal": f"{start_date}T00:00:00Z,{end_date}T23:59:59Z",
"page_size": count,
}
response = requests.get(CMR_GRANULES_URL, params=params, timeout=60)
response.raise_for_status()
rows = []
for item in response.json().get("items", []):
umm = item["umm"]
dates = umm.get("TemporalExtent", {}).get("RangeDateTime", {})
# CMR lists each file twice (HTTPS and S3); keep the HTTPS copy.
urls = {}
for link in umm.get("RelatedUrls", []):
url = link.get("URL", "")
match = BAND_PATTERN.search(url)
if url.startswith("https://") and match:
urls.setdefault(match.group(1), url)
rows.append(
{
"granule_ur": umm["GranuleUR"],
"concept_id": item["meta"]["collection-concept-id"],
"begin_date": dates.get("BeginningDateTime", ""),
"bands": ",".join(urls),
"urls": urls,
"geometry": granule_geometry(umm),
}
)
columns = ["granule_ur", "concept_id", "begin_date", "bands", "urls", "geometry"]
gdf = gpd.GeoDataFrame(rows, columns=columns, geometry="geometry", crs="EPSG:4326")
return gdf.sort_values("begin_date").reset_index(drop=True)
TITILER_CMR = "https://titiler-cmr.opengeos.org"
def titiler_params(granule, band, assets_regex, **render):
"""Build the titiler-cmr query that pins one granule and one band.
Args:
granule: A row of the GeoDataFrame returned by ``search_granules``.
band: Band (asset) name, e.g. ``"B01_WTR"`` or ``"VV"``.
assets_regex: Regex titiler-cmr uses to find band names in file names.
**render: Extra titiler parameters such as ``rescale``,
``colormap_name``, ``colormap`` (a dict), or ``expression``.
Returns:
A dict of query parameters.
"""
params = {
"collection_concept_id": granule["concept_id"],
"granule_ur": granule["granule_ur"],
"assets": band,
"assets_regex": assets_regex,
}
for key, value in render.items():
if value is None:
continue
if key == "colormap":
value = json.dumps(value)
params[key] = value
if "expression" in params:
# Expressions refer to the requested asset as b1.
params["asset_as_band"] = "true"
return params
def opera_tiles(granule, band, assets_regex, **render):
"""Ask titiler-cmr for an XYZ tile URL that renders one OPERA granule.
titiler-cmr reads the Cloud-Optimized GeoTIFF and handles Earthdata
authentication server-side, so no NASA login is needed here.
Args:
granule: A row of the GeoDataFrame returned by ``search_granules``.
band: Band (asset) name to render.
assets_regex: Regex titiler-cmr uses to find band names in file names.
**render: Extra titiler parameters, see ``titiler_params``.
Returns:
A ``(tile_url, bounds)`` tuple for ``Map.add_tile_layer``.
"""
url = f"{TITILER_CMR}/rasterio/WebMercatorQuad/tilejson.json"
params = titiler_params(granule, band, assets_regex, **render)
response = requests.get(url, params=params, timeout=120)
response.raise_for_status()
tilejson = response.json()
return tilejson["tiles"][0], tilejson.get("bounds")
Create the map¶
DSWx-HLS water classification¶
B01_WTR is categorical, and titiler-cmr does not apply the file's embedded color table, so
pass an explicit colormap that maps each class to RGBA. "Not water" and fill are transparent
so the water sits on top of the basemap. This is the colormap the plugin uses.
HLS is optical, so pick a scene with little cloud. The first clear look after the flood is on 2 November 2024.
DSWX_HLS_COLORMAP = {
"0": [0, 0, 0, 0], # not water
"1": [0, 0, 255, 255], # open water
"2": [135, 206, 250, 255], # partial surface water
"252": [0, 255, 255, 255], # snow/ice
"253": [128, 128, 128, 255], # cloud/cloud shadow
"254": [0, 0, 128, 255], # ocean masked
"255": [0, 0, 0, 0], # fill
}
dswx_hls = search_granules("OPERA_L3_DSWX-HLS_V1", bbox, "2024-11-02", "2024-11-02")
dswx_hls[["granule_ur", "begin_date"]]
for _, granule in dswx_hls.iterrows():
url, bounds = opera_tiles(
granule, "B01_WTR", r"B[0-9]{2}_[A-Z]+", colormap=DSWX_HLS_COLORMAP
)
m.add_tile_layer(
url,
name=f"DSWx-HLS {granule['begin_date'][:10]} {granule['granule_ur'].split('_')[3]}",
bounds=bounds,
attribution="NASA OPERA",
)
m.add_legend(
"DSWx-HLS water",
legend_dict={
"Open water": "#0000ff",
"Partial surface water": "#87cefa",
"Snow/ice": "#00ffff",
"Cloud/cloud shadow": "#808080",
"Ocean masked": "#000080",
},
)
DSWx-S1 water classification¶
DSWx-S1 comes from Sentinel-1 radar, so it sees through cloud. Its classes differ from DSWx-HLS: class 3 is inundated vegetation, and 250/251 mask steep terrain and radar layover/shadow. Here is the descending pass on 31 October 2024, two days after the flood.
DSWX_S1_COLORMAP = {
"0": [0, 0, 0, 0], # not water
"1": [0, 0, 255, 255], # open water
"3": [0, 168, 132, 255], # inundated vegetation
"250": [0, 0, 0, 0], # HAND masked
"251": [0, 0, 0, 0], # layover/shadow masked
"255": [0, 0, 0, 0], # fill
}
dswx_s1 = search_granules("OPERA_L3_DSWX-S1_V1", bbox, "2024-10-31", "2024-10-31")
dswx_s1_granule = dswx_s1[dswx_s1["granule_ur"].str.contains("_T30SYJ_")].iloc[0]
url, bounds = opera_tiles(
dswx_s1_granule, "B01_WTR", r"B[0-9]{2}_[A-Z]+", colormap=DSWX_S1_COLORMAP
)
dswx_s1_layer = m.add_tile_layer(
url, name="DSWx-S1 2024-10-31", bounds=bounds, attribution="NASA OPERA"
)
m.add_legend(
"DSWx-S1 water",
legend_dict={"Open water": "#0000ff", "Inundated vegetation": "#00a884"},
position="bottom-right",
)
RTC-S1 backscatter in decibels¶
RTC-S1 stores backscatter as linear power. A rio-tiler expression converts it to decibels on
the fly (b1 is the requested band), and rescale stretches -30 to 0 dB to grey. Smooth open
water reflects the radar away from the sensor, so flooded fields show up dark.
rtc = search_granules("OPERA_L2_RTC-S1_V1", bbox, "2024-10-31", "2024-10-31")
rtc[["granule_ur", "begin_date"]]
rtc_layers = []
for _, granule in rtc.iterrows():
url, bounds = opera_tiles(
granule,
"VV",
r"(VV|VH|HH|HV)",
expression="10*log10(b1)",
rescale="-30,0",
colormap_name="gray",
)
rtc_layers.append(
m.add_tile_layer(
url,
name=f"RTC-S1 VV dB {granule['granule_ur'].split('_')[3]}",
bounds=bounds,
attribution="NASA OPERA",
)
)
DIST-ALERT vegetation disturbance¶
VEG-DIST-STATUS tracks vegetation loss, from a first detection through confirmed and finished
disturbance. DIST file names have no B01_ prefix, so the band regex matches the VEG- and
GEN- layer names instead.
DIST_STATUS_COLORMAP = {
"0": [0, 0, 0, 0], # no disturbance
"1": [255, 255, 178, 255], # first detection, <50% loss
"2": [254, 204, 92, 255], # provisional, <50% loss
"3": [253, 141, 60, 255], # confirmed, <50% loss
"4": [252, 187, 161, 255], # first detection, >=50% loss
"5": [251, 106, 74, 255], # provisional, >=50% loss
"6": [203, 24, 29, 255], # confirmed, >=50% loss
"7": [189, 189, 189, 255], # finished, <50% loss
"8": [99, 99, 99, 255], # finished, >=50% loss
"255": [0, 0, 0, 0], # no data
}
dist = search_granules("OPERA_L3_DIST-ALERT-HLS_V1", bbox, "2024-11-05", "2024-11-05")
granule = dist[dist["granule_ur"].str.contains("_T30SYJ_")].iloc[0]
url, bounds = opera_tiles(
granule, "VEG-DIST-STATUS", r"(VEG|GEN)-[A-Z-]+", colormap=DIST_STATUS_COLORMAP
)
dist_layer = m.add_tile_layer(
url, name="DIST-ALERT 2024-11-05", bounds=bounds, attribution="NASA OPERA"
)
Continuous bands¶
Continuous bands such as DSWx-HLS B10_DEM render flat without a stretch. Pass rescale and a
named colormap_name (any titiler colormap, such as terrain, viridis, or magma).
granule = dswx_hls.iloc[0]
url, bounds = opera_tiles(
granule, "B10_DEM", r"B[0-9]{2}_[A-Z]+", rescale="0,1500", colormap_name="terrain"
)
dem_layer = m.add_tile_layer(url, name="DSWx-HLS DEM", bounds=bounds)
m.set_layer_opacity(dem_layer, 0.7)
m.add_colorbar(
colormap="terrain",
vmin=0,
vmax=1500,
label="Elevation",
units="m",
position="top-right",
)
Compare radar and water¶
Hide the layers drawn on top, then drag the slider to compare the RTC-S1 backscatter (left) with the DSWx-S1 water map derived from it (right). You can also toggle layers in the Layers panel.
m.hide_layer(dem_layer)
m.hide_layer(dist_layer)
m.split_map(left_layers=rtc_layers, right_layers=dswx_s1_layer)
Query a pixel¶
titiler-cmr's /point endpoint returns the band value at a location. This is the value under
L'Albufera lagoon in the DSWx-S1 scene (1 = open water).
lon, lat = -0.35, 39.33
params = titiler_params(dswx_s1_granule, "B01_WTR", r"B[0-9]{2}_[A-Z]+")
response = requests.get(
f"{TITILER_CMR}/rasterio/point/{lon},{lat}", params=params, timeout=60
)
response.raise_for_status()
response.json()["assets"][0]["values"]
Save the map¶
Save the map as a GeoLibre project. Open it later in the GeoLibre web or desktop app with the
OPERA plugin, or reload it here with m.load_project(...).