4

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

enter image description here

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?

  • 4
    It will only take a huge amount of memory if your raster or polygons are very high resolution. You might get some speedup if you crop the raster, which reduces its size, rather than mask it which just sets NA outside the polygon but creates a raster of the same (large) size. Any chance you can make a reproducible example by giving us your data or using sample data from packages? – Spacedman Dec 14 '16 at 11:47
  • 1
    You error has nothing to do with memory. Memory errors is the one thing that R is very clear on. Besides, there is no reason to expect that an apply like function would clear up memory issues. You need to be more mindful in defining objects. The objects; "df", "c" and "reclass" are R functions and could produce very unexpected results in your code. This is likely the source of your error. – Jeffrey Evans Dec 14 '16 at 17:20
  • Thank you all for comments. @Jeffrey, you are right, actually this is a simplified code I posted to help understand the issue, in my script I didn't use "df" but I did use "c" and "reclass". In fact, changing the code resulted in no errors. – Quechua Dec 15 '16 at 15:15
  • @Spacedman, I solved the issue by cropping first and subsequently masking the cropped raster, which saved a lot of time. The DEM has a 20x20 m cell size and covers more than 10000 km2 (~30000x20000 cells), pretty big. I'm further commenting on the answer below – Quechua Dec 15 '16 at 15:15
4

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

imageplot

  • Thank very much Guzman, cropping definitely helped saving time and memory, supposing that memory was an issue. Actually, only cropping was not enough, since crop works on the extent of the polygon, therefore applying the (wrong) reclassification interval on parts adjacent polygons that fall within the target polygon extent. But masking the cropped raster with the polygon resulted in a great amount of time saved (still 5-6 hours for the full loop on 25 polygons though). If you wish to edit your code above adding mask I'll mark it as accepted. – Quechua Dec 15 '16 at 15:23
  • Btw, I'll try the memory usage code you posted on my script – Quechua Dec 15 '16 at 15:25
  • @Quechua you are right about crop! I didn't realise that! I edited the answer as you suggested. Please, check that if you want to improve the speed of the process you can parallelize your code using multiple CPUs. – Guzmán Dec 15 '16 at 15:43
  • I know one day I ought to go to parallel computing... Anyway, code is fine now I guess, maybe consider improving the code by changing the object names as @Jeffrey suggested – Quechua Dec 15 '16 at 16:16

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.