rasterio backend example: HLS¶

This notebook demonstrates how to use titiler-cmr's rasterio backend to render tiles for collections with granules that contain two-dimensional raster data.

In addition to showing how to parameterize tilesson endpoint requests this guide also demonstrates how to use the earthaccess package to search for collections and how the relevant details for tiling can be extracted from the search results.

Requirements¶

To run some of the chunks in this notebook you will need to install a few packages:

• earthaccess
• folium
• httpx

!pip install folium httpx earthaccess

In [1]:

import earthaccess
import httpx
import json

from folium import Map, TileLayer

import earthaccess import httpx import json from folium import Map, TileLayer

In [2]:

# titiler_endpoint = "http://localhost:8081"  # docker network endpoint
titiler_endpoint = "https://staging.openveda.cloud/api/titiler-cmr"  # deployed endpoint

# titiler_endpoint = "http://localhost:8081" # docker network endpoint titiler_endpoint = "https://staging.openveda.cloud/api/titiler-cmr" # deployed endpoint

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 [3]:

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

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

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

Concept-Id:  C2021957295-LPCLOUD
Abstract:  The Harmonized Landsat Sentinel-2 (HLS) project provides consistent surface reflectance data from the Operational Land Imager (OLI) aboard the joint NASA/USGS Landsat 8 satellite and the Multi-Spectral Instrument (MSI) aboard Europe’s Copernicus Sentinel-2A, Sentinel-2B, and Sentinel-2C satellites. The combined measurement enables global observations of the land every 2–3 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 HLSS30 product provides 30-m Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (NBAR) and is derived from Sentinel-2A, Sentinel-2B, and Sentinel-2C MSI data products. The HLSS30 and [HLSL30](https://doi.org/10.5067/HLS/HLSL30.002) 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 HLSS30 product is provided in Cloud Optimized GeoTIFF (COG) format, and each band is distributed as a separate COG. There are 13 bands included in the HLSS30 product along with four angle bands and a quality assessment (QA) band. See the User Guide for a more detailed description of the individual bands provided in the HLSS30 product.

Known Issues

* Unrealistically high aerosol and low surface reflectance over bright areas: The atmospheric correction over bright targets occasionally retrieves unrealistically high aerosol and thus makes the surface reflectance too low. High aerosol retrievals, both false high aerosol and realistically high aerosol, are masked when quality bits 6 and 7 are both set to 1 (see Table 9 in the [User Guide](https://lpdaac.usgs.gov/documents/1698/HLS_User_Guide_V2.pdf)); the corresponding spectral data should be discarded from analysis.

* Issues over high latitudes: For scenes greater than or equal to 80 degrees north, multiple overpasses can be gridded into a single MGRS tile resulting in an L30 granule with data sensed at two different times. In this same area, it is also possible that Landsat overpasses that should be gridded into a single MGRS tile are actually written as separate data files. Finally, for scenes with a latitude greater than or equal to 65 degrees north, ascending Landsat scenes may have a slightly higher error in the BRDF correction because the algorithm is calibrated using descending scenes.

* Fmask omission errors: There are known issues regarding the Fmask band of this data product that impacts HLSL30 data prior to April of 2022. The HLS Fmask data band may have omission errors in water detection for cases where water detection using spectral data alone is difficult, and omission and commission errors in cloud shadow detection for areas with great topographic relief. This issue does not impact other bands in the dataset.

* Inconsistent snow surface reflectance between Landsat and Sentinel-2: The HLS snow surface reflectance can be highly inconsistent between Landsat and Sentinel-2. When assessed on same-day acquisitions from Landsat and Sentinel-2, Landsat reflectance is generally higher than Sentinel-2 reflectance in the visible bands.

* Unrealistically high snow surface reflectance in the visible bands: By design, the Land Surface Reflectance Code (LaSRC) atmospheric correction does not attempt aerosol retrieval over snow; instead, a default aerosol optical thickness (AOT) is used to drive the snow surface reflectance. If the snow detection fails, the full LaSRC is used in both AOT retrieval and surface reflectance derivation over snow, which produces surface reflectance values as high as 1.6 in the visible bands. This is a common problem for spring images at high latitudes.

* Unrealistically low surface reflectance surrounding snow/ice: Related to the above, the AOT retrieval over snow/ice is generally too high. When this artificially high AOT is used to derive the surface reflectance of the neighboring non-snow pixels, very low surface reflectance will result. These pixels will appear very dark in the visible bands. If the surface reflectance value of a pixel is below -0.2, a NO_DATA value of -9999 is used. In Figure 1, the pixels in front of the glaciers have surface reflectance values that are too low.  

* Unrealistically low reflectance surrounding clouds: Like for snow, the HLS atmospheric correction does not attempt aerosol retrieval over clouds and a default AOT is used instead. But if the cloud detection fails, an artificially high AOT will be retrieved over clouds. If the high AOT is used to derive the surface reflectance of the neighboring cloud-free pixels, very low surface reflectance values will result. If the surface reflectance value of a pixel is below -0.2, a NO_DATA value of -9999 is used. 

* Unusually low reflectance around other bright land targets: While the HLS atmospheric correction retrieves AOT over non-cloud, non-snow bright pixels, the retrieved AOT over bright targets can be unrealistically high in some cases, similar to cloud or snow. If this unrealistically high AOT is used to derive the surface reflectance of the neighboring pixels, very low surface reflectance values can result as shown in Figure 2. If the surface reflectance value of a pixel is below -0.2, a NO_DATA value of -9999 is used. These types of bright targets are mostly man-made, such as buildings, parking lots, and roads.  

* Dark plumes over water: The HLS atmospheric correction does not attempt aerosol retrieval over water. For water pixels, the AOT retrieved from the nearest land pixels is used to derive the surface reflectance, but if the retrieval is incorrect, e.g. from a cloud pixel, this high AOT will create dark stripes over water, as shown in Figure 3. This happens more often over large water bodies, such as lakes and bays, than over narrow rivers.  

* Landsat WRS-2 Path/Row boundary in L30 reflectance: HLS performs atmospheric correction on Landsat Level 1 images in the original Worldwide Reference System 2 (WRS2) path/row before the derived surface reflectance is reprojected into Military Grid Reference System (MGRS) tiles. If a WRS-2 Landsat image is very cloudy, the AOT from a few remaining clear pixels might be used for the atmospheric correction of the entire image. The AOT that is used can be quite different from the value for the adjacent row in the same path, which results in an artificial abrupt change from one row to the next, as shown in Figure 4. This occurrence is very rare.  
 
* Landsat WRS2 path/row boundary in cloud masks: The cloud mask algorithm Fmask creates mask labels by applying thresholds to the histograms of some metrics for each path/row independently. If two adjacent rows in the same path have distinct distributions within the metrics, abrupt changes in masking patterns can appear across the row boundary, as shown in Figure 5. This occurrence is very rare.  

* Fmask configuration was deficient for 2-3 months in 2021: The HLS installation of Fmask failed to include auxiliary digital elevation model (DEM) and European Space Agency (ESA) Global Surface Water Occurrence data for a 2-3 month run in 2021. This impacted the masking results over water and in mountainous regions. 

* The reflectance “scale_factor” and “offset” for some L30 and S30 bands were not set: The HLS reflectance scaling factor is 0.0001 and offset is 0. However, this information was not set in the Cloud Optimized GeoTIFF (COG) files of some bands for a small number of granules. The lack of this information creates a problem for automatic conversion of the reflectance data, requiring explicit scaling in applications. The problem has been corrected, but the affected granules have not been reprocessed. 

* Incomplete map projection information: For a time, HLS imagery was produced with an incomplete coordinate reference system (CRS). The metadata contains the Universal Transverse Mercator (UTM) zone and coordinates necessary to geolocate pixels within the image but might not be in a standard form, especially for granules produced early in the HLS mission. As a result, an error will occur in certain image processing packages due to the incomplete CRS. The simplest solution is to update to the latest version of Geospatial Data Abstraction Library (GDAL) and/or rasterio, which use the available information without error. 

* False northing of 10^7 for the L30 angle data: The L30 and S30 products do not use a false northing for the UTM projection, and the angle data are supposed to follow the same convention. However, the L30 angle data incorrectly uses a false northing of 10^7. There is no problem with the angle data itself, but the false northing needs to be set to 0 for it to be aligned with the reflectance.

* L30 from Landsat L1GT scenes: Landsat L1GT scenes were not intended for HLS due to their poor geolocation. However, some scenes made it through screening for a short period of HLS production. L1GT L30 scenes mainly consist of extensive cloud or snow that can be eliminated using the Fmask quality bits layer. Users can also identify an L1GT-originated L30 granule by examining the HLS cmr.xml metadata file.

* The UTC dates in the L30/S30 filenames may not be the local dates: UTC dates are used by ESA and the U.S. Geological Survey (USGS) in naming their Level 1 images, and HLS processing retains this information to name the L30 and S30 products. Landsat and Sentinel-2 overpass eastern Australia and New Zealand around 10AM local solar time, but this area is in either UTC+10:00 or +11:00 zone; therefore, the UTC time for some orbits is in fact near the end of the preceding UTC day. For example, HLS.S30.T59HQS.2016117T221552.v2.0 was acquired in the 22nd hour of day 117 of year 2016 in UTC, but the time was 10:15:52 of day 118 locally. Approximately 100 minutes later HLS.S30.T56JML.2016117T235252.v2.0 was acquired in the next orbit in eastern Australia. 

     This issue also occurs for Landsat. For example, HLS.L30.T59HQS.2016117T221209.v2.0 was acquired on the same day as the 
     first S30 example given above, but both on day 118 of 2016 locally. Adding to the confusion for L30, in the same region, Landsat 
     8 and 9 can each overpass once in one of the two adjacent WRS-2 Paths (91/92/93) over a two-day period on a local calendar, 
     but based on UTC time, the two overpasses can appear to be on the same day. For example, in the following seemingly same- 
    day pair, the second L30 is actually for day 168 locally:  
     HLS.L30.T55GCN.2023167T000407.v2.0  
     HLS.L30.T55GCN.2023167T235747.v2.0  
     Bear in mind, the date peculiarity for the data occurs when the overpass time is during the late hours of a UTC day. 

* The atmospheric ancillary data from the wrong date was used for LaSRC: Related to the above, for eastern Australia and New Zealand, L30 and S30 surface reflectance on certain days was created using the atmospheric ancillary data from a date that was one day too early. The exact geographic extent of the affected HLS products and the impact on the surface reflectance quality are under investigation. Practice caution when using data with overpass times during the late hours of a UTC day.

* Duplicates in L30: The Landsat 9 acquisitions from October 2021 to March 2023 in Landsat Collection 2 were reprocessed by USGS in March 2023. This reprocessing updated the overpass time by a fraction of a second for some scenes. Since HLS uses overpass time as part of the L30 filename, the older L30 granules were not automatically overwritten due to the different filenames. For example, the first L30 granule in the following pair originated from an older version of L1TP of Landsat 9 with the second granule originating from the reprocessed version.  
HLS.L30.T11SLC.2022166T182646.v2.0  
HLS.L30.T11SLC.2022166T182645.v2.0  
There are other causes of duplicate L30 granules, but the overall number of duplicates is very small.

* Poor Geolocation: A large amount of granules that were processed for May through July 2023 were created with L1GT input scenes which were deemed undesirable due to a poor geolocation issue. These granules were removed from the archive. (see the full list of removed [granules](https://lpdaac.usgs.gov/documents/2161/L30_L1GT_granules_May_July_2023.csv))

Improvements/Changes from Previous Versions

* Aerosol QA bits from the USGS Land Surface Reflectance Code (LaSRC) model output have been added into the Function of Mask (Fmask) data layer. The added two bits indicate the aerosol levels: high, medium, low, and climatology aerosol.

Examine a granule¶

Each granule contains the data for a single point in time for an MGRS tile.

In [4]:

import earthaccess
import morecantile

tms = morecantile.tms.get("WebMercatorQuad")

bounds = tms.bounds(62, 44, 7)
xmin, ymin, xmax, ymax = (round(n, 8) for n in bounds)
concept_id = "C2021957295-LPCLOUD"

results = earthaccess.search_data(
    bounding_box=(xmin, ymin, xmax, ymax),
    count=1,
    concept_id=concept_id,
    temporal=("2024-02-11", "2024-02-13"),
)
print("Granules:")
print(results)
print()
print("Example of COGs URL: ")
for link in results[0].data_links(access="direct"):
    print(link)

import earthaccess import morecantile tms = morecantile.tms.get("WebMercatorQuad") bounds = tms.bounds(62, 44, 7) xmin, ymin, xmax, ymax = (round(n, 8) for n in bounds) concept_id = "C2021957295-LPCLOUD" results = earthaccess.search_data( bounding_box=(xmin, ymin, xmax, ymax), count=1, concept_id=concept_id, temporal=("2024-02-11", "2024-02-13"), ) print("Granules:") print(results) print() print("Example of COGs URL: ") for link in results[0].data_links(access="direct"): print(link)

Granules:
[Collection: {'EntryTitle': 'HLS Sentinel-2 Multi-spectral Instrument Surface Reflectance Daily Global 30m v2.0'}
Spatial coverage: {'HorizontalSpatialDomain': {'Geometry': {'GPolygons': [{'Boundary': {'Points': [{'Longitude': -5.71462895, 'Latitude': 48.63303181}, {'Longitude': -5.04818328, 'Latitude': 48.64735841}, {'Longitude': -4.57847134, 'Latitude': 49.64196234}, {'Longitude': -5.76906794, 'Latitude': 49.61959687}, {'Longitude': -5.71462895, 'Latitude': 48.63303181}]}}]}}}
Temporal coverage: {'RangeDateTime': {'BeginningDateTime': '2024-02-11T11:37:33.939Z', 'EndingDateTime': '2024-02-11T11:37:33.939Z'}}
Size(MB): 134.71638679504395
Data: ['https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B12.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B11.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.Fmask.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.SAA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B07.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B04.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B03.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B05.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B02.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B8A.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.VAA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B01.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.SZA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B06.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.VZA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B08.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B09.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B10.tif']]

Example of COGs URL: 
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B12.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B11.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.Fmask.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.SAA.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B07.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B04.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B03.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B05.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B02.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B8A.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.VAA.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B01.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.SZA.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B06.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.VZA.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B08.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B09.tif
s3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B10.tif

Demonstrate assets_for_tile method¶

While rendering xyz tile images, titiler-cmr searches for assets using the assets_for_tile method which converts the xyz tile extent into a bounding box.

In [5]:

from titiler.cmr.backend import CMRBackend
from titiler.cmr.reader import MultiFilesBandsReader

with CMRBackend(reader=MultiFilesBandsReader) as backend:
    assets = backend.assets_for_tile(
        x=62,
        y=44,
        z=7,
        bands_regex="B[0-9][0-9]",
        concept_id=concept_id,
        temporal=("2024-02-11", "2024-02-13"),
    )

print(assets[0])

from titiler.cmr.backend import CMRBackend from titiler.cmr.reader import MultiFilesBandsReader with CMRBackend(reader=MultiFilesBandsReader) as backend: assets = backend.assets_for_tile( x=62, y=44, z=7, bands_regex="B[0-9][0-9]", concept_id=concept_id, temporal=("2024-02-11", "2024-02-13"), ) print(assets[0])

{'url': {'B12': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B12.tif', 'B11': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B11.tif', 'B07': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B07.tif', 'B04': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B04.tif', 'B03': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B03.tif', 'B05': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B05.tif', 'B02': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B02.tif', 'B01': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B01.tif', 'B06': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B06.tif', 'B08': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B08.tif', 'B09': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B09.tif', 'B10': 's3://lp-prod-protected/HLSS30.020/HLS.S30.T30UUV.2024042T113311.v2.0/HLS.S30.T30UUV.2024042T113311.v2.0.B10.tif'}, 'provider': 'LPCLOUD'}

titiler.cmr API documentation¶

In [6]:

from IPython.display import IFrame

IFrame(f"{titiler_endpoint}/api.html", 900, 500)

from IPython.display import IFrame IFrame(f"{titiler_endpoint}/api.html", 900, 500)

Out[6]:

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 [7]:

r = httpx.get(
    f"{titiler_endpoint}/WebMercatorQuad/tilejson.json",
    params=(
        ("concept_id", concept_id),
        # Datetime in form of `start_date/end_date`
        ("datetime", "2024-10-01T00:00:00Z/2024-10-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
        ("bands_regex", "B[0-9][0-9]"),
        # titiler-cmr can work with both Zarr and COG dataset
        # but we need to tell the endpoints in advance which backend
        # to use
        ("backend", "rasterio"),
        # True Color Image B04,B03,B02
        ("bands", "B04"),
        ("bands", "B03"),
        ("bands", "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),
    ),
).json()

print(r)

r = httpx.get( f"{titiler_endpoint}/WebMercatorQuad/tilejson.json", params=( ("concept_id", concept_id), # Datetime in form of `start_date/end_date` ("datetime", "2024-10-01T00:00:00Z/2024-10-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 ("bands_regex", "B[0-9][0-9]"), # titiler-cmr can work with both Zarr and COG dataset # but we need to tell the endpoints in advance which backend # to use ("backend", "rasterio"), # True Color Image B04,B03,B02 ("bands", "B04"), ("bands", "B03"), ("bands", "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), ), ).json() print(r)

{'tilejson': '2.2.0', 'version': '1.0.0', 'scheme': 'xyz', 'tiles': ['https://staging.openveda.cloud/api/titiler-cmr/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?concept_id=C2021957295-LPCLOUD&datetime=2024-10-01T00%3A00%3A00Z%2F2024-10-10T23%3A59%3A59Z&bands_regex=B%5B0-9%5D%5B0-9%5D&backend=rasterio&bands=B04&bands=B03&bands=B02&color_formula=Gamma+RGB+3.5+Saturation+1.7+Sigmoidal+RGB+15+0.35'], 'minzoom': 8, 'maxzoom': 13, 'bounds': [-180.0, -90.0, 180.0, 90.0], 'center': [0.0, 0.0, 8]}

In [8]:

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[8]:

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 [9]:

r = httpx.get(
    f"{titiler_endpoint}/WebMercatorQuad/tilejson.json",
    params=(
        ("concept_id", concept_id),
        # Datetime in form of `start_date/end_date`
        ("datetime", "2024-06-20T00:00:00Z/2024-06-27T23: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
        ("bands_regex", "B[0-9][0-9]"),
        # titiler-cmr can work with both Zarr and COG dataset
        # but we need to tell the endpoints in advance which backend
        # to use
        ("backend", "rasterio"),
        # NDVI
        ("expression", "(B05-B04)/(B05+B04)"),
        # Need red (B04) and nir (B05) for NDVI
        ("bands", "B05"),
        ("bands", "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),
    ),
).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}/WebMercatorQuad/tilejson.json", params=( ("concept_id", concept_id), # Datetime in form of `start_date/end_date` ("datetime", "2024-06-20T00:00:00Z/2024-06-27T23: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 ("bands_regex", "B[0-9][0-9]"), # titiler-cmr can work with both Zarr and COG dataset # but we need to tell the endpoints in advance which backend # to use ("backend", "rasterio"), # NDVI ("expression", "(B05-B04)/(B05+B04)"), # Need red (B04) and nir (B05) for NDVI ("bands", "B05"), ("bands", "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), ), ).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[9]:

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 [10]:

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 [11]:

r = httpx.post(
    f"{titiler_endpoint}/statistics",
    params=(
        ("concept_id", concept_id),
        # Datetime in form of `start_date/end_date`
        ("datetime", "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
        ("bands_regex", "B[0-9][0-9]"),
        # titiler-cmr can work with both Zarr and COG dataset
        # but we need to tell the endpoints in advance which backend
        # to use
        ("backend", "rasterio"),
        # NDVI
        ("expression", "(B05-B04)/(B05+B04)"),
        # Need red (B04) and nir (B05) for NDVI
        ("bands", "B05"),
        ("bands", "B04"),
    ),
    json=geojson,
    timeout=30,
).json()

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

r = httpx.post( f"{titiler_endpoint}/statistics", params=( ("concept_id", concept_id), # Datetime in form of `start_date/end_date` ("datetime", "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 ("bands_regex", "B[0-9][0-9]"), # titiler-cmr can work with both Zarr and COG dataset # but we need to tell the endpoints in advance which backend # to use ("backend", "rasterio"), # NDVI ("expression", "(B05-B04)/(B05+B04)"), # Need red (B04) and nir (B05) for NDVI ("bands", "B05"), ("bands", "B04"), ), json=geojson, timeout=30, ).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": {
          "(B05-B04)/(B05+B04)": {
            "min": -111.0,
            "max": 1.7976931348623157e+308,
            "mean": null,
            "count": 57304.8046875,
            "sum": null,
            "std": null,
            "median": 0.28023598820059,
            "majority": 0.3333333333333333,
            "minority": -111.0,
            "unique": 42086.0,
            "histogram": [
              [
                57796,
                0,
                0,
                0,
                0,
                0,
                0,
                0,
                0,
                2
              ],
              [
                -111.0,
                1.7976931348623158e+307,
                3.5953862697246315e+307,
                5.393079404586948e+307,
                7.190772539449263e+307,
                8.988465674311579e+307,
                1.0786158809173895e+308,
                1.258385194403621e+308,
                1.4381545078898526e+308,
                1.6179238213760842e+308,
                1.7976931348623157e+308
              ]
            ],
            "valid_percent": 100.0,
            "masked_pixels": 0.0,
            "valid_pixels": 57798.0,
            "percentile_2": 0.001119263773162542,
            "percentile_98": 0.5649122807017544
          }
        }
      }
    }
  ]
}

In [ ]: