rasterstats

Raster Data Support

rasterstats works with any raster data source supported by rasterio. It supports both:

• Categorical data (e.g., vegetation types)
• Continuous values (e.g., elevation)

Vector Data Support

rasterstats provides flexible support for vector features including Point, LineString, Polygon, or Multi* geometries. Supported inputs include:

• Any fiona data source
• GeoJSON-like mappings
• Objects implementing the geo_interface
• GeoJSON strings
• Well-Known Text/Binary (WKT/WKB) geometries

Any format readable by rasterio is supported. You can provide:

1. File Paths: Direct path to the raster file.
2. Bands: For multi-band rasters, specify the band using the band argument (1-indexed).
3. Numpy Arrays: Pass a numpy.ndarray along with its affine transform mapping dimensions to a coordinate reference system.

import rasterio
from rasterstats import zonal_stats

# Specifying a band
zs = zonal_stats('polygons.shp', 'raster.tif', band=1)

# Using a numpy array and affine transform
with rasterio.open('raster.tif') as src:
    affine = src.transform
    array = src.read(1)
zs = zonal_stats('polygons.shp', array, affine=affine)

You can install rasterstats using uv or pip.

uv add rasterstats

# or

pip install rasterstats

You can install rasterstats using uv (recommended) or pip. Note that rasterstats depends on libgdal, rasterio, fiona, shapely, and numpy.

# Using uv (recommended)
uv add rasterstats

# Using pip
pip install rasterstats

Beyond basic usage, you can find real-world implementations of rasterstats in the following notebooks:

• Integrating with GeoPandas and Numpy: Demonstrates how to combine raster statistics with geospatial dataframes and numerical arrays.
• Bioclimatic Envelope Modeling: Shows usage in ecological modeling, specifically involving False Negative Masking.
• Agricultural zone climate summary: Demonstrates calculating climate summaries for specific agricultural zones.

The point_query function extracts raster values at specific point locations. It accepts a single point geometry (as a GeoJSON-like dictionary) and a raster band (file path).

It returns a list of the raster values found at the provided coordinates.

from rasterstats import point_query
point = {'type': 'Point', 'coordinates': (245309.0, 1000064.0)}
values = point_query(point, "tests/data/slope.tif")
# values is a list of floats, e.g., [74.09817594635244]

For rasters where values represent discrete classes (e.g., land cover), use categorical=True. The output for each feature will be a dictionary where the keys are the unique raster values and the values are the pixel counts.

You can map these raw pixel values to human-readable labels using the category_map argument.

# Basic categorical counts
# Returns: {1.0: 1, 2.0: 9, 5.0: 40}
counts = zonal_stats('polygons.shp', 'classes.tif', categorical=True)

# Using a category map
cmap = {1.0: 'low', 2.0: 'med', 5.0: 'high'}
counts = zonal_stats('polygons.shp', 'classes.tif', categorical=True, category_map=cmap)
# Returns: {'high': 40, 'med': 9, 'low': 1}

The rasterstats command-line interface provides rio subcommands designed to work with GeoJSON features. You can pipe GeoJSON data (e.g., from fio) into these commands to perform spatial analysis.

Available subcommands:

• rio zonalstats: Calculates zonal statistics for the input features using the specified raster.
• rio pointquery: Performs point queries for the input features using the specified raster.

$ fio cat polygon.shp | rio zonalstats -r elevation.tif

$ fio cat points.shp | rio pointquery -r elevation.tif

To inspect the actual pixel data used for a calculation, set raster_out=True. This adds three keys to the resulting dictionary for each feature:

• mini_raster_array: The clipped and masked numpy array.
• mini_raster_affine: The transformation as an Affine object.
• mini_raster_nodata: The nodata value.

Warning: Including large numpy arrays in your output makes it difficult to serialize the results to JSON or other text formats.

stats = zonal_stats('polygons.shp', 'raster.tif', raster_out=True)
# stats[0] will contain 'mini_raster_array', 'mini_raster_affine', etc.

For a fundamental demonstration of how to use rasterstats, refer to the simple.py script located in the examples/ directory of the repository.

# Refer to examples/simple.py for basic usage

If you need to retain the original geometries and properties of your input features, set geojson_out=True. The resulting list will contain GeoJSON Feature objects where the zonal statistics are added to the properties dictionary.

stats = zonal_stats("polygons.shp", "raster.tif", geojson_out=True)
# stats[0] is now a GeoJSON Feature
print(stats[0]['type']) # 'Feature'
print(stats[0]['properties'].keys()) # includes 'id', 'count', 'max', etc.