I have a 3-bands raster and a set of 30 polygons. I would like to get the median value of all pixels inside each polygon (zonal statistics). For that I'm using rasterstats package for zonal statistics.

This is how I do it:

stats=zonal_stats(gdf, url,stats="median")

The result is 1 list with 38 dictionary, each dictionary has one value of median. It looks like this:

[{'median': -6.876018239606002},
 {'median': -4.680290606391631},
 {'median': -6.939485923255775},
 {'median': -2.833199109587449},
 {'median': -4.4440896899689815},
 {'median': -4.47542384471649},...

I was expecting to get 3 values of median, one per band, and instead, I'm getting only 1 value per polygon, and I can't really know from which band the 1 value is taken from.

I couldn't find any information about it in the documentation.

My end goal is to get the median value for each polygon based on the pixels that are inside the polygon. Open also for other packages/method.

1 Answer 1


Found a way to do this, using the "band_num" variable, and running once per bands:

stats1=zonal_stats(gdf, url,stats="median",band_num=1)
stats2=zonal_stats(gdf, url,stats="median",band_num=2)
stats3=zonal_stats(gdf, url,stats="median",band_num=3)
  • 1
    The doc you linked to states "band_num ... defaults to 1" so your original code would always return the median for band 1.
    – user2856
    Jun 26, 2023 at 12:46
  • yes, I missed this part :)
    – ReutKeller
    Jun 26, 2023 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.