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~= 5.0
• matplotlib

In [1]:

# !pip install rio-tiler matplotlib

# !pip install rio-tiler matplotlib

In [2]:

import morecantile
from rio_tiler.io import Reader
from rio_tiler.profiles import img_profiles
from rio_tiler.models import ImageData

from matplotlib.pyplot import plot, imshow, subplots

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

In [3]:

# 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 [4]:

?Reader

?Reader

Info¶

Read GDAL/Rasterio dataset metadata

In [5]:

# 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().json(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().json(exclude_none=True)) print(src.dataset.closed)

rasterio dataset:
<open DatasetReader name='https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif' mode='r'>

metadata from rasterio:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -9999.0, 'width': 2000, 'height': 2000, 'count': 1, 'crs': CRS.from_epsg(2056), 'transform': Affine(0.5, 0.0, 2573000.0,
       0.0, -0.5, 1086000.0)}

rio-tiler dataset info:
{"bounds": [7.090624928537461, 45.91605844102823, 7.1035698381384185, 45.925093000254144], "minzoom": 15, "maxzoom": 18, "band_metadata": [["b1", {"STATISTICS_COVARIANCES": "10685.98787505646", "STATISTICS_EXCLUDEDVALUES": "-9999", "STATISTICS_MAXIMUM": "2015.0944824219", "STATISTICS_MEAN": "1754.471184271", "STATISTICS_MINIMUM": "1615.8128662109", "STATISTICS_SKIPFACTORX": "1", "STATISTICS_SKIPFACTORY": "1", "STATISTICS_STDDEV": "103.37305197708"}]], "band_descriptions": [["b1", ""]], "dtype": "float32", "nodata_type": "Nodata", "colorinterp": ["gray"], "count": 1, "height": 2000, "overviews": [2, 4, 8], "width": 2000, "driver": "GTiff", "nodata_value": -9999.0}
True

Statistics¶

Return basic data statistics

In [6]:

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

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

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

['b1']
{'min': 1615.81982421875, 'max': 2015.094482421875, 'mean': 1754.59130859375, 'count': 65536.0, 'sum': 114988896.0, 'std': 103.58233071843753, 'median': 1721.3946533203125, 'majority': 1957.414794921875, 'minority': 1615.81982421875, 'unique': 61645.0, 'histogram': [[10417.0, 15877.0, 9360.0, 6441.0, 5490.0, 4938.0, 4231.0, 3141.0, 3532.0, 2109.0], [1615.81982421875, 1655.747314453125, 1695.6748046875, 1735.6021728515625, 1775.5296630859375, 1815.4571533203125, 1855.3846435546875, 1895.3121337890625, 1935.239501953125, 1975.1669921875, 2015.094482421875]], 'valid_percent': 100.0, 'masked_pixels': 0.0, 'valid_pixels': 65536.0, 'percentile_2': 1626.7143310546876, 'percentile_98': 1987.7415161132812}

Plot Histogram values¶

In [7]:

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

Out[7]:

[<matplotlib.lines.Line2D at 0x1686ec190>]

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

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)

(1, 1024, 1024)

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

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"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"assets: {data.assets}") print(f"dataset stats: {data.dataset_statistics}") # If stored in the original dataset print(type(data.data)) print(type(data.mask))

width: 1024
height: 1024
bands: 1
crs: EPSG:2056
bounds: BoundingBox(left=2573000.0, bottom=1085000.0, right=2574000.0, top=1086000.0)
metadata: {'AREA_OR_POINT': 'Area'}
assets: ['https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif']
dataset stats: [(1615.8128662109, 2015.0944824219)]
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>

Display the data¶

In [10]:

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

<class 'rio_tiler.models.ImageData'>
(1, 1024, 1024)
<class 'numpy.ndarray'>
(1024, 1024, 1)

Out[10]:

<matplotlib.image.AxesImage at 0x16855e4f0>

Multi Spectral Data¶

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

In [11]:

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

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

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

rio-tiler dataset info:
{"bounds": [-74.3095632062702, 40.603994417539994, -74.29151245384847, 40.61775082944064], "minzoom": 14, "maxzoom": 19, "band_metadata": [["b1", {}], ["b2", {}], ["b3", {}], ["b4", {}]], "band_descriptions": [["b1", ""], ["b2", ""], ["b3", ""], ["b4", ""]], "dtype": "uint16", "nodata_type": "None", "colorinterp": ["red", "green", "blue", "undefined"], "count": 4, "height": 5000, "overviews": [2, 4, 8, 16], "width": 5000, "driver": "GTiff"}

In [13]:

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

['b1', 'b2', 'b3', 'b4']

Out[13]:

[<matplotlib.lines.Line2D at 0x168bfe4c0>]

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

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

(1, 1024, 1024)
uint16

Out[14]:

<matplotlib.image.AxesImage at 0x168d35220>