rio-tiler.mosaic¶

In This notebook you'll learn how to:

• use create mercator tiles from multiple observations (assets) using rio_tiler.mosaic submodule
• create custom pixel_selection methods
• look for sentinel-2-cogs data
• create custom tiler using STACReader

Requirements¶

To be able to run this notebook you'll need the following requirements:

• rasterio
• ipyleaflet
• rio-tiler~=7.0
• matplotlib

In [1]:

# !pip install rio-tiler
# !pip install ipyleaflet
# !pip install matplotlib

# !pip install rio-tiler # !pip install ipyleaflet # !pip install matplotlib

In [2]:

import datetime
import json

import httpx
import morecantile
import numpy
from ipyleaflet import GeoJSON, Map, basemaps
from matplotlib.pyplot import figure
from rasterio.features import bounds as featureBounds

from rio_tiler.io import Reader, STACReader
from rio_tiler.models import ImageData
from rio_tiler.mosaic import mosaic_reader
from rio_tiler.mosaic.methods import defaults
from rio_tiler.mosaic.methods.base import MosaicMethodBase

import datetime import json import httpx import morecantile import numpy from ipyleaflet import GeoJSON, Map, basemaps from matplotlib.pyplot import figure from rasterio.features import bounds as featureBounds from rio_tiler.io import Reader, STACReader from rio_tiler.models import ImageData from rio_tiler.mosaic import mosaic_reader from rio_tiler.mosaic.methods import defaults from rio_tiler.mosaic.methods.base import MosaicMethodBase

Data¶

For this demo we will use the Sentinel-2 data stored as COGs on AWS.

Sentinel 2 COGs¶

Thanks to Digital Earth Africa and in collaboration with Sinergise, Element 84, Amazon Web Services (AWS) and the Committee on Earth Observation Satellites (CEOS), Sentinel 2 (Level 2) data over Africa, usually stored as JPEG2000, has been translated to COG more important a STAC database and API has been setup.

https://www.digitalearthafrica.org/news/operational-and-ready-use-satellite-data-now-available-across-africa

The API is provided by @element84 and follows the latest specification: https://earth-search.aws.element84.com/v0

In [3]:

stac_endpoint = "https://earth-search.aws.element84.com/v0/search"

stac_endpoint = "https://earth-search.aws.element84.com/v0/search"

Search for data¶

Define Area of Interest (AOI)¶

In [4]:

# use geojson.io
geojson = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {},
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [
                        [30.810813903808594, 29.454247067148533],
                        [30.88600158691406, 29.454247067148533],
                        [30.88600158691406, 29.51879923863822],
                        [30.810813903808594, 29.51879923863822],
                        [30.810813903808594, 29.454247067148533],
                    ]
                ],
            },
        }
    ],
}

bounds = featureBounds(geojson)

# use geojson.io geojson = { "type": "FeatureCollection", "features": [ { "type": "Feature", "properties": {}, "geometry": { "type": "Polygon", "coordinates": [ [ [30.810813903808594, 29.454247067148533], [30.88600158691406, 29.454247067148533], [30.88600158691406, 29.51879923863822], [30.810813903808594, 29.51879923863822], [30.810813903808594, 29.454247067148533], ] ], }, } ], } bounds = featureBounds(geojson)

STAC Search¶

Use STAC API to search for data over our AOI

Doc: https://github.com/radiantearth/stac-api-spec

In [ ]:

# Date filter
date_min = "2019-06-01"
date_max = "2019-09-01"

start = datetime.datetime.strptime(date_min, "%Y-%m-%d").strftime("%Y-%m-%dT00:00:00Z")
end = datetime.datetime.strptime(date_max, "%Y-%m-%d").strftime("%Y-%m-%dT23:59:59Z")

query = {
    "collections": [
        "sentinel-s2-l2a-cogs"
    ],  # Make sure to query only sentinel-2 COGs collection
    "datetime": f"{start}/{end}",
    "query": {"eo:cloud_cover": {"lt": 5}},  # Use low cloud cover
    "intersects": geojson["features"][0]["geometry"],
    "limit": 1000,
    "fields": {
        "include": [
            "id",
            "properties.datetime",
            "properties.eo:cloud_cover",
        ],  # Make returned response ligth
        "exclude": ["links"],
    },
}

headers = {
    "Content-Type": "application/json",
    "Accept-Encoding": "gzip",
    "Accept": "application/geo+json",
}


data = httpx.post(stac_endpoint, headers=headers, json=query).json()
print(data["context"])
print()
print("Example:")
print(json.dumps(data["features"][0], indent=4))

sceneid = [f["id"] for f in data["features"]]
cloudcover = [f["properties"]["eo:cloud_cover"] for f in data["features"]]
dates = [f["properties"]["datetime"][0:10] for f in data["features"]]

# For this demo we will use the True color image `TCI` asset
assets = [f["assets"]["visual"]["href"] for f in data["features"]]

