How to use the API with the Rasterio Backend¶

Since most of the API documentation uses the xarray backend, this notebook provides an example of how to use most endpoints with the rasterio backend.

See also the API docs.

Setup¶

In [1]:

import json
import os

import earthaccess
import httpx

from folium import Map, TileLayer

titiler_endpoint = os.getenv(
    "TITILER_CMR_ENDPOINT", "https://openveda.cloud/api/titiler-cmr"
)

import json import os import earthaccess import httpx from folium import Map, TileLayer titiler_endpoint = os.getenv( "TITILER_CMR_ENDPOINT", "https://openveda.cloud/api/titiler-cmr" )

Identify the dataset¶

The Harmonized Landsat Sentinel-2 dataset is available in two collections in CMR. This example will use data from the HLSL30.002 (Landsat) dataset. You can find the HLSL30.002 dataset using the earthaccess.search_datasets function.

In [2]:

datasets = earthaccess.search_datasets(keyword="HLSL30")
ds = datasets[0]

collection_concept_id = ds["meta"]["concept-id"]
print("Collection Concept-Id: ", collection_concept_id)
print("Abstract: ", ds["umm"]["Abstract"])

datasets = earthaccess.search_datasets(keyword="HLSL30") ds = datasets[0] collection_concept_id = ds["meta"]["concept-id"] print("Collection Concept-Id: ", collection_concept_id) print("Abstract: ", ds["umm"]["Abstract"])

Collection Concept-Id:  C2021957657-LPCLOUD
Abstract:  The Harmonized Landsat Sentinel-2 (HLS) project provides consistent surface reflectance (SR) and top of atmosphere (TOA) brightness data from a virtual constellation of satellite sensors. The Operational Land Imager (OLI) is housed aboard the joint NASA/USGS Landsat 8 and Landsat 9 satellites, while the Multi-Spectral Instrument (MSI) is mounted aboard Europe’s Copernicus Sentinel-2A, Sentinel-2B, and Sentinel-2C satellites. The combined measurement enables global observations of the land every 1.6 days at 30-meter (m) spatial resolution. The HLS project uses a set of algorithms to obtain seamless products from OLI and MSI that include atmospheric correction, cloud and cloud-shadow masking, spatial co-registration and common gridding, illumination and view angle normalization, and spectral bandpass adjustment.

The HLSL30 product provides 30-m Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (NBAR) and is derived from Landsat 8/9 OLI data products. The [HLSS30](https://doi.org/10.5067/HLS/HLSS30.002) and HLSL30 products are gridded to the same resolution and Military Grid Reference System ([MGRS](https://hls.gsfc.nasa.gov/products-description/tiling-system/)) tiling system and thus are “stackable” for time series analysis.

The HLSL30 product is provided in Cloud Optimized GeoTIFF (COG) format, and each band is distributed as a separate file. There are 11 bands included in the HLSL30 product along with one quality assessment (QA) band and four angle bands. See the User Guide for a more detailed description of the individual bands provided in the HLSL30 product.

The HLS project is funded by NASA’s Satellite Needs Working Group (SNWG) which provides data products developed to meet the needs of stakeholders from US government agencies.

For information on [known issues](https://lpdaac.usgs.gov/documents/2434/HLS_v2.0_L30_known_issues_April2026.pdf) see the anomalies link under the Documents and Resources tab.

Explore the collection using the /compatibility endpoint¶

In [3]:

compatibility_response = httpx.get(
    f"{titiler_endpoint}/compatibility",
    params={"collection_concept_id": collection_concept_id},
    timeout=None,
).json()

print(json.dumps(compatibility_response, indent=2))

compatibility_response = httpx.get( f"{titiler_endpoint}/compatibility", params={"collection_concept_id": collection_concept_id}, timeout=None, ).json() print(json.dumps(compatibility_response, indent=2))

{
  "concept_id": "C2021957657-LPCLOUD",
  "backend": "rasterio",
  "datetime": [
    {
      "RangeDateTimes": [
        {
          "BeginningDateTime": "2013-04-11T00:00:00.000Z"
        }
      ],
      "EndsAtPresentFlag": true,
      "PrecisionOfSeconds": 0,
      "TemporalResolution": {
        "Value": 1,
        "Unit": "Day"
      }
    }
  ],
  "variables": null,
  "dimensions": null,
  "coordinates": null,
  "example_assets": {
    "0": "s3://lp-prod-protected/HLSL30.020/HLS.L30.T59WPT.2013101T001445.v2.0/HLS.L30.T59WPT.2013101T001445.v2.0.B09.tif"
  },
  "sample_asset_raster_info": {
    "bounds": [
      600000.0,
      7690200.0,
      709800.0,
      7800000.0
    ],
    "crs": "PROJCS[\"UTM Zone 59, Northern Hemisphere\",GEOGCS[\"Unknown datum based upon the WGS 84 ellipsoid\",DATUM[\"Not specified (based on WGS 84 spheroid)\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",171],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]",
    "band_metadata": [
      [
        "b1",
        {}
      ]
    ],
    "band_descriptions": [
      [
        "b1",
        "Cirrus"
      ]
    ],
    "dtype": "int16",
    "nodata_type": "Nodata",
    "colorinterp": [
      "gray"
    ],
    "scales": [
      0.0001
    ],
    "offsets": [
      0.0
    ],
    "colormap": null,
    "driver": "GTiff",
    "count": 1,
    "width": 3660,
    "height": 3660,
    "overviews": [
      2,
      4,
      8,
      16
    ],
    "nodata_value": -9999.0
  },
  "links": [
    {
      "rel": "tilejson",
      "href": "https://openveda.cloud/api/titiler-cmr/rasterio/WebMercatorQuad/tilejson.json?collection_concept_id=C2021957657-LPCLOUD&temporal={temporal}",
      "title": "TileJSON",
      "type": "application/json"
    },
    {
      "rel": "map",
      "href": "https://openveda.cloud/api/titiler-cmr/rasterio/WebMercatorQuad/map.html?collection_concept_id=C2021957657-LPCLOUD&temporal={temporal}",
      "title": "Map viewer",
      "type": "text/html"
    },
    {
      "rel": "tile",
      "href": "https://openveda.cloud/api/titiler-cmr/rasterio/tiles/WebMercatorQuad/{z}/{x}/{y}?collection_concept_id=C2021957657-LPCLOUD&temporal={temporal}",
      "title": "Map tile",
      "type": "image/png"
    }
  ]
}

The details from the sample granule show that it has 18 assets (you would need to look more into what each of the assets represents). To properly configure the assets for titiler-cmr we will need to use the bands_regex parameter to identify the bands that we want to be available for visualizations. The datetime key shows the reported temporal range from CMR which indicates that the dataset has granules from 2015-11-28 to present.

Display tiles in an interactive map¶

The /tilejson.json endpoint will provide a parameterized xyz tile URL that can be added to an interactive map.

In [4]:

r = httpx.get(
    f"{titiler_endpoint}/rasterio/WebMercatorQuad/tilejson.json",
    params=(
        ("collection_concept_id", collection_concept_id),
        # Temporal range in form of `start_date/end_date`
        ("temporal", "2024-10-01T00:00:00Z/2024-10-31T23:59:59Z"),
        # Sort by cloud cover in ascending order to render less cloudy images!
        ("sort_key", "cloud_cover"),
        # We know that the HLS collection dataset is stored as File per Band
        # so we need to pass a `asset_regex` option to assign `assets` to each URL
        ("assets_regex", "B[0-9][0-9]"),
        # True Color Image B04,B03,B02
        ("assets", "B04"),
        ("assets", "B03"),
        ("assets", "B02"),
        # The data is in type of Uint16 so we need to apply some
        # rescaling/color_formula in order to create PNGs
        ("color_formula", "Gamma RGB 3.5 Saturation 1.7 Sigmoidal RGB 15 0.35"),
        # We need to set min/max zoom because we don't want to use lowerzoom level (e.g 0)
        # which will results in useless large scale query
        ("minzoom", 8),
        ("maxzoom", 13),
    ),
    timeout=None,
).json()

print(r)

