HLS cloud/water-masked max seasonal NDVI¶

TiTiler-CMR provides a really powerful query interface for collections in CMR. In addition to being able to apply nuanced queries during the granule filter phase we can also use rio-tiler's expression and pixel_selection parameters to do complex pixel-wise rendering.

For the Harmonized Landsat Sentinel-2 collection this is really useful because it allows us to make use of the Fmask band which encodes lots of information about every pixel's condition (cloud/cloud shadow/snow/ice/water/etc.).

The following example demonstrates how to use numpy.where and a custom colormap to mask out low quality and water pixels and take the maximum NDVI value from the pixels that remain for the summer of 2024.

At low zoom levels it needs to process a lot of granules and pixels so be patient waiting for the tiles to load!

In [1]:

import json

import httpx
from folium import Map, TileLayer

titiler_cmr_endpoint = "https://v4jec6i5c0.execute-api.us-west-2.amazonaws.com"
collection_concept_id = "C2021957657-LPCLOUD"

# assign this value to pixels that we want to mask out
mask_out = -9999

# identify Fmask value we can use to mask out clouds/shadows/open water
hls_mask_bitfields = [
    1,
    2,
    3,
    5,
]  # cloud shadow, adjacent to cloud shadow, cloud, water
hls_bitmask = 0
for field in hls_mask_bitfields:
    hls_bitmask |= 1 << field

not_cloudy_or_water = f"(Fmask & {hls_bitmask}) == 0"
ndvi = "(B05-B04)/(B05+B04)"

ndvi_colormap = [
    # Water/Snow/Non-vegetated (-1.0 to 0.0)
    ((-10.0, 0.0), (10, 10, 150, 255)),  # Deep Blue
    # Bare Soil/Sand (0.0 to 0.1)
    ((0.0, 0.1), (160, 120, 90, 255)),  # Tan/Brown
    # Sparse Vegetation (0.1 to 0.3)
    ((0.1, 0.3), (190, 210, 100, 255)),  # Yellowish Green
    # Moderate Vegetation (0.3 to 0.6)
    ((0.3, 0.6), (60, 160, 40, 255)),  # Medium Green
    # Dense Vegetation (0.6 to 0.8)
    ((0.6, 0.8), (20, 100, 20, 255)),  # Dark Green
    # Very Dense Canopy (0.8 to 1.0)
    ((0.8, 2.0), (0, 60, 0, 255)),  # Forest Green
]

r = httpx.get(
    f"{titiler_cmr_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-08-30T23:59:59Z"),
        ("sort_key", "-start_date"),
        ("cloud_cover", "0,50"),  # skip the cloudiest granules
        ("exitwhenfull", False),
        ("skipcovered", False),
        ("items_limit", 150),
        ("pixel_selection", "highest"),  # highest NDVI
        ("assets_regex", "B[0-9][0-9]|Fmask"),
        # give me NDVI where it is not cloudy or water based on Fmask
        ("expression", f"where({not_cloudy_or_water},{ndvi}, {mask_out})"),
        ("colormap", json.dumps(ndvi_colormap)),
        # 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", 6),
        ("maxzoom", 13),
    ),
    timeout=None,
).json()

m = Map(location=(40, -110), zoom_start=r["minzoom"] + 1)


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

m

import json import httpx from folium import Map, TileLayer titiler_cmr_endpoint = "https://v4jec6i5c0.execute-api.us-west-2.amazonaws.com" collection_concept_id = "C2021957657-LPCLOUD" # assign this value to pixels that we want to mask out mask_out = -9999 # identify Fmask value we can use to mask out clouds/shadows/open water hls_mask_bitfields = [ 1, 2, 3, 5, ] # cloud shadow, adjacent to cloud shadow, cloud, water hls_bitmask = 0 for field in hls_mask_bitfields: hls_bitmask |= 1 << field not_cloudy_or_water = f"(Fmask & {hls_bitmask}) == 0" ndvi = "(B05-B04)/(B05+B04)" ndvi_colormap = [ # Water/Snow/Non-vegetated (-1.0 to 0.0) ((-10.0, 0.0), (10, 10, 150, 255)), # Deep Blue # Bare Soil/Sand (0.0 to 0.1) ((0.0, 0.1), (160, 120, 90, 255)), # Tan/Brown # Sparse Vegetation (0.1 to 0.3) ((0.1, 0.3), (190, 210, 100, 255)), # Yellowish Green # Moderate Vegetation (0.3 to 0.6) ((0.3, 0.6), (60, 160, 40, 255)), # Medium Green # Dense Vegetation (0.6 to 0.8) ((0.6, 0.8), (20, 100, 20, 255)), # Dark Green # Very Dense Canopy (0.8 to 1.0) ((0.8, 2.0), (0, 60, 0, 255)), # Forest Green ] r = httpx.get( f"{titiler_cmr_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-08-30T23:59:59Z"), ("sort_key", "-start_date"), ("cloud_cover", "0,50"), # skip the cloudiest granules ("exitwhenfull", False), ("skipcovered", False), ("items_limit", 150), ("pixel_selection", "highest"), # highest NDVI ("assets_regex", "B[0-9][0-9]|Fmask"), # give me NDVI where it is not cloudy or water based on Fmask ("expression", f"where({not_cloudy_or_water},{ndvi}, {mask_out})"), ("colormap", json.dumps(ndvi_colormap)), # 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", 6), ("maxzoom", 13), ), timeout=None, ).json() m = Map(location=(40, -110), zoom_start=r["minzoom"] + 1) TileLayer( tiles=r["tiles"][0], opacity=1, attr="NASA", ).add_to(m) m

Out[1]:

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

In [ ]: