# Speeding up extraction time of proportion of land cover types in buffer from raster using R?

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()` ? Commented Dec 9, 2015 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 ` Commented Dec 9, 2015 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? Commented Dec 9, 2015 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 ? Commented Dec 10, 2015 at 19:58

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)
``````

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)
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). Commented Dec 9, 2015 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. Commented Dec 9, 2015 at 19:32
• Is there an alternative function for `rasterize` in R ? Commented Dec 11, 2015 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. Commented Dec 11, 2015 at 21:47
• Update 4/18/2019, you can use fasterize for much faster rasterization. Commented Apr 18, 2019 at 22:25