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
# !pip install rio-tiler matplotlib
import morecantile
from matplotlib.pyplot import imshow, plot, subplots
from rio_tiler.io import Reader
from rio_tiler.models import ImageData
# 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' id='WebMercatorQuad' crs='http://www.opengis.net/def/crs/EPSG/0/3857>, 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`. 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: 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://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":[2573000.0,1085000.0,2574000.0,1086000.0],"crs":"http://www.opengis.net/def/crs/EPSG/0/2056","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"],"scales":[1.0],"offsets":[0.0],"driver":"GTiff","count":1,"width":2000,"height":2000,"overviews":[2,4,8],"nodata_value":-9999.0} 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"].model_dump())
['b1'] {'min': 1615.81982421875, 'max': 2015.094482421875, 'mean': 1754.5912663564086, 'count': 65536.0, 'sum': 114988893.2319336, 'std': 103.58233071842892, 'median': 1721.393310546875, 'majority': 1957.414794921875, 'minority': 1615.81982421875, 'unique': 61645.0, 'histogram': [[10417, 15877, 9360, 6441, 5490, 4938, 4231, 3141, 3532, 2109], [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.7139892578125, 'percentile_98': 1987.7581787109375}
Plot Histogram values¶
# Band 1
plot(meta["b1"].histogram[1][0:-1], meta["b1"].histogram[0])
[<matplotlib.lines.Line2D at 0x138633f10>]
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)
(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
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¶
# 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.ma.MaskedArray'> (1024, 1024, 1)
<matplotlib.image.AxesImage at 0x138a2e550>
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":[544999.99999999,645000.0,549999.99999999,650000.0],"crs":"http://www.opengis.net/def/crs/EPSG/0/6527","band_metadata":[["b1",{}],["b2",{}],["b3",{}],["b4",{}]],"band_descriptions":[["b1",""],["b2",""],["b3",""],["b4",""]],"dtype":"uint16","nodata_type":"None","colorinterp":["red","green","blue","undefined"],"scales":[1.0,1.0,1.0,1.0],"offsets":[0.0,0.0,0.0,0.0],"driver":"GTiff","count":4,"width":5000,"height":5000,"overviews":[2,4,8,16]}
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 0x138bcd5e0>]
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 0x138c95160>