# rasterstats: zonal statistics for geospatial rasters, with rio subcommands

> rasterstats is a Python module for summarizing a raster band over vector geometries, the operation most GIS people describe as zonal statistics, plus interpolated point queries for the case where you have coordinates instead of polygons. It ships as a library, as two rio command line subcommands that read and write GeoJSON features, and it runs on rasterio rather than on GDAL bindings directly.

**perrygeo/python-rasterstats** — Summary statistics of geospatial raster datasets based on vector geometries.

- Repository: https://github.com/perrygeo/python-rasterstats
- Stars: 564 · Forks: 117
- Language: Python
- License: BSD-3-Clause
- Published: 2026-10-22 · Updated: 2026-10-22 · Language: en
- Canonical page: https://hysenlabs.com/projects/perrygeo-python-rasterstats

## One raster band, many vector geometries, one summary each

rasterstats is a Python module for summarizing geospatial raster datasets based on vector geometries. The README states the operation in a single sentence: given a vector layer and a raster band, calculate the summary statistics of each vector geometry.

The example it uses is the canonical one. A polygon vector layer and a digital elevation model raster, and the output is the mean elevation of each polygon. That is zonal statistics, the workhorse of land cover and terrain analysis, and it is tedious to write by hand because the interesting part is the geometry: clipping a raster to an arbitrary polygon, handling pixels that only partly fall inside it, and doing that in the right coordinate reference system.

The library covers two operations. Zonal statistics, which summarises a band over polygons, and interpolated point queries, which sample a band at coordinates. The second one exists because point data is the other half of the workflow, and the README's example of it is a plain GeoJSON Point dictionary with coordinates and a call that returns the value at that location.

The audience is analysts and researchers who already have the rasters and the vectors. The package classifiers say Development Status :: 4 - Beta, with the intended audiences Developers and Science/Research, and the topic keywords are gis, geospatial, geographic, raster, vector and zonal statistics. It is a library and a command line tool, not an application, and it deliberately leaves the visualisation and the report to other tools.

## rasterio underneath, numpy in the middle

The dependency list is short and it explains the implementation.

rasterio at version 1.0 or newer is the raster access layer, which means the library speaks GeoTIFF and the other raster formats through rasterio rather than through GDAL bindings directly. numpy at 1.9 or newer does the arithmetic on the array of pixel values. affine handles the transformation from pixel coordinates to world coordinates, which is what makes the point queries land in the right place. shapely does the vector geometry work, and pyogrio handles vector reading, with fiona offered as an alternative for Python versions below 3.14.

The command line layer is click, with cligj and simplejson for the output formats. The project metadata pins click above 7.1 with version 8.2.1 explicitly excluded, which is the kind of detail that tells you a dependency regression was found and worked around rather than accidentally left in place.

Two optional extras are worth knowing. The progress extra adds tqdm, which gives you a progress bar on long running zonal operations over many geometries, and the fiona extra installs fiona only when python_version is below 3.14. The Python requirement itself is 3.9 or newer, with classifiers for 3.9 through 3.14.

The build system is hatchling, with the version taken dynamically rather than written into the metadata, and the repository keeps the source under src/ with tests/ and docs/ alongside it.

## The rio subcommands, and the pipe-friendly CLI

The command line interfaces to zonalstats and point_query are rio subcommands which read and write geojson features. That is the design decision that makes this usable in a pipeline rather than only in a script.

The README's quick start pipes a vector file into a rio subcommand:

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

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

Two things to note about those lines. The output goes to standard output as GeoJSON features, so a second command can consume it, which is the interoperability the README mentions. And the first command in the pipe is spelled fio in the README, where a GeoJSON handling command would normally be rio cat, so that line reads as a typo in the documentation and is worth checking against the rio tooling you have installed.

The subcommands are registered as plugins of rasterio's own CLI, through the entry point group rasterio.rio_plugins in the project metadata, with zonalstats pointing at rasterstats.cli:zonalstats and pointquery at rasterstats.cli:pointquery. Registering under an existing plugin group rather than shipping its own executable is why the commands appear as rio zonalstats and rio pointquery rather than as a separate rasterstats program.

