Introduction to rio-tiler¶

The goal of this notebook is to give a quick introduction of the main rio-tiler features.

Requirements¶

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

• rio-tiler~=7.0
• matplotlib

In [ ]:

# !pip install rio-tiler matplotlib

# !pip install rio-tiler matplotlib

In [ ]:

import morecantile
from matplotlib.pyplot import imshow, plot, subplots

from rio_tiler.io import Reader
from rio_tiler.models import ImageData

import morecantile from matplotlib.pyplot import imshow, plot, subplots from rio_tiler.io import Reader from rio_tiler.models import ImageData

In [ ]:

# For this DEMO we will use this file
src_path = "https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif"

# For this DEMO we will use this file src_path = "https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif"

rio_tiler.io.COGReader¶

In rio-tiler 2.0 we introduced COGReader (renamed Reader in 4.0), which is a python class providing usefull methods to read and inspect any GDAL/rasterio raster dataset.

Docs: https://cogeotiff.github.io/rio-tiler/readers/#cogreader

In [ ]:

?Reader

?Reader

Info¶

Read GDAL/Rasterio dataset metadata

In [ ]:

# As for Rasterio, using context manager is a good way to
# make sure the dataset are closed when we exit.
with Reader(src_path) as src:
    print("rasterio dataset:")
    print(src.dataset)
    print()
    print("metadata from rasterio:")
    print(src.dataset.meta)
    print()
    # Using rio-tiler Info() method
    info = src.info()
    print("rio-tiler dataset info:")
    print(src.info().model_dump(exclude_none=True))

print(src.dataset.closed)

# As for Rasterio, using context manager is a good way to # make sure the dataset are closed when we exit. with Reader(src_path) as src: print("rasterio dataset:") print(src.dataset) print() print("metadata from rasterio:") print(src.dataset.meta) print() # Using rio-tiler Info() method info = src.info() print("rio-tiler dataset info:") print(src.info().model_dump(exclude_none=True)) print(src.dataset.closed)

Statistics¶

Return basic data statistics

In [ ]:

with Reader(src_path) as src:
    meta = src.statistics(max_size=256)

    assert isinstance(meta, dict)
    print(list(meta))
    print(meta["b1"].model_dump())

with Reader(src_path) as src: meta = src.statistics(max_size=256) assert isinstance(meta, dict) print(list(meta)) print(meta["b1"].model_dump())

Plot Histogram values¶

In [ ]:

# Band 1
plot(meta["b1"].histogram[1][0:-1], meta["b1"].histogram[0])

# Band 1 plot(meta["b1"].histogram[1][0:-1], meta["b1"].histogram[0])

Preview¶

Read a low resolution version of the data (useful when working with COG, because this method will only fetch the overview layer it needs)

In [ ]:

with Reader(src_path) as src:
    # By default `preview()` will return an array with its longest dimension lower or equal to 1024px
    data = src.preview()
    print(data.data.shape)
    assert isinstance(data, ImageData)

with Reader(src_path) as src: # By default `preview()` will return an array with its longest dimension lower or equal to 1024px data = src.preview() print(data.data.shape) assert isinstance(data, ImageData)

The ImageData class¶

To ease data manipulation, rio-tiler version 2.0 uses a new ImageData class that holds the arrays returned by rio-tiler/rasterio low level functions.

Docs: https://cogeotiff.github.io/rio-tiler/models/#imagedata

In [ ]:

print(f"width: {data.width}")
print(f"height: {data.height}")
print(f"bands: {data.count}")
print(f"crs: {data.crs}")
print(f"bounds: {data.bounds}")
print(f"metadata: {data.metadata}")
print(f"descriptions: {data.band_descriptions}")
print(f"assets: {data.assets}")
print(f"dataset stats: {data.dataset_statistics}")  # If stored in the original dataset

print(type(data.data))
print(type(data.mask))

print(f"width: {data.width}") print(f"height: {data.height}") print(f"bands: {data.count}") print(f"crs: {data.crs}") print(f"bounds: {data.bounds}") print(f"metadata: {data.metadata}") print(f"descriptions: {data.band_descriptions}") print(f"assets: {data.assets}") print(f"dataset stats: {data.dataset_statistics}") # If stored in the original dataset print(type(data.data)) print(type(data.mask))

Display the data¶

In [ ]:

# Rasterio doesn't use the same axis order than visualization libraries (e.g matplotlib, PIL)
# in order to display the data we need to change the order (using rasterio.plot.array_to_image).
# the ImageData class wraps the rasterio function in the `data_as_image()` method.
print(type(data))
print(data.data.shape)

image = data.data_as_image()
# data_as_image() returns a numpy.ndarray
print(type(image))
print(image.shape)

imshow(image)

# Rasterio doesn't use the same axis order than visualization libraries (e.g matplotlib, PIL) # in order to display the data we need to change the order (using rasterio.plot.array_to_image). # the ImageData class wraps the rasterio function in the `data_as_image()` method. print(type(data)) print(data.data.shape) image = data.data_as_image() # data_as_image() returns a numpy.ndarray print(type(image)) print(image.shape) imshow(image)

Multi Spectral Data¶

For this demo we will use some High resolution RGB-Nir data hosted on AWS.

In [ ]:

src_path = "https://njogis-imagery.s3.amazonaws.com/2020/cog/I7D16.tif"

src_path = "https://njogis-imagery.s3.amazonaws.com/2020/cog/I7D16.tif"

In [ ]:

with Reader(src_path) as src:
    info = src.info()
    print("rio-tiler dataset info:")
    print(info.model_dump(exclude_none=True))

with Reader(src_path) as src: info = src.info() print("rio-tiler dataset info:") print(info.model_dump(exclude_none=True))

In [ ]:

with Reader(src_path) as src:
    meta = src.statistics()

print(list(meta))

fig, axs = subplots(1, 4, sharey=True, tight_layout=True, dpi=150)
# Red (index 1)
axs[0].plot(meta["b1"].histogram[1][0:-1], meta["b1"].histogram[0])

# Green (index 2)
axs[1].plot(meta["b2"].histogram[1][0:-1], meta["b2"].histogram[0])

# Blue (index 3)
axs[2].plot(meta["b3"].histogram[1][0:-1], meta["b3"].histogram[0])

# Nir (index 3)
axs[3].plot(meta["b4"].histogram[1][0:-1], meta["b4"].histogram[0])

with Reader(src_path) as src: meta = src.statistics() print(list(meta)) fig, axs = subplots(1, 4, sharey=True, tight_layout=True, dpi=150) # Red (index 1) axs[0].plot(meta["b1"].histogram[1][0:-1], meta["b1"].histogram[0]) # Green (index 2) axs[1].plot(meta["b2"].histogram[1][0:-1], meta["b2"].histogram[0]) # Blue (index 3) axs[2].plot(meta["b3"].histogram[1][0:-1], meta["b3"].histogram[0]) # Nir (index 3) axs[3].plot(meta["b4"].histogram[1][0:-1], meta["b4"].histogram[0])

Using Expression¶

rio-tiler reader methods accept indexes option to select the bands you want to read, but also expression to perform band math.

In [ ]:

with Reader(src_path) as src:
    # Return only the third band
    nir_band = src.preview(indexes=4)
    print(nir_band.data.shape)
    print(nir_band.data.dtype)

imshow(nir_band.data_as_image())

with Reader(src_path) as src: # Return only the third band nir_band = src.preview(indexes=4) print(nir_band.data.shape) print(nir_band.data.dtype) imshow(nir_band.data_as_image())

In [ ]:

with Reader(src_path) as src:
    # Return only the third band
    nrg = src.preview(indexes=(4, 3, 1))

    # Data is in Uint16 so we need to rescale
    nrg.rescale(((nrg.data.min(), nrg.data.max()),))

imshow(nrg.data_as_image())

with Reader(src_path) as src: # Return only the third band nrg = src.preview(indexes=(4, 3, 1)) # Data is in Uint16 so we need to rescale nrg.rescale(((nrg.data.min(), nrg.data.max()),)) imshow(nrg.data_as_image())

In [ ]:

with Reader(src_path) as src:
    # Apply NDVI band math
    # (NIR - RED) / (NIR + RED)
    ndvi = src.preview(expression="(b4-b1)/(b4+b1)")
    print(ndvi.data.shape)
    print(ndvi.data.dtype)
    print(ndvi.band_descriptions)
    print("NDVI range: ", ndvi.data.min(), ndvi.data.max())