r = httpx.get( f"{titiler_endpoint}/rasterio/WebMercatorQuad/tilejson.json", params=( ("collection_concept_id", collection_concept_id), # Temporal range in form of `start_date/end_date` ("temporal", "2024-10-01T00:00:00Z/2024-10-31T23:59:59Z"), # Sort by cloud cover in ascending order to render less cloudy images! ("sort_key", "cloud_cover"), # We know that the HLS collection dataset is stored as File per Band # so we need to pass a `asset_regex` option to assign `assets` to each URL ("assets_regex", "B[0-9][0-9]"), # True Color Image B04,B03,B02 ("assets", "B04"), ("assets", "B03"), ("assets", "B02"), # The data is in type of Uint16 so we need to apply some # rescaling/color_formula in order to create PNGs ("color_formula", "Gamma RGB 3.5 Saturation 1.7 Sigmoidal RGB 15 0.35"), # We need to set min/max zoom because we don't want to use lowerzoom level (e.g 0) # which will results in useless large scale query ("minzoom", 8), ("maxzoom", 13), ), timeout=None, ).json() print(r)

{'tilejson': '3.0.0', 'version': '1.0.0', 'scheme': 'xyz', 'tiles': ['https://openveda.cloud/api/titiler-cmr/rasterio/tiles/WebMercatorQuad/{z}/{x}/{y}?collection_concept_id=C2021957657-LPCLOUD&temporal=2024-10-01T00%3A00%3A00Z%2F2024-10-31T23%3A59%3A59Z&sort_key=cloud_cover&assets_regex=B%5B0-9%5D%5B0-9%5D&assets=B04&assets=B03&assets=B02&color_formula=Gamma+RGB+3.5+Saturation+1.7+Sigmoidal+RGB+15+0.35&tilesize=512'], 'minzoom': 8, 'maxzoom': 13, 'bounds': [-180.0, -90.0, 180.0, 90.0], 'center': [0.0, 0.0, 8]}

In [5]:

bounds = r["bounds"]
m = Map(location=(47.590266824611675, -91.03729840730689), zoom_start=r["maxzoom"] - 2)

TileLayer(
    tiles=r["tiles"][0],
    opacity=1,
    attr="NASA",
).add_to(m)
m

bounds = r["bounds"] m = Map(location=(47.590266824611675, -91.03729840730689), zoom_start=r["maxzoom"] - 2) TileLayer( tiles=r["tiles"][0], opacity=1, attr="NASA", ).add_to(m) m

Out[5]:

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

Render NDVI using the expression parameter¶

The expression parameter can be used to render images from an expression of a combination of the individual bands.

In [6]:

r = httpx.get(
    f"{titiler_endpoint}/rasterio/WebMercatorQuad/tilejson.json",
    params=(
        ("collection_concept_id", collection_concept_id),
        # Temporal range in form of `start_date/end_date`
        ("temporal", "2024-06-01T00:00:00Z/2024-06-30T23:59:59Z"),
        ("sort_key", "cloud_cover"),
        # We know that the HLS collection dataset is stored as File per Band
        # so we need to pass a `assets_regex` option to assign `assets` to each URL
        ("assets_regex", "B[0-9][0-9]"),
        ("asset_as_band", True),
        # NDVI
        ("expression", "(B05-B04)/(B05+B04)"),
        # Need red (B04) and nir (B05) for NDVI
        ("assets", "B05"),
        ("assets", "B04"),
        # The data is in type of Uint16 so we need to apply some
        # rescaling/color_formula in order to create PNGs
        ("colormap_name", "greens"),
        ("rescale", "-1,1"),
        # We need to set min/max zoom because we don't want to use lowerzoom level (e.g 0)
        # which will results in useless large scale query
        ("minzoom", 8),
        ("maxzoom", 13),
    ),
    timeout=None,
).json()

m = Map(location=(47.9221313337365, -91.65432884883238), zoom_start=r["maxzoom"] - 1)


TileLayer(
    tiles=r["tiles"][0],
    opacity=1,
    attr="NASA",
).add_to(m)

m

r = httpx.get( f"{titiler_endpoint}/rasterio/WebMercatorQuad/tilejson.json", params=( ("collection_concept_id", collection_concept_id), # Temporal range in form of `start_date/end_date` ("temporal", "2024-06-01T00:00:00Z/2024-06-30T23:59:59Z"), ("sort_key", "cloud_cover"), # We know that the HLS collection dataset is stored as File per Band # so we need to pass a `assets_regex` option to assign `assets` to each URL ("assets_regex", "B[0-9][0-9]"), ("asset_as_band", True), # NDVI ("expression", "(B05-B04)/(B05+B04)"), # Need red (B04) and nir (B05) for NDVI ("assets", "B05"), ("assets", "B04"), # The data is in type of Uint16 so we need to apply some # rescaling/color_formula in order to create PNGs ("colormap_name", "greens"), ("rescale", "-1,1"), # We need to set min/max zoom because we don't want to use lowerzoom level (e.g 0) # which will results in useless large scale query ("minzoom", 8), ("maxzoom", 13), ), timeout=None, ).json() m = Map(location=(47.9221313337365, -91.65432884883238), zoom_start=r["maxzoom"] - 1) TileLayer( tiles=r["tiles"][0], opacity=1, attr="NASA", ).add_to(m) m

Out[6]:

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

Visualize a specific granule¶

TiTiler-CMR allows you to visualize a specific granule with the granule_ur query parameter.

In [7]:

r = httpx.get(
    f"{titiler_endpoint}/rasterio/WebMercatorQuad/tilejson.json",
    params=(
        ("collection_concept_id", "C2021957657-LPCLOUD"),
        # Unique granule id (granule_ur)
        ("granule_ur", "HLS.L30.T44RQR.2026061T045429.v2.0"),
        # We know that the HLS collection dataset is stored as File per Band
        # so we need to pass a `asset_regex` option to assign `assets` to each URL
        ("assets_regex", "B[0-9][0-9]"),
        # True Color Image B04,B03,B02
        ("assets", "B04"),
        ("assets", "B03"),
        ("assets", "B02"),
        # The data is in type of Uint16 so we need to apply some
        # rescaling/color_formula in order to create PNGs
        ("color_formula", "Gamma RGB 3.5 Saturation 1.7 Sigmoidal RGB 15 0.35"),
    ),
    timeout=None,
).json()

print(r)

m = Map(location=(r["center"][1], r["center"][0]), zoom_start=9)

TileLayer(
    tiles=r["tiles"][0],
    opacity=1,
    attr="NASA",
).add_to(m)

m

r = httpx.get( f"{titiler_endpoint}/rasterio/WebMercatorQuad/tilejson.json", params=( ("collection_concept_id", "C2021957657-LPCLOUD"), # Unique granule id (granule_ur) ("granule_ur", "HLS.L30.T44RQR.2026061T045429.v2.0"), # We know that the HLS collection dataset is stored as File per Band # so we need to pass a `asset_regex` option to assign `assets` to each URL ("assets_regex", "B[0-9][0-9]"), # True Color Image B04,B03,B02 ("assets", "B04"), ("assets", "B03"), ("assets", "B02"), # The data is in type of Uint16 so we need to apply some # rescaling/color_formula in order to create PNGs ("color_formula", "Gamma RGB 3.5 Saturation 1.7 Sigmoidal RGB 15 0.35"), ), timeout=None, ).json() print(r) m = Map(location=(r["center"][1], r["center"][0]), zoom_start=9) TileLayer( tiles=r["tiles"][0], opacity=1, attr="NASA", ).add_to(m) m

{'tilejson': '3.0.0', 'version': '1.0.0', 'scheme': 'xyz', 'tiles': ['https://openveda.cloud/api/titiler-cmr/rasterio/tiles/WebMercatorQuad/{z}/{x}/{y}?collection_concept_id=C2021957657-LPCLOUD&granule_ur=HLS.L30.T44RQR.2026061T045429.v2.0&assets_regex=B%5B0-9%5D%5B0-9%5D&assets=B04&assets=B03&assets=B02&color_formula=Gamma+RGB+3.5+Saturation+1.7+Sigmoidal+RGB+15+0.35&tilesize=512'], 'minzoom': 0, 'maxzoom': 18, 'bounds': [83.01552348, 26.99956652, 84.14924457, 28.0103867], 'center': [83.582384025, 27.50497661, 0]}

Out[7]:

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

GeoJSON Statistics¶

The /statistics endpoint can be used to get summary statistics for a geojson Feature or FeatureCollection.

In [8]:

geojson = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {},
            "geometry": {
                "coordinates": [
                    [
                        [-91.65432884883238, 47.9221313337365],
                        [-91.65432884883238, 47.86503396133904],
                        [-91.53842043960762, 47.86503396133904],
                        [-91.53842043960762, 47.9221313337365],
                        [-91.65432884883238, 47.9221313337365],
                    ]
                ],
                "type": "Polygon",
            },
        }
    ],
}

geojson = { "type": "FeatureCollection", "features": [ { "type": "Feature", "properties": {}, "geometry": { "coordinates": [ [ [-91.65432884883238, 47.9221313337365], [-91.65432884883238, 47.86503396133904], [-91.53842043960762, 47.86503396133904], [-91.53842043960762, 47.9221313337365], [-91.65432884883238, 47.9221313337365], ] ], "type": "Polygon", }, } ], }

In [9]:

r = httpx.post(
    f"{titiler_endpoint}/rasterio/statistics",
    params=(
        ("collection_concept_id", collection_concept_id),
        # Temporal range in form of `start_date/end_date`
        ("temporal", "2024-07-01T00:00:00Z/2024-07-10T23:59:59Z"),
        # We know that the HLS collection dataset is stored as File per Band
        # so we need to pass a `band_regex` option to assign `bands` to each URL
        ("assets_regex", "B[0-9][0-9]"),
        # NDVI
        ("expression", "(b1-b2)/(b1+b2)"),
        # Need red (B04) and nir (B05) for NDVI
        ("assets", "B05"),
        ("assets", "B04"),
    ),
    json=geojson,
    timeout=None,
).json()

print(json.dumps(r, indent=2))

r = httpx.post( f"{titiler_endpoint}/rasterio/statistics", params=( ("collection_concept_id", collection_concept_id), # Temporal range in form of `start_date/end_date` ("temporal", "2024-07-01T00:00:00Z/2024-07-10T23:59:59Z"), # We know that the HLS collection dataset is stored as File per Band # so we need to pass a `band_regex` option to assign `bands` to each URL ("assets_regex", "B[0-9][0-9]"), # NDVI ("expression", "(b1-b2)/(b1+b2)"), # Need red (B04) and nir (B05) for NDVI ("assets", "B05"), ("assets", "B04"), ), json=geojson, timeout=None, ).json() print(json.dumps(r, indent=2))

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -91.65432884883238,
              47.9221313337365
            ],
            [
              -91.65432884883238,
              47.86503396133904
            ],
            [
              -91.53842043960762,
              47.86503396133904
            ],
            [
              -91.53842043960762,
              47.9221313337365
            ],
            [
              -91.65432884883238,
              47.9221313337365
            ]
          ]
        ]
      },
      "properties": {
        "statistics": {
          "b1": {
            "min": -75.4,
            "max": 26.6,
            "mean": 0.5241462608197205,
            "count": 57155.80078125,
            "sum": 29957.999263649046,
            "std": 0.6059345886116388,
            "median": 0.6050372066399542,
            "majority": 0.75,
            "minority": -75.4,
            "unique": 47491.0,
            "histogram": [
              [
                1,
                0,
                2,
                1,
                0,
                0,
                16,
                57615,
                12,
                2
              ],
              [
                -75.4,
                -65.2,
                -55.00000000000001,
                -44.80000000000001,
                -34.60000000000001,
                -24.400000000000006,
                -14.20000000000001,
                -4.000000000000014,
                6.199999999999989,
                16.39999999999999,
                26.6
              ]
            ],
            "valid_percent": 99.74,
            "masked_pixels": 149.0,
            "valid_pixels": 57649.0,
            "description": "(B05_NIR-B04_Red)/(B05_NIR+B04_Red)",
            "percentile_2": 0.04376706768027161,
            "percentile_98": 0.8688275512788047
          }
        },
        "used_assets": [
          "HLS.L30.T15UXP.2024186T165705.v2.0",
          "HLS.L30.T15UWP.2024186T165705.v2.0"
        ]
      }
    }
  ]
}

In [ ]: