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~= 4.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 [17]:

import json
import datetime
import httpx

import morecantile

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

from rasterio.features import bounds as featureBounds

import numpy
from matplotlib.pyplot import figure
from ipyleaflet import Map, basemaps, GeoJSON

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

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

# 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"]]

{'page': 1, 'limit': 1000, 'matched': 23, 'returned': 23}

Example:
{
    "assets": {
        "overview": {
            "proj:shape": [
                343,
                343
            ],
            "proj:transform": [
                320,
                0,
                199980,
                0,
                -320,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/L2A_PVI.tif",
            "title": "True color image",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "thumbnail": {
            "href": "https://roda.sentinel-hub.com/sentinel-s2-l1c/tiles/36/R/TT/2019/8/30/0/preview.jpg",
            "title": "Thumbnail",
            "type": "image/png"
        },
        "metadata": {
            "href": "https://roda.sentinel-hub.com/sentinel-s2-l2a/tiles/36/R/TT/2019/8/30/0/metadata.xml",
            "title": "Original XML metadata",
            "type": "application/xml"
        },
        "B11": {
            "proj:shape": [
                5490,
                5490
            ],
            "proj:transform": [
                20,
                0,
                199980,
                0,
                -20,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B11.tif",
            "title": "Band 11 (swir16)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B01": {
            "proj:shape": [
                1830,
                1830
            ],
            "proj:transform": [
                60,
                0,
                199980,
                0,
                -60,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B01.tif",
            "title": "Band 1 (coastal)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B12": {
            "proj:shape": [
                5490,
                5490
            ],
            "proj:transform": [
                20,
                0,
                199980,
                0,
                -20,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B12.tif",
            "title": "Band 12 (swir22)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B02": {
            "proj:shape": [
                10980,
                10980
            ],
            "proj:transform": [
                10,
                0,
                199980,
                0,
                -10,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B02.tif",
            "title": "Band 2 (blue)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B03": {
            "proj:shape": [
                10980,
                10980
            ],
            "proj:transform": [
                10,
                0,
                199980,
                0,
                -10,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B03.tif",
            "title": "Band 3 (green)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B04": {
            "proj:shape": [
                10980,
                10980
            ],
            "proj:transform": [
                10,
                0,
                199980,
                0,
                -10,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B04.tif",
            "title": "Band 4 (red)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "AOT": {
            "proj:shape": [
                1830,
                1830
            ],
            "proj:transform": [
                60,
                0,
                199980,
                0,
                -60,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/AOT.tif",
            "title": "Aerosol Optical Thickness (AOT)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B05": {
            "proj:shape": [
                5490,
                5490
            ],
            "proj:transform": [
                20,
                0,
                199980,
                0,
                -20,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B05.tif",
            "title": "Band 5",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B06": {
            "proj:shape": [
                5490,
                5490
            ],
            "proj:transform": [
                20,
                0,
                199980,
                0,
                -20,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B06.tif",
            "title": "Band 6",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B07": {
            "proj:shape": [
                5490,
                5490
            ],
            "proj:transform": [
                20,
                0,
                199980,
                0,
                -20,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B07.tif",
            "title": "Band 7",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B08": {
            "proj:shape": [
                10980,
                10980
            ],
            "proj:transform": [
                10,
                0,
                199980,
                0,
                -10,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B08.tif",
            "title": "Band 8 (nir)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B8A": {
            "proj:shape": [
                5490,
                5490
            ],
            "proj:transform": [
                20,
                0,
                199980,
                0,
                -20,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B8A.tif",
            "title": "Band 8A",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "B09": {
            "proj:shape": [
                1830,
                1830
            ],
            "proj:transform": [
                60,
                0,
                199980,
                0,
                -60,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/B09.tif",
            "title": "Band 9",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "WVP": {
            "proj:shape": [
                10980,
                10980
            ],
            "proj:transform": [
                10,
                0,
                199980,
                0,
                -10,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/WVP.tif",
            "title": "Water Vapour (WVP)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "visual": {
            "proj:shape": [
                10980,
                10980
            ],
            "proj:transform": [
                10,
                0,
                199980,
                0,
                -10,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/TCI.tif",
            "title": "True color image",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "SCL": {
            "proj:shape": [
                5490,
                5490
            ],
            "proj:transform": [
                20,
                0,
                199980,
                0,
                -20,
                3300000,
                0,
                0,
                1
            ],
            "href": "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/SCL.tif",
            "title": "Scene Classification Map (SCL)",
            "type": "image/tiff; application=geotiff; profile=cloud-optimized"
        },
        "info": {
            "href": "https://roda.sentinel-hub.com/sentinel-s2-l2a/tiles/36/R/TT/2019/8/30/0/tileInfo.json",
            "title": "Original JSON metadata",
            "type": "application/json"
        }
    },
    "bbox": [
        29.896473859714554,
        28.804454491507947,
        31.003792553495717,
        29.81537540150385
    ],
    "geometry": {
        "coordinates": [
            [
                [
                    29.92628490999458,
                    28.804454491507947
                ],
                [
                    29.896473859714554,
                    29.793998377705638
                ],
                [
                    31.003792553495717,
                    29.81537540150385
                ],
                [
                    30.721048631911938,
                    28.8202241961706
                ],
                [
                    29.92628490999458,
                    28.804454491507947
                ]
            ]
        ],
        "type": "Polygon"
    },
    "id": "S2B_36RTT_20190830_0_L2A",
    "collection": "sentinel-s2-l2a-cogs",
    "type": "Feature",
    "properties": {
        "datetime": "2019-08-30T08:52:09Z",
        "eo:cloud_cover": 1.96
    }
}

In [6]:

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

Out[6]:

Map(center=[29.486523152893376, 30.848407745361328], controls=(ZoomControl(options=['position', 'zoom_in_text'…

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

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

4

FirstMethod: Fill with the first value available¶

In [10]:

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)

Out[10]:

<matplotlib.image.AxesImage at 0x176988970>

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

In [11]:

print(len(assets_used))
print(assets_used)

print(len(assets_used)) print(assets_used)

1
['https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/36/R/TT/2019/8/S2B_36RTT_20190830_0_L2A/TCI.tif']

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

In [12]:

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)

Out[12]:

<matplotlib.image.AxesImage at 0x177a1c700>