rasterio is a Python package which aims to provide a friendlier API to GDAL than GDAL’s own Python API (which feels very C-like). It is an open source project on GitHub that is created and maintained by mapbox.
Most of the code examples below assume you have imported
rasterio into the current module with:
rasterio's API documentation can be found at https://rasterio.readthedocs.io/en/stable/index.html. Be warned that it is very incomplete (as of November 2019) — there is missing documentation for many
Reading A GeoTIFF
There are two common ways to do this, with or without a context manager.
With a context manager:
with rasterio.open('example.tif') as dataset: pixels = dataset.read() # This will read all bands # Dataset (the file) is closed automatically once you leave the context
Without a context manager:
ds = rasterio.open('example.tif') pixels = ds.read() # You have to remember to close the dataset yourself ds.close()
read() function as used above will read all bands of data from the
.tif file. You can read a specific band by providing a band index to
read(). The band indexes start from
1, just as they do in GDAL. The following example just reads the first band:
pixels = ds.read(1)
You can also open a raster by passing in a
Path object to
open() (Python v3.4 or higher only):
from pathlib import Path file_path = Path('example.tif') ds = rasterio.open(file_path)
Getting Projection Info
The projection information in obtained through the
ds = rasterio.open('example.tif') ds.crs # EPSG:4326
You can get also get the “Well Known Text” (WKT) syntax:
ds.crs.wkt # GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
To get the Affine transformation:
ds.transform # | 0.00, 0.00, 26.04| # | 0.00,-0.00,-15.29| # | 0.00, 0.00, 1.00|
Converting Coordinates To Pixels
rasterio provides the
index() function in the
Dataset class to convert coordinates from the projection space (e.g.
(latitude, longitude) if in WGS 84) of a dataset to
(x, y) pixel coordinates in the image.
lat = [ 10.0, 20.0 ] lng = [ -120.0, -110.0 ] x, y = ds.index(lat, lng)
reproject() does not create the destination array for you, you have to create the array yourself and pass it into the function.
rasterio.reproject( src_array, dst_array, src_transform, src_crs, dst_transform, dst_crs, resampling)
More on masking can be found at https://rasterio.readthedocs.io/en/latest/api/rasterio.mask.html.
rasterio._err.CPLE_AppDefinedError: Too many points (10201 out of 10201) failed to transform, unable to compute output bounds.
This error usually occurs if you are trying to reproject an image into a projection space that does not contain the image (e.g. images are in completely different UTM zones).
The documentation for the latest version of
rasterio can be found at https://rasterio.readthedocs.io/en/latest/.