# Date filter date_min = "2019-06-01" date_max = "2019-09-01" start = datetime.datetime.strptime(date_min, "%Y-%m-%d").strftime("%Y-%m-%dT00:00:00Z") end = datetime.datetime.strptime(date_max, "%Y-%m-%d").strftime("%Y-%m-%dT23:59:59Z") query = { "collections": [ "sentinel-s2-l2a-cogs" ], # Make sure to query only sentinel-2 COGs collection "datetime": f"{start}/{end}", "query": {"eo:cloud_cover": {"lt": 5}}, # Use low cloud cover "intersects": geojson["features"][0]["geometry"], "limit": 1000, "fields": { "include": [ "id", "properties.datetime", "properties.eo:cloud_cover", ], # Make returned response ligth "exclude": ["links"], }, } headers = { "Content-Type": "application/json", "Accept-Encoding": "gzip", "Accept": "application/geo+json", } data = httpx.post(stac_endpoint, headers=headers, json=query).json() print(data["context"]) print() print("Example:") print(json.dumps(data["features"][0], indent=4)) sceneid = [f["id"] for f in data["features"]] cloudcover = [f["properties"]["eo:cloud_cover"] for f in data["features"]] dates = [f["properties"]["datetime"][0:10] for f in data["features"]] # For this demo we will use the True color image `TCI` asset assets = [f["assets"]["visual"]["href"] for f in data["features"]]

In [ ]:

m = Map(
    basemap=basemaps.OpenStreetMap.Mapnik,
    center=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2),
    zoom=8,
)

# add scenes
geo_json = GeoJSON(
    data=data,
    style={"opacity": 1, "dashArray": "1", "fillOpacity": 0, "weight": 1},
)
m.add_layer(geo_json)

# add AOI
geo_json = GeoJSON(
    data=geojson,
    style={"opacity": 1, "dashArray": "1", "fillOpacity": 1, "weight": 1},
)
m.add_layer(geo_json)

m

m = Map( basemap=basemaps.OpenStreetMap.Mapnik, center=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2), zoom=8, ) # add scenes geo_json = GeoJSON( data=data, style={"opacity": 1, "dashArray": "1", "fillOpacity": 0, "weight": 1}, ) m.add_layer(geo_json) # add AOI geo_json = GeoJSON( data=geojson, style={"opacity": 1, "dashArray": "1", "fillOpacity": 1, "weight": 1}, ) m.add_layer(geo_json) m

Define the tiler¶

In [7]:

def tiler(asset, *args, **kwargs):
    with Reader(asset) as src:
        return src.tile(*args, **kwargs)

def tiler(asset, *args, **kwargs): with Reader(asset) as src: return src.tile(*args, **kwargs)

In [ ]:

# List of z12 mercator tiles
tms = morecantile.tms.get("WebMercatorQuad")

tiles = list(tms.tiles(*bounds, 12))
print(len(tiles))

# List of z12 mercator tiles tms = morecantile.tms.get("WebMercatorQuad") tiles = list(tms.tiles(*bounds, 12)) print(len(tiles))

FirstMethod: Fill with the first value available¶

In [ ]:

tile = tiles[0]

img, assets_used = mosaic_reader(
    assets,
    tiler,
    tile.x,
    tile.y,
    tile.z,
    threads=1,
)

fig = figure(figsize=(30, 10))

ax = fig.add_subplot(1, 2, 1)
ax.imshow(img.data_as_image())

ax = fig.add_subplot(1, 2, 2)
ax.imshow(img.mask)

tile = tiles[0] img, assets_used = mosaic_reader( assets, tiler, tile.x, tile.y, tile.z, threads=1, ) fig = figure(figsize=(30, 10)) ax = fig.add_subplot(1, 2, 1) ax.imshow(img.data_as_image()) ax = fig.add_subplot(1, 2, 2) ax.imshow(img.mask)

Print the number and list of assets used to construct the image

In [ ]:

print(len(assets_used))
print(assets_used)

print(len(assets_used)) print(assets_used)

MeanMethod: Get the mean from all the stack of data¶

In [ ]:

tile = tiles[0]

img, assets_used = mosaic_reader(
    assets,
    tiler,
    tile.x,
    tile.y,
    tile.z,
    pixel_selection=defaults.MeanMethod(),
)

fig = figure(figsize=(30, 10))

ax = fig.add_subplot(1, 2, 1)
ax.imshow(img.data_as_image())

ax = fig.add_subplot(1, 2, 2)
ax.imshow(img.mask)

tile = tiles[0] img, assets_used = mosaic_reader( assets, tiler, tile.x, tile.y, tile.z, pixel_selection=defaults.MeanMethod(), ) fig = figure(figsize=(30, 10)) ax = fig.add_subplot(1, 2, 1) ax.imshow(img.data_as_image()) ax = fig.add_subplot(1, 2, 2) ax.imshow(img.mask)

Print the number and list of assets used to construct the image

In [ ]:

print(len(assets_used))
print(assets_used)

print(len(assets_used)) print(assets_used)

Create Custom pixel_selection class¶

1. Using a 4th band for decision: highest NDVI¶

In [13]:

stac_item = "https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a-cogs/items/{sceneid}"
stac_assets = [stac_item.format(sceneid=scene) for scene in sceneid]

stac_item = "https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a-cogs/items/{sceneid}" stac_assets = [stac_item.format(sceneid=scene) for scene in sceneid]

In [ ]:

# Fisrt, let's checkout the STDEV of the NDVI values
# Because we need to use multiple STAC assets, it's easier to use the STACReader


def custom_tiler(asset, *args, **kwargs):
    with STACReader(asset) as stac:
        return stac.tile(*args, expression="(B08_b1-B04_b1)/(B08_b1+B04_b1)")


tile = tiles[0]

img, assets_used = mosaic_reader(
    stac_assets,
    custom_tiler,
    tile.x,
    tile.y,
    tile.z,
    pixel_selection=defaults.StdevMethod(),
)

fig = figure(figsize=(30, 10))

ax = fig.add_subplot(1, 2, 1)
ax.imshow(img.data_as_image())

ax = fig.add_subplot(1, 2, 2)
ax.imshow(img.mask)

# Fisrt, let's checkout the STDEV of the NDVI values # Because we need to use multiple STAC assets, it's easier to use the STACReader def custom_tiler(asset, *args, **kwargs): with STACReader(asset) as stac: return stac.tile(*args, expression="(B08_b1-B04_b1)/(B08_b1+B04_b1)") tile = tiles[0] img, assets_used = mosaic_reader( stac_assets, custom_tiler, tile.x, tile.y, tile.z, pixel_selection=defaults.StdevMethod(), ) fig = figure(figsize=(30, 10)) ax = fig.add_subplot(1, 2, 1) ax.imshow(img.data_as_image()) ax = fig.add_subplot(1, 2, 2) ax.imshow(img.mask)

In [ ]:

print(len(assets_used))
print(assets_used)

print(len(assets_used)) print(assets_used)

In [16]:

# We create a custom tiler function that will read the TCI asset and add a 4th representing the NDVI


def custom_tiler(asset, *args, **kwargs):
    with STACReader(asset) as stac:
        img = stac.tile(*args, assets="visual")
        ndvi = stac.tile(*args, expression="(B08_b1-B04_b1)/(B08_b1+B04_b1)")
        return ImageData(
            numpy.concatenate((img.data, ndvi.data)),
            img.mask,
            crs=img.crs,
            bounds=img.bounds,
        )

# We create a custom tiler function that will read the TCI asset and add a 4th representing the NDVI def custom_tiler(asset, *args, **kwargs): with STACReader(asset) as stac: img = stac.tile(*args, assets="visual") ndvi = stac.tile(*args, expression="(B08_b1-B04_b1)/(B08_b1+B04_b1)") return ImageData( numpy.concatenate((img.data, ndvi.data)), img.mask, crs=img.crs, bounds=img.bounds, )

In [17]:

class CustomFourthBandH(MosaicMethodBase):
    """Feed the mosaic tile with the Mean pixel value."""

    @property
    def data(self):
        """Return data and mask."""
        if self.mosaic is not None:
            return self.mosaic[:-1].copy()

        return None

    def feed(self, array):
        """Add data to mosaic."""

        if self.mosaic is None:
            self.mosaic = array
            return

        pidex = (
            numpy.bitwise_and(array.data[-1] > self.mosaic.data[-1], ~array.mask)
            | self.mosaic.mask
        )

        mask = numpy.where(pidex, array.mask, self.mosaic.mask)
        self.mosaic = numpy.ma.where(pidex, array, self.mosaic)
        self.mosaic.mask = mask

class CustomFourthBandH(MosaicMethodBase): """Feed the mosaic tile with the Mean pixel value.""" @property def data(self): """Return data and mask.""" if self.mosaic is not None: return self.mosaic[:-1].copy() return None def feed(self, array): """Add data to mosaic.""" if self.mosaic is None: self.mosaic = array return pidex = ( numpy.bitwise_and(array.data[-1] > self.mosaic.data[-1], ~array.mask) | self.mosaic.mask ) mask = numpy.where(pidex, array.mask, self.mosaic.mask) self.mosaic = numpy.ma.where(pidex, array, self.mosaic) self.mosaic.mask = mask

In [ ]:

tile = tiles[0]

img, assets_used = mosaic_reader(
    stac_assets,
    custom_tiler,
    tile.x,
    tile.y,
    tile.z,
    pixel_selection=CustomFourthBandH(),
)

fig = figure(figsize=(30, 10))

ax = fig.add_subplot(1, 2, 1)

# NOTE: because we are using NDVI + Visual, the output array, will be in float32
ax.imshow(img.data_as_image().astype("uint8"))

ax = fig.add_subplot(1, 2, 2)
ax.imshow(img.mask)

tile = tiles[0] img, assets_used = mosaic_reader( stac_assets, custom_tiler, tile.x, tile.y, tile.z, pixel_selection=CustomFourthBandH(), ) fig = figure(figsize=(30, 10)) ax = fig.add_subplot(1, 2, 1) # NOTE: because we are using NDVI + Visual, the output array, will be in float32 ax.imshow(img.data_as_image().astype("uint8")) ax = fig.add_subplot(1, 2, 2) ax.imshow(img.mask)

In [ ]:

print(len(assets_used))
print(assets_used)

print(len(assets_used)) print(assets_used)

In [ ]: