rasterio backend example: HLS¶

The Harmonized Landsat Sentinel-2 dataset is available in two collections in CMR. This example will use data from the HLSL30.002 (Landsat) dataset.

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 geojson_pydantic
import httpx
import json


from folium import GeoJson, Map, TileLayer

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

/home/henry/workspace/devseed/titiler-cmr/.venv/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

In [2]:

# titiler_endpoint = "http://localhost:8081"  # docker network endpoint
titiler_endpoint = "https://dev-titiler-cmr.delta-backend.com"  # deployed endpoint

# titiler_endpoint = "http://localhost:8081" # docker network endpoint titiler_endpoint = "https://dev-titiler-cmr.delta-backend.com" # deployed endpoint

Identify the dataset¶

You can find the HLSL30.002 dataset using the earthaccess.search_datasets function.

In [3]:

datasets = earthaccess.search_datasets(doi="10.5067/HLS/HLSL30.002")
ds = datasets[0]

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

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

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 and Sentinel-2B 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 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 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.

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)

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) 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 Landsat Operational Land Imager Surface Reflectance and TOA Brightness Daily Global 30m v2.0'}
Spatial coverage: {'HorizontalSpatialDomain': {'Geometry': {'GPolygons': [{'Boundary': {'Points': [{'Longitude': -2.64743819, 'Latitude': 48.6644919}, {'Longitude': -2.21521695, 'Latitude': 49.65006328}, {'Longitude': -3.00027708, 'Latitude': 49.65272281}, {'Longitude': -3.00027162, 'Latitude': 48.66503141}, {'Longitude': -2.64743819, 'Latitude': 48.6644919}]}}]}}}
Temporal coverage: {'RangeDateTime': {'BeginningDateTime': '2024-02-12T11:05:26.302Z', 'EndingDateTime': '2024-02-12T11:05:50.181Z'}}
Size(MB): 56.62721920013428
Data: ['https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B02.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B06.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B01.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.SAA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B07.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.SZA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B03.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.Fmask.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B04.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B05.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.VAA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.VZA.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B11.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B10.tif', 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B09.tif']]

Example of COGs URL: 
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B02.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B06.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B01.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.SAA.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B07.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.SZA.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B03.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.Fmask.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B04.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B05.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.VAA.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.VZA.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B11.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B10.tif
s3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B09.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': {'B02': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B02.tif', 'B06': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B06.tif', 'B01': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B01.tif', 'B07': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B07.tif', 'B03': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B03.tif', 'B04': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B04.tif', 'B05': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B05.tif', 'B11': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B11.tif', 'B10': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B10.tif', 'B09': 's3://lp-prod-protected/HLSL30.020/HLS.L30.T30UWV.2024043T110526.v2.0/HLS.L30.T30UWV.2024043T110526.v2.0.B09.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://dev-titiler-cmr.delta-backend.com/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?concept_id=C2021957657-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 [10]:

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", "viridis"),
        ("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", "viridis"), ("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[10]:

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

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

import json

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))

import json 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": -75.4,
            "max": 26.6,
            "mean": 0.5238783261952482,
            "count": 57304.8046875,
            "sum": 30020.745162633113,
            "std": 0.6052277569586431,
            "median": 0.6041512231282431,
            "majority": 0.75,
            "minority": -75.4,
            "unique": 47613.0,
            "histogram": [
              [
                1,
                0,
                2,
                1,
                0,
                0,
                16,
                57764,
                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": 100.0,
            "masked_pixels": 0.0,
            "valid_pixels": 57798.0,
            "percentile_2": 0.04382638010956595,
            "percentile_98": 0.8685282140779523
          }
        }
      }
    }
  ]
}

In [ ]: