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 ?


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

  • I used your method, it worked but I had to split my rasters and my dictionaries to avoid memory error.
    – Delforge
    May 6 '15 at 9:32

Your Answer

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

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