Most efficient way to retrieve pixel distributions inside many polygones (custom zonal stats)

Question

I want to get the distribution of pixels from a 1 band raster that fall in each polygon of a shapefile to proceed to some statistics. I guess a zonal statistics does it in some kind of way but I need other statitsics that the ones provided by the tools.

Specifically, I want to sort the pixels inside the polygons and retrieve the 50th value (if exists).

Converting my raster to points and joining them spatially to my polygons and exporting my table for the statistical analysis would do it but I will obtain more than 80M points and I estimate that it will take more than 4 days if no issues.

Does anyone know how to perform this in a more efficient way ?

Answer 1

Convert polygons to integer raster, same cell size and extent as value raster, snap raster=v.raster. Use rastertonumyarray on both. Create dictionary. Go through both arrays using integer as key of dictionary and updating list of values from main raster. When done go through dictionary, pick associated list, sort, etc. I've done similar thing on 5000×6000 rasters, think can be done on 9000×9000