Lapply function on raster iterating through SPDF polygons in R

Question

I have a raster DTM and a SPDF. I'd need to mask a different elevation interval of the DTM within each overlaid polygon.

I tried a for loop to iterate through the polygons, but as one might expect this process requires a huge amount of memory (not to mention the hours) and after 2 or 3 iterations the loop stops for errors

 # open DEM and SPDF
 dem <- raster('dem.tif')
 SPDF <- readOGR(dsn='foo', layer='SPDF')

 # open elevation df
 df <- read_excel('elevation_intervals.xls') # the df is structured as the SPDF@data

 # loop
 for(c in 1:nrow(SPDF)) {
   dem_masked <- mask(dem, SPDF[c,])
   min_elev <- df[c,]$elevation-150
   max_elev <- df[c,]$elevation+150
   reclass <- data.frame(from=as.numeric(c(minValue(dem_masked), min_elev, max_elev)), 
                         to=as.numeric(c(min_elev, max_elev, maxValue(dem_masked))), 
                         becomes=as.numeric(c(NA, 1, NA)))
   dem_masked <- reclassify(dem_masked, reclass)
   assign(paste('DEM_', c, sep=''), dem_masked)
   rm(dem_masked) # I tried to gain memory
   gc() # idem
 } 

This is the error message

 Error in matrix(unlist(ini), ncol = 2, byrow = TRUE) : 
    'data' must be of a vector type, was 'NULL'

and probably relates to memory issues

I guess that using an *apply function would be much effective but I can't find a way how write properly the function since the iteration relates to the argument of mask and not to the x, i.e. the DEM

Otherwise, any other suggestions?

Answer 1

You can try crop and mask instead of only masking over the Raster Layers. Try the reproducible and commented code below. In this example, using crop reduced ~ 92% the size of the objects in the R environment (memory usage). To measure the memory usage I ran the function in this post: Tricks to manage the available memory in an R session. And if you want to improve the speed in the loop you can check this: Processing vector to raster faster with R.

# Load libraries
library('raster')
library('rgdal')

# Get SpatialPolygonsDataFrame object example
SPDF <- getData('GADM', country = 'URY', level = 1)

# Get DEM data example
dem <- getData('alt', country = 'URY')

# Plot data
plot(dem)
plot(SPDF, add = TRUE)

# Make an elevation df as an example

elevation <- numeric()

for (i in 1:nrow(SPDF)) {

    elevation[i] <- sum(range(values(crop(dem, SPDF[i,])), na.rm = TRUE))/2

}

# Add elevation to SPDF
SPDF$elevation <- elevation


# Measure memory usage
# Note: function from https://stackoverflow.com/questions/1358003/tricks-to-manage-the-available-memory-in-an-r-session?rq=1

# improved list of objects
.ls.objects <- function(pos = 1, pattern, order.by, decreasing=FALSE, head=FALSE, n=5) {
  napply <- function(names, fn) sapply(names, function(x)
    fn(get(x, pos = pos)))
  names <- ls(pos = pos, pattern = pattern)
  obj.class <- napply(names, function(x) as.character(class(x))[1])
  obj.mode <- napply(names, mode)
  obj.type <- ifelse(is.na(obj.class), obj.mode, obj.class)
  obj.size <- napply(names, object.size)
  obj.dim <- t(napply(names, function(x)
    as.numeric(dim(x))[1:2]))
  vec <- is.na(obj.dim)[, 1] & (obj.type != "function")
  obj.dim[vec, 1] <- napply(names, length)[vec]
  out <- data.frame(obj.type, obj.size, obj.dim)
  names(out) <- c("Type", "Size", "Rows", "Columns")
  if (!missing(order.by))
    out <- out[order(out[[order.by]], decreasing = decreasing), ]
  if (head)
    out <- head(out, n)
  out
}

# shorthand
lsos <- function(..., n=100) {
  .ls.objects(..., order.by = "Size", decreasing = TRUE, head = TRUE, n = n)
}


TotalMemoryInit <- lsos()
(TotalMemoryInit <- sum(TotalMemoryInit$Size)/1000000)  
# Initial total memory in MB    
# [1] 0.69324

# Masking Loop
for (c in 1:nrow(SPDF)) {

  print(c)
  dem_masked <- mask(dem, SPDF[c,])
  min_elev <- SPDF@data[c,]$elevation - 20
  max_elev <-  SPDF@data[c,]$elevation + 20
  reclass <- data.frame(from = as.numeric(c(minValue(dem_masked), min_elev, max_elev)),
                        to = as.numeric(c(min_elev, max_elev, maxValue(dem_masked))),
                        becomes = as.numeric(c(NA, 1, NA)))
  dem_masked <- reclassify(dem_masked, reclass)
  assign(paste('DEM_mask_', c, sep = ''), dem_masked)
  rm(dem_masked) # I tried to gain memory
  gc() # idem

}


# Memory usage using masking in MB
maskedTotalMemoryMB <- lsos()
(maskedTotalMemoryMB <- sum(maskedTotalMemoryMB$Size)/1000000 - TotalMemoryInit)
# maskedTotalMemoryMB
# [1] 60.50129

# Crop and Mask Loop
for (c in 1:nrow(SPDF)) {

  print(c)
  dem_crop <- crop(dem, SPDF[c,])
  dem_crop_mask <- mask(dem_crop, SPDF[c,])
  min_elev <- SPDF@data[c,]$elevation - 20
  max_elev <-  SPDF@data[c,]$elevation + 20
  reclass <- data.frame(from = as.numeric(c(minValue(dem_crop_mask), min_elev, max_elev)),
                        to = as.numeric(c(min_elev, max_elev, maxValue(dem_crop_mask))),
                        becomes = as.numeric(c(NA, 1, NA)))
  dem_crop <- reclassify(dem_crop_mask, reclass)
  assign(paste('DEM_crop_mask', c, sep = ''), dem_crop_mask)
  rm(dem_crop) # I tried to gain memory
  rm(dem_crop_mask) # I tried to gain memory
  gc() # idem

}

# Memory usage using Crop and Mask in MB
cropAndMaskTotalMemoryMB <- lsos()
(cropAndMaskTotalMemoryMB <- sum(cropAndMaskTotalMemoryMB$Size)/1000000 - TotalMemoryInit - maskedTotalMemoryMB)
# cropAndMaskTotalMemoryMB    
# [1] 2.420832