ndvi.rescale(in_range=((-1, 1),))
imshow(ndvi.data_as_image())

with Reader(src_path) as src: # Apply NDVI band math # (NIR - RED) / (NIR + RED) ndvi = src.preview(expression="(b4-b1)/(b4+b1)") print(ndvi.data.shape) print(ndvi.data.dtype) print(ndvi.band_descriptions) print("NDVI range: ", ndvi.data.min(), ndvi.data.max()) ndvi.rescale(in_range=((-1, 1),)) imshow(ndvi.data_as_image())

Tile¶

Read data for a specific slippy map tile coordinates

In [ ]:

with Reader(src_path) as src:
    print(f"Bounds in dataset CRS: {src.bounds}")
    print(f"Bounds WGS84: {src.get_geographic_bounds('epsg:4326')}")
    print(f"MinZoom (WebMercator): {src.minzoom}")
    print(f"MaxZoom (WebMercator): {src.maxzoom}")

with Reader(src_path) as src: print(f"Bounds in dataset CRS: {src.bounds}") print(f"Bounds WGS84: {src.get_geographic_bounds('epsg:4326')}") print(f"MinZoom (WebMercator): {src.minzoom}") print(f"MaxZoom (WebMercator): {src.maxzoom}")

In [ ]:

# rio-tiler defaults to the WebMercator Grids. The grid definition is provided by the morecantile module
# Docs: https://github.com/developmentseed/morecantile
tms = morecantile.tms.get("WebMercatorQuad")
print(repr(tms))

# Get the list of tiles for the COG minzoom
with Reader(src_path) as cog:
    tile_cover = list(
        tms.tiles(*cog.get_geographic_bounds("epsg:4326"), zooms=cog.minzoom)
    )

print(f"Nb of Z{cog.minzoom} Mercator tiles: {len(tile_cover)}")
print(tile_cover)

# rio-tiler defaults to the WebMercator Grids. The grid definition is provided by the morecantile module # Docs: https://github.com/developmentseed/morecantile tms = morecantile.tms.get("WebMercatorQuad") print(repr(tms)) # Get the list of tiles for the COG minzoom with Reader(src_path) as cog: tile_cover = list( tms.tiles(*cog.get_geographic_bounds("epsg:4326"), zooms=cog.minzoom) ) print(f"Nb of Z{cog.minzoom} Mercator tiles: {len(tile_cover)}") print(tile_cover)

In [ ]:

with Reader(src_path) as src:
    img_1 = src.tile(*tile_cover[0])
    img_1.rescale(((0, 40000),))
    print(img_1.data.shape)

    img_2 = src.tile(*tile_cover[1])
    img_2.rescale(((0, 40000),))

    print(img_2.data.shape)

with Reader(src_path) as src: img_1 = src.tile(*tile_cover[0]) img_1.rescale(((0, 40000),)) print(img_1.data.shape) img_2 = src.tile(*tile_cover[1]) img_2.rescale(((0, 40000),)) print(img_2.data.shape)

In [ ]:

# Show the first 3 bands (RGB)
imshow(img_1.data_as_image()[:, :, 0:3])

# Show the first 3 bands (RGB) imshow(img_1.data_as_image()[:, :, 0:3])

In [ ]:

imshow(img_2.data_as_image()[:, :, 0:3])

imshow(img_2.data_as_image()[:, :, 0:3])

In [ ]:

with Reader(src_path) as src:
    ndvi = src.tile(*tile_cover[0], expression="(b4-b1)/(b4+b1)")
    print(ndvi.data.shape)

ndvi.rescale(in_range=((-1, 1),))
imshow(ndvi.data[0])

with Reader(src_path) as src: ndvi = src.tile(*tile_cover[0], expression="(b4-b1)/(b4+b1)") print(ndvi.data.shape) ndvi.rescale(in_range=((-1, 1),)) imshow(ndvi.data[0])

Part¶

Read data for a given bounding box

In [ ]:

