Can we optimize zonal_statistics by spatially grouping polygon sets?

[emlys]

The bulk of the runtime of zonal_statistics is taken up iterating through raster blocks:

for each raster:
    for each disjoint polygon set:
        for each block offset:
            calculate stats

Could we reduce the number of block offsets by only working on the area that's the intersection of the raster and each set of polygons? (Currently, we iterate through the area that's the intersection of the raster and all the polygons).

If so, the difference might be small, because each disjoint polygon set could be spread over the whole area. Would it help to use more, smaller, disjoint polygon sets that each have a smaller bounding box?

I'd need to do some testing to determine if this optimization would work.