isciences / exactextractr

R package for fast and accurate raster zonal statistics
https://isciences.gitlab.io/exactextractr/
274 stars 26 forks source link

Memory error with coverage_fraction #13

Closed herryATcsiro closed 4 years ago

herryATcsiro commented 4 years ago

Hi there, I am getting the following: Error in CPP_coverage_fraction(x, wkb, crop) : Evaluation error: cannot allocate vector of size 13.6 Mb. Any way of calculating this in chunks to avoid memory issues? Thanks, Herry

dbaston commented 4 years ago

Can you share the code you're running and/or the data used?

The exact_extract function works in chunks to keep a small memory footprint, but coverage_fraction by definition returns an entire raster. If you're computing coverage fraction rasters for an entire set of polygons, you may want to manually loop over the rows in your data frame, computing the coverage fraction for each feature and then writing it to a file.

herryATcsiro commented 4 years ago

Thanks, I'll try more memory (13.6 Mb is not really much additional need) and failing this I'll work on polygon chunks.

This is the code I was using:

poly2rasterArea<-function(r, p){

require(exactextractr)

require(sf)

requires p=polygon as sf, and r = raster

returns a raster with area values for each raster cell

equivalent to the area of the polygon that the pixel covers

cf<-exactextractr::coverage_fraction(r,p) s1<-lapply(1:length(cf), function(x) cf[[x]]res(cf[[x]])[1]res(cf[[x]])[2]) rs<-calc(stack(s1), fun=sum,na.rm=T) rs[rs==0]<-NA return(rs) }

dbaston commented 4 years ago

This is memory inefficient because it will materialize an entire raster for each row in p. An alternative is to combine the features in p into a single geometry and provide that to coverage_fraction:

rast <- raster('gpw_v4_population_count_rev10_2010_30_sec.tif')
poly <- st_read('ne_50m_admin_0_countries.shp')

cell_area <- prod(res(rast))
tot_cov <- coverage_fraction(rast, st_combine(poly))[[1]]*cell_area

Of course, the area calculated in this example is only meaningful if the raster is in projected coordinates (the raster in my example is not.) For geographic coordinates, you could set cell_area <- raster::area(raster), although this will double the memory requirements.