I had related situation with 85000 non-overlapping polygons, and a raster dataset with dimensions ~ 10000x10000. In that case I wanted the full distribution of pixel values inside each polygon, rather than just the mean.
An efficient approach (much faster than extract) was: 1) Rasterize the polygons (to a single raster with the same dimensions as the input raster dataset). My polygons had a unique ID value that was burned into the raster. I did this with 'gdal_rasterize' 2) Use the 'crosstab' function in the 'raster' package to cross-tabulate the values from the ID raster with the values on the input raster dataset. This gives the full range of pixel values associated with each polygon ID. 3) From there, it is straightforward to apply whichever function to the pixel values for each ID. -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Extract-function-too-much-data-for-R-tp7584509p7584521.html Sent from the R-sig-geo mailing list archive at Nabble.com. _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo