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]:
Copied!
# !pip install rio-tiler
# !pip install ipyleaflet
# !pip install matplotlib
# !pip install rio-tiler
# !pip install ipyleaflet
# !pip install matplotlib
In [17]:
Copied!
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.
The API is provided by @element84 and follows the latest specification: https://earth-search.aws.element84.com/v0
In [3]:
Copied!
stac_endpoint = "https://earth-search.aws.element84.com/v0/search"
stac_endpoint = "https://earth-search.aws.element84.com/v0/search"
In [4]:
Copied!
# 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
In [5]:
Copied!
# 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]:
Copied!
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]:
Copied!
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]:
Copied!
# 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]:
Copied!
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]:
Copied!
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]:
Copied!
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>