I'm trying to extract data out of many raster images with different bounding boxes. Using the great rasterstats.zonal_stats it's a good approach, but it takes a considerable amount of time, especially for a large list of geometries. I wrote a small solution that works really fine for a set of raster files with equal bounds and resolution. Nonetheless, extending this to more complex scenarios, where rasters do not necessarily are equal (see image below)

bounding boxes

Case A: raster images with equal bounds and resolution:

Looking at the source code behind rasterstats, the code iterates through each of the geometries passed to zonal_stats and use rasterio.features.rasterize, to mask the raster to the selected geometry and make the needed calculations. Despite being a robust approach it might be overacting when extracting > 100 rasters with similar metadata.

When a raster has the same metadata: same resolution, same boundaries (bounding box), height, width, etc. It is possible to do the following:

1. Create a crosswalk dictionary:

Store all the indices for each of the geometries using one of the rasters sources (again, here I have a list of identical rasters). Hence, the indices will be equal across all files.

with rasterio('path/to/raster/.tif') as src:
    with fiona.open(/path/to/shape, 'r') as shp:
        geoms = [feature['geometry'] for feature in shp]
        index = [feature['properties'][id_var] for feature in shp]

        crosswalk_dict = {}
        for idx, geom in zip(index, geoms):
            geom_rasterize = rasterize([(geom, 1)],

             crosswalk_dict[idx] = np.where(geom_rasterize == 1)

2. Extract data

Loop through all the files and geometries in the crosswalk dictionary and calculate the mean using the source data for each of the files.

for path_raster in list_paths:
    with rasterio.open('/path/to/raster/.tif') as src:
     r_array = src.read(1)
     r_array[r_array == src.nodata] = np.nan

     mean_dict = {}
     for key, value in crosswalk_dict.items():
      mean_dict[key] = np.nanmean(r_array[value])

      df = pd.DataFrame.from_dict(mean_dict,

This process outputs similar results that rasterstats.zonal_stats means and it is quite fast since the rasterization is done only once. Nonetheless, in the image shown above, images have no similar meta. Each of the images has its own affine and boundaries, although they have the same resolution.

Case B: raster images with different bounds:

In this case, I made a similar approach. All the functions in brackets are in this gist. I didn't write them here for space readability issues, they're quite verbose.

  1. Calculate an envelope boundary (the blue box in the image) (draw_envelope)
  2. Give all files the bounds from the envelope built in step 1, but keeping the same data (transform_raster_bounds)
  3. Now with "equal" files, create a crosswalk list using either a meta dict, or a file to create the dictionary (get_crosswalk)
  4. Extract!

This is fine -I believe-, but somehow the results are incorrect, or not comparable to the ones made by rasterstats. I know this is a huge question, but I'm out of ideas on what can be wrong with the transform_raster_bounds function. Also, is this an overkill? There's a way to optimize rasterstats?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.