The documentation for the flags lives at the CLI docs page referenced from the README, and the examples directory in the repository holds scripts you can read instead: simple.py, multiproc.py and benchmark.py.

## Installing the package and taking a first zonal summary

The README points at the hosted documentation for installation and usage details, and the package name in the project metadata is rasterstats, so the install is a pip install of that name:

```bash
pip install rasterstats
```

That pulls in rasterio, numpy, shapely, affine, pyogrio and the command line dependencies. If you want the progress bar on long jobs, or the documentation extras, the names to ask for are the ones declared in the metadata: progress, fiona, docs, test and dev.

The first real use from Python is four lines. The zonal statistics example from the README:

```python
>>> from rasterstats import zonal_stats
>>> stats = zonal_stats("tests/data/polygons.shp", "tests/data/slope.tif")
>>> stats[0].keys()
dict_keys(['min', 'max', 'mean', 'count'])
>>> [f['mean'] for f in stats]
[14.660084635416666, 56.60576171875]
```

So you pass a vector path and a raster path, you get back a list of dictionaries, and the keys in the example are min, max, mean and count. The two means printed are the average value of the raster band inside each of the two polygons, which is the number a user came for.

The point query path is the same shape, with a GeoJSON geometry instead of a layer:

```python
>>> from rasterstats import point_query
>>> point = {'type': 'Point', 'coordinates': (245309.0, 1000064.0)}
>>> point_query(point, "tests/data/slope.tif")
[74.09817594635244]
```

The coordinate pair is in the raster's own coordinate reference system, which is the detail to get right on a first run, since a point in latitude and longitude against a projected raster will land nowhere sensible.

The test data used in these examples lives in the repository under tests/data, so you can reproduce the exact output before pointing the library at your own layers.

## Multiprocessing and benchmarks are in the examples folder

The examples directory is worth a look, because it answers the performance question that the README leaves open.

There are three scripts. simple.py is the obvious one. multiproc.py exists to do the work in parallel, which is the natural question for anyone zonal-stats-ing a large polygon layer, since the per-geometry operation is independent and embarrassingly parallel. And benchmark.py exists to measure it.

That the maintainer shipped a multiprocessing example rather than only a benchmark tells you something about how the library is expected to be used at scale. The parallelism lives in the example rather than in the library itself, which is a defensible choice: a library that spawned its own process pool would be awkward to use inside somebody else's pipeline, whereas an example can be copied and adapted.

The other thing the examples make clear is that the raw output is plain dictionaries, so a parallel version is easy to assemble. Each geometry produces a dictionary of summary values, and there is nothing stateful between them, so the results can be gathered in any order and matched back to the input features by index.

The README's issue guidance is also worth noting, since it tells you what the maintainer considers a complete bug report: a link to the smallest possible raster and vector dataset needed to reproduce the error, the Python code or command that reproduces it, and information about the environment with versions of Python, GDAL and numpy and the system memory. The request for the smallest possible dataset is the sensible one, and the request for GDAL versions tells you that the underlying geospatial stack is part of what gets diagnosed.

## Where rasterstats is the wrong tool

The scope is narrow, and the narrowness is the main thing to understand before you adopt it.

The module covers zonal statistics and interpolated point queries. It is not an analysis package in the wider sense, and it does not attempt to classify, interpolate surfaces, compare dates or model anything. If your question is what the statistics of a raster over these polygons are, it is the right tool. If your question is what changed between two rasters over these polygons, you are outside it.

The summary keys are the second boundary. The example output shows min, max, mean and count, and the function is documented to produce a set of summary statistics rather than a full distribution. A workflow that needs a histogram, a percentile, a majority class or a variance has to compute it differently, and while the library may support additional statistics, the visible example does not enumerate them, so check the documentation for the version you install before designing around a particular key.

The dependency picture carries its own limits. The fiona extra is declared as fiona; python_version < '3.14', so on Python 3.14 you are on pyogrio alone for vector reading. The development status classifier is Beta. The repository has no homepage. And the release history has an odd shape: 0.17.0 in July 2022, 0.18.0 in February 2023, then a long gap before 0.21.0 on 2026-05-23, the same day as the last push to the repository. The code moved recently, but the version line before that was dormant for three years, which is worth weighing before you plan around a release schedule.

And there is the coordinate system question again. The library does not describe a reprojection step, so the raster and the vector layer have to be in the same reference system, and the failure mode when they are not is a set of plausible-looking numbers computed from the wrong pixels rather than an error.

## The alternative: masking the raster yourself with rasterio and numpy

There is a real alternative here, and it is not another package. It is to skip the abstraction and do the work directly.

The direct approach is to open the raster with rasterio, read the band as a numpy array with the window covering the geometry's bounding box, build a mask from the polygon with rasterio's own mask or with shapely containment tests, and then compute the statistics on the masked array. That is roughly what rasterstats does, and you can write it in about twenty lines.

The difference in approach is where the complexity lands. Written by hand, the mask handling, the pixel counting rules for partially covered pixels, the coordinate transformation for point queries and the handling of nodata values all become your problem, and every script that needs the numbers has to reimplement them. Through rasterstats, those rules are in one place and the result is a list of dictionaries with consistent keys, which is what makes the parallel example and the command line subcommands possible in the first place.

The trade is that rasterstats adds a dependency layer between you and rasterio, and it offers a fixed set of statistics rather than whatever array operation you would have written. For a one-off calculation on a small layer, the direct approach is less code and more freedom. For a pipeline, a repeated analysis, or a set of layers that must be summarised consistently, the library is the smaller total amount of work, and the GeoJSON in and out of the rio subcommands is what makes it composable with the rest of a command line workflow.

## Conclusion

rasterstats suits a GIS or remote sensing analyst who needs per-polygon summaries of a raster without writing the geometry masking by hand, and it suits pipeline work because the rio subcommands read and write GeoJSON features that compose with other command line tools. It does not suit a project that needs more than summary statistics per zone, since the module's scope is zonal statistics and interpolated point queries, and it does not suit a Python 3.14 environment that needs Fiona, because the fiona extra is declared only for Python versions below 3.14. Check three things before you rely on it: that the summary keys you need are among those the function returns, that your raster and vector layers are in the same coordinate reference system, since the module does not describe a reprojection step for you, and that the version you install matches your rasterio, because the project depends on rasterio >=1.0 and pins click above 7.1 with 8.2.1 excluded.

## FAQ

### What does rasterstats do?

It summarizes a geospatial raster dataset based on vector geometries, providing zonal statistics for polygons and interpolated point queries for coordinates. Given a vector layer and a raster band, it returns the summary statistics of each vector geometry.

### How do I install rasterstats?

The package name is rasterstats and it installs with pip install rasterstats, pulling in rasterio, numpy, shapely, affine and pyogrio. Python 3.9 or newer is required, and optional extras are declared for progress, fiona, docs, test and dev.

### How do I run zonal statistics from the command line?

The command line interfaces are rio subcommands that read and write GeoJSON features, so you pipe a vector file into rio zonalstats with a raster given by the -r flag, and rio pointquery does the same for point features. They are registered as rasterio.rio_plugins entry points.

### Which statistics does rasterstats return?

The README's example shows the keys min, max, mean and count for the first geometry, with the mean values printed for both polygons. The function returns a list of dictionaries, one per input geometry, so the results can be matched back by index.

### Does rasterstats support parallel processing?

The repository ships examples/multiproc.py alongside examples/benchmark.py and examples/simple.py, so the parallel approach is documented as an example to copy rather than built into the library. Each geometry produces an independent dictionary, which is what makes the parallel version straightforward.

## Sources

- [Issues](https://github.com/perrygeo/python-rasterstats/issues)
- [License: BSD-3-Clause](https://github.com/perrygeo/python-rasterstats/blob/master/LICENSE)
- [perrygeo/python-rasterstats on GitHub](https://github.com/perrygeo/python-rasterstats)
- [README](https://github.com/perrygeo/python-rasterstats/blob/master/README.md)
- [Releases](https://github.com/perrygeo/python-rasterstats/releases)

---

Hysen Labs editorial analysis, written from the project's own repository and release notes. Cite the canonical page: https://hysenlabs.com/projects/perrygeo-python-rasterstats