I previously asked Combine several raster files with heterogeneous bands about an error I was getting trying to combine non-overlapping rasters, and that Q&A surfaced a better approach that warrants a new question.

@JeffreyEvans wrote a very helpful answer to the original question that outlined a better strategy for merging rasters and calculating a mean raster value by polygon. It works well except that for very large rasters and lots of small polygons the memory requirements are intense. Jeffrey suggested a new question thread to outline a different way of looping through the processing steps more efficiently.

Here are the rasters:

1 Answer 1


Extracting raster data using vector features can be complicated by rasters being in tiles and the use of numerous or large features (long lines or large polygons). At first blush it would seem necessary to mosaic your raster tiles into a single raster. However, this can be time consuming with large rasters and, if you do not need the contigious raster for other analysis, not really necessary. It is possible to extract data across multiple raster tiles and combine the data into an identical object type that a single raster extraction would produce.

However, one can run into issues if the extracted data will not fit into RAM. Here is a simple solution that allows for an iterative approach that uses subsets of polygons. The one caveat is that it is possible to have a single polygon whose extracted raster data will not fit into RAM but in this case, one would ask if a polygon that large is necessary and cannot be subset into a smaller regions? If this problem arises with very high resolution raster data, then an entirely different solution needs to be found.

Add requirements


geo.prj <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'  

r <- list.files(getwd(), pattern="tif$")

Here we simulate some polygon data, using a random sample of each raster. Since the raster data is in a geographic projection, I had to use a function that performed on-the-fly reprojections of the data to get the buffers that will act as the polygon data. This particular raster data has a background value of NA rather than zero so, we do not get anything near a n=5000 sample size but, it is still large enough to illustrate the solution.

s <- list() 
  for(i in 1:4){
    rs <- sampleRandom(raster(r[i]), n, sp=TRUE)
                names(rs) <- "value"
    s[[i]] <- rs 
  s <- do.call(rbind, s)
    sp::proj4string(s) <- geo.prj

geo.buffer <- function(x, r) {
  results <- list()  
    for(i in 1:length(x)){
                l <- x[i,]
      p <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
                   l@coords[[2]], l@coords[[1]])
      b <- rgeos::gBuffer(spTransform(l, CRS(p)), width=r, byid=TRUE)
    results[[i]] <- spTransform(b, x@proj4string)
  return( do.call(rbind, results) ) 
b <- geo.buffer(s, 1000)
  b <- as(p, "sf")

Now, here is where the code becomes more relevant. What I am doing here is creating an index list object size that represents groups of polygons, to run at each iteration. The size of the polygon groups is defined by n. In this case we are running 100 polygons per iteration. This allows for controlling RAM but also speeds up the process, especially compared to a one polygon at a time approach in a i in 1:nrow(p) loop. Based on specific data, you can find the sweet spot between speed and memory allocation and change n accordingly.

The Map function allow us to combine the multiple resulting lists created for each raster. The code is written so it can easily take a different statistic than mean. For ~ 1000 polygons with a 5000m buffer it is clocking at ~20 secs

n = 100             # number of observations per iteration
means <- vector()   # empty vector to hold statistic

# create row indexes for iterator
( size <- split(1:nrow(b), ceiling(seq_along(1:nrow(b)) / n)) )

  for(i in 1:length(size)) {
    p.sub <- b[size[[i]],] # subset polygons
    d <- lapply(Map(rbind, exact_extract(raster(r[1]), p.sub), 
                    exact_extract(raster(r[2]), p.sub),
                    exact_extract(raster(r[3]), p.sub),
                    exact_extract(raster(r[4]), p.sub)), na.omit)
    means <- append(means, unlist(lapply(d, FUN=function(x) 
                                  mean(x[,1], na.rm=TRUE))))
Sys.time() - t1

means[1:100]       # first 100 obs
length(means)      # number of obs, should match number of polygons
b$rmean <- means   # assign results to a new column in the polygon data 
  • 1
    From the standpoint of memory efficiency, I'd recommend working with the "summary operations" for exact_extract instead of R functions (fun='mean', etc.). If this is done, data will be extracted for each polygon separately, so there is no need to manually break up the features into groups. Also, the case of a single polygon whose extracted raster data would not fit into RAM is handled by the max_cells_in_memory argument, which lets you process polygons of any size within a fixed amount of RAM.
    – dbaston
    Commented May 18, 2020 at 18:08
  • I am assuming that the subset of "summary operations" are written in C++, correct? So, I imagine that using one of them speeds things up considerable. Normally, I pass custom functions (aside from standard moments) so, I often gloss over these types of built in functions but, this is a nice addition to exact_extract. I will modify my answer accordingly. Commented May 18, 2020 at 21:18
  • 1
    They tend to be faster, and in my opinion have less awkwardness around how to handle cell coverage fractions. For example, exact_extract(r, p, 'sum')) is equivalent to exact_extract(r, p, function(x,c) sum(x*c, na.rm=TRUE))). It's also twice as fast and doesn't have memory limitations. But as you point out, that's not helpful if you need to pass a custom function.
    – dbaston
    Commented May 18, 2020 at 23:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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