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~= 4.0
- matplotlib
# !pip install rio-tiler matplotlib
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
# 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
?Reader
Init signature: Reader( input: str, dataset: Union[rasterio.io.DatasetReader, rasterio.io.DatasetWriter, rasterio.io.MemoryFile, rasterio.vrt.WarpedVRT] = None, tms: morecantile.models.TileMatrixSet = <TileMatrixSet title='Google Maps Compatible for the World' identifier='WebMercatorQuad'>, geographic_crs: rasterio.crs.CRS = CRS.from_epsg(4326), colormap: Dict = None, options: rio_tiler.reader.Options = NOTHING, ) -> None Docstring: Rasterio Reader. Attributes: input (str): dataset path. dataset (rasterio.io.DatasetReader or rasterio.io.DatasetWriter or rasterio.vrt.WarpedVRT, optional): Rasterio dataset. tms (morecantile.TileMatrixSet, optional): TileMatrixSet grid definition. Defaults to `WebMercatorQuad`. geographic_crs (rasterio.crs.CRS, optional): CRS to use as geographic coordinate system. Defaults to WGS84. colormap (dict, optional): Overwrite internal colormap. options (dict, optional): Options to forward to low-level reader methods. Examples: >>> with Reader(src_path) as src: src.tile(...) >>> # Set global options with Reader(src_path, options={"unscale": True, "nodata": 0}) as src: src.tile(...) >>> with rasterio.open(src_path) as src_dst: with WarpedVRT(src_dst, ...) as vrt_dst: with Reader(None, dataset=vrt_dst) as src: src.tile(...) >>> with rasterio.open(src_path) as src_dst: with Reader(None, dataset=src_dst) as src: src.tile(...) Init docstring: Method generated by attrs for class Reader. File: ~/Dev/CogeoTiff/rio-tiler/rio_tiler/io/rasterio.py Type: ABCMeta Subclasses: GCPCOGReader, ImageReader
Info¶
Read GDAL/Rasterio dataset metadata
# 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://njogis-imagery.s3.amazonaws.com/2020/cog/I7D16.tif' mode='r'> metadata from rasterio: {'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 5000, 'height': 5000, 'count': 4, 'crs': CRS.from_epsg(6527), 'transform': Affine(1.0, 0.0, 544999.99999999, 0.0, -1.0, 650000.0)} 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"], "overviews": [2, 4, 8, 16], "height": 5000, "width": 5000, "count": 4, "driver": "GTiff"} True
Statistics¶
Return basic data statistics
with Reader(src_path) as src:
meta = src.statistics(max_size=256)
assert isinstance(meta, dict)
print(list(meta))
print(meta["b1"].dict())
['b1', 'b2', 'b3', 'b4'] {'min': 2074.0, 'max': 63290.0, 'mean': 18474.233474731445, 'count': 65536.0, 'sum': 1210727365.0, 'std': 6623.984653299175, 'median': 17826.0, 'majority': 19834.0, 'minority': 2074.0, 'unique': 23890.0, 'histogram': [[3289.0, 13676.0, 26643.0, 14944.0, 4876.0, 1585.0, 411.0, 67.0, 35.0, 10.0], [2074.0, 8195.6, 14317.2, 20438.800000000003, 26560.4, 32682.0, 38803.600000000006, 44925.200000000004, 51046.8, 57168.4, 63290.0]], 'valid_percent': 100.0, 'masked_pixels': 0.0, 'valid_pixels': 65536.0, 'percentile_2': 6801.7, 'percentile_98': 34989.89999999999}
Plot Histogram values¶
# Band 1
plot(meta["b1"].histogram[1][0:-1], meta["b1"].histogram[0])
[<matplotlib.lines.Line2D at 0x14f996ca0>]
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)
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)
(4, 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
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: {} 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¶
# 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)
<matplotlib.image.AxesImage at 0x14faaad00>
src_path = "https://njogis-imagery.s3.amazonaws.com/2020/cog/I7D16.tif"
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"], "overviews": [2, 4, 8, 16], "height": 5000, "width": 5000, "count": 4, "driver": "GTiff"}
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']
[<matplotlib.lines.Line2D at 0x28f5dc850>]
Using Expression¶
rio-tiler
reader methods accept indexes
option to select the bands you want to read, but also expression
to perform band math.
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
<matplotlib.image.AxesImage at 0x28f6a1af0>