Get zonal statistics given polygon from 3-bands raster using Rasterstats

Question

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:

stats
>>>
[{'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.

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)