with Reader(src_path) as src:
    # By default `part()` will read the highest resolution. We can limit this by using the `max_size` option.
    img = src.part(
        (-74.30680274963379, 40.60748547709819, -74.29478645324707, 40.61567903099978),
        max_size=1024,
    )
    print("data shape: ", img.data.shape)
    print("bounds: ", img.bounds)
    print("CRS: ", img.crs)

with Reader(src_path) as src: # By default `part()` will read the highest resolution. We can limit this by using the `max_size` option. img = src.part( (-74.30680274963379, 40.60748547709819, -74.29478645324707, 40.61567903099978), max_size=1024, ) print("data shape: ", img.data.shape) print("bounds: ", img.bounds) print("CRS: ", img.crs)

In [ ]:

img.rescale(((0, 40000),))

imshow(img.data_as_image()[:, :, 0:3])

img.rescale(((0, 40000),)) imshow(img.data_as_image()[:, :, 0:3])

Point¶

Read the pixel value for a specific lon/lat coordinate

In [ ]:

with Reader(src_path) as src:
    pt = src.point(-74.30680274963379, 40.60748547709819)

    print("RGB-Nir values:")
    print([(b, pt.data[ii]) for ii, b in enumerate(pt.band_names)])

    print("NDVI values:")
    ndvi = pt.apply_expression("(b4-b1)/(b4+b1)")
    print([(b, ndvi.data[ii]) for ii, b in enumerate(ndvi.band_names)])

with Reader(src_path) as src: pt = src.point(-74.30680274963379, 40.60748547709819) print("RGB-Nir values:") print([(b, pt.data[ii]) for ii, b in enumerate(pt.band_names)]) print("NDVI values:") ndvi = pt.apply_expression("(b4-b1)/(b4+b1)") print([(b, ndvi.data[ii]) for ii, b in enumerate(ndvi.band_names)])

Feature/GeoJSON¶

Read value for a geojson feature defined area

In [ ]:

feat = {
    "type": "Feature",
    "properties": {},
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [-74.30384159088135, 40.614245638811646],
                [-74.30680274963379, 40.61121586776988],
                [-74.30590152740477, 40.608967884350946],
                [-74.30272579193115, 40.60748547709819],
                [-74.29875612258911, 40.60786015456402],
                [-74.2960524559021, 40.61012446497514],
                [-74.29478645324707, 40.61390357476733],
                [-74.29882049560547, 40.61515780103489],
                [-74.30294036865233, 40.61567903099978],
                [-74.3035626411438, 40.61502749290829],
                [-74.30384159088135, 40.614245638811646],
            ]
        ],
    },
}

feat = { "type": "Feature", "properties": {}, "geometry": { "type": "Polygon", "coordinates": [ [ [-74.30384159088135, 40.614245638811646], [-74.30680274963379, 40.61121586776988], [-74.30590152740477, 40.608967884350946], [-74.30272579193115, 40.60748547709819], [-74.29875612258911, 40.60786015456402], [-74.2960524559021, 40.61012446497514], [-74.29478645324707, 40.61390357476733], [-74.29882049560547, 40.61515780103489], [-74.30294036865233, 40.61567903099978], [-74.3035626411438, 40.61502749290829], [-74.30384159088135, 40.614245638811646], ] ], }, }

In [ ]:

with Reader(src_path) as src:
    # we use the feature to define the bounds and the mask
    # but we use `dst_crs` options to keep the projection from the input dataset
    # By default `feature()` will read the highest resolution. We can limit this by using the `max_size` option.
    img = src.feature(feat, dst_crs=src.crs, max_size=1024)
    print("data shape: ", img.data.shape)
    print("bounds: ", img.bounds)
    print("CRS: ", img.crs)

with Reader(src_path) as src: # we use the feature to define the bounds and the mask # but we use `dst_crs` options to keep the projection from the input dataset # By default `feature()` will read the highest resolution. We can limit this by using the `max_size` option. img = src.feature(feat, dst_crs=src.crs, max_size=1024) print("data shape: ", img.data.shape) print("bounds: ", img.bounds) print("CRS: ", img.crs)

In [ ]:

img.rescale(((0, 40000),))
imshow(img.data_as_image()[:, :, 0:3])

img.rescale(((0, 40000),)) imshow(img.data_as_image()[:, :, 0:3])

In [ ]: