3

I would like to extract land cover data in a buffer of 10 km around objects of class SpatialLines and calculate proportion of each land cover type. So, I used the function extract (package raster). I also used the function crop to crop my raster. In my raster, I have 10 land cover types. Here is my code:

buf_line <- gBuffer(seg_line, width=10000) # seg_line = Lines objects
ha <-extract(x=data_raster,y=buf_line)
# The proportion of each land cover type must be in columns (one column = one land cover type)
ha_1 <-length(ha[[1]][ha[[1]]==1])/length(ha[[1]])
ha_2 <-length(ha[[1]][ha[[1]]==2])/length(ha[[1]])
ha_3 <-length(ha[[1]][ha[[1]]==3])/length(ha[[1]])
ha_4 <-length(ha[[1]][ha[[1]]==4])/length(ha[[1]])
ha_5 <-length(ha[[1]][ha[[1]]==5])/length(ha[[1]])
ha_6 <-length(ha[[1]][ha[[1]]==6])/length(ha[[1]])
ha_7 <-length(ha[[1]][ha[[1]]==7])/length(ha[[1]])
ha_8 <-length(ha[[1]][ha[[1]]==8])/length(ha[[1]])
ha_9 <-length(ha[[1]][ha[[1]]==9])/length(ha[[1]])
ha_10 <-length(ha[[1]][ha[[1]]==10])/length(ha[[1]])

I used this code in a loop lapply to extract land cover data in a buffer of 10 km around 30 000 lines. The problem is that the extraction around 30 000 lines is very slow.

How can I speed up the processing time ?

  • which part is the slow one? ha <- extract() ? – wittich Dec 9 '15 at 13:47
  • Yes, it's the part ha <-extract(x=data_raster,y=buf_line). Here is the output of the function system.time(): user system elapsed 6.66 0.34 7.00 – Pierre Dec 9 '15 at 15:28
  • How many pixels is data_raster? Can you reduce the resolution without unacceptable loss of accuracy? Or tile it and farm it out to a parallel cluster? – Spacedman Dec 9 '15 at 16:59
  • There are 222293148 cells. dimensions : 14903, 14916, 222293148 (nrow, ncol, ncell) resolution : 25, 25 (x, y) How can I do a parallell cluster ? – Pierre Dec 10 '15 at 19:58
1

I find that, for problems like this, one can speed up things up a bit if you crop the raster to the buffer, coerce to vector/matrix and then perform any calculations.

Also, take a look at table and prop.table for calculating your landcover proportions. Here is a polygon example, which is what you are after given the line buffer. The result here will be a list "prop" that is ordered the same as the source polygon data. You can use lapply or do.call on the list to summarize or collapse results to a data.frame where columns are the classes.

First, create some example data and plot it

library(raster)
library(sp)   
p <- raster(nrow=10, ncol=10)
  p[] <- runif(ncell(p)) * 10
    p <- rasterToPolygons(p, fun=function(x){x > 9})
      r <- raster(nrow=100, ncol=100)
        r[] <- round(runif(ncell(r), 1,5),0) 
plot(r)
  plot(p, add=TRUE, lwd=4) 

Now, we create an empty list "prop" to store results and write a for loop that iterates through the polygons, crops the raster to each subset polygon, uses prop.table to get proportions of the nominal values (eg., landcover class) and writes results to the list object.

prop <- list()  
  for(i in 1:nrow(p)){
     lsub <- p[i, ]
     cr <- raster::crop(r, raster::extent(lsub), snap = "out")
     fr <- raster::rasterize(lsub, cr)
     r.sub <- raster::mask(x = cr, mask = fr)
    prop[[i]] <- prop.table(table(getValues(r.sub))) 
  }
as.data.frame( do.call("rbind", prop) ) 

Add you would need to modify for your analysis would be passing the code a line object and using gBuffer to create the polygon for each line.

If you have the RAM available to process the rasters, an alternative would to be the use the sp class "SpatialPixlesDataFrame" for raster objects and the over function. You can import a raster, to a sp class raster, from disk using readGDAL in the rgdal package. Starting with our previous example here is what code would look like.

r <- as(r, "SpatialPixelsDataFrame") # coerce r into a sp object  
prop <- list()  
  for(i in 1:nrow(p)){
     lsub <- p[i, ]
     r.sub <- r[!is.na(sp::over(r, sp::geometry(lsub))), ]
     prop[[i]] <- prop.table(table(r.sub@data[,1])) 
  } 
  • Thank you very much Jeffrey for your help. I tested your code. For one line object, the elapsed time is 5.29 s (instead of 7.00). user system elapsed 5.26 0.03 5.29. It's still a little too long. It's the function rasterize that takes more time (4.24). – Pierre Dec 9 '15 at 18:14
  • I think that this is just a reality that you have to live with until everything is moved into a compiled language. If you have the computing power available you could play with threading the problem. If you read the raster vignette there are examples of writing functions that use multiple cores. I came from a period where we had to write all of our own code and what would now be considered very small problems took hours, if not days, to solve. What can I say, proper analysis still takes time. – Jeffrey Evans Dec 9 '15 at 19:32
  • Is there an alternative function for rasterize in R ? – Pierre Dec 11 '15 at 20:59
  • If you have the RAM avaliable, one thing you could try would be to read your raster as a SpatialPixelsDataFrame (using readGDAL in rgdal) and use "over" on each subset polygon. This would avoid rasterizing. – Jeffrey Evans Dec 11 '15 at 21:47
  • Update 4/18/2019, you can use fasterize for much faster rasterization. – Jeffrey Evans Apr 18 at 22:25

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.