3

Problem

I'm looking for an efficient, memory friendly to calculate total line (vector data) length per raster cell.

The question of line length per cell has already been posed here in 2014 and a few solutions were proposed. I implemented the fastest solution proposed there (kudos to Robert Hijmans) and updated it a little bit to use sf instead of rgeos (see the example below).

My problem is, that my real life data contains a lot of lines (around 40 000) and a fine regular raster (around 700 000 cells), so I ran into the following error:

Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func) :
  rgeos_binpredfunc_prepared: maximum returned dense matrix size exceeded

I have around 200GB RAM available for this task...

Possible Solutions

  1. I guess an efficient solution via terra::intersect or stars::st_intersects is possible, but I didn't manage to make them work. Maybe they could replace raster::intersect. Btw, I'm a big fan of both stars and terra packages...
  2. Some clever way to do it in chunks? Maybe chop up the raster into multiple small rasters and loop through these, combine the resulting small rasters into a big one after that.
  3. Getting somehow rid of step 3.) Make Polygons from raster could help a lot, too.

Combining the lines into one via sf::st_combine did not make it faster...

Example

library("raster")
library("sf")
library("stars") # not used, but potentially helpful for solution
library("terra") # not used, but potentially helpful for solution


# 1.) Create "SpatialLines" object 
# (not sf Multilinestring, because raster::intersect() requires "Spatial")
nc_lines <- system.file("shape/nc.shp", package = "sf") %>%
  st_read() %>%
  st_geometry() %>%
  st_cast("MULTILINESTRING") %>%
  as("Spatial")

# plot(nc_lines)

calc_lengths_per_cell <-  function(n_cell_rows = 10,
                                   n_cell_cols = 20) {

  # 2.) Create raster with n_cell_rows and n_cell_cols on extent of lines
  # (we need "id" layer for matching later, "values" layer will store the lengths later)
  nc_extent <- sf::st_bbox(nc_lines) %>%
    raster::extent()
  nc_raster <- raster::raster(
    nc_extent,
    nrows = n_cell_rows,
    ncols = n_cell_cols,
    crs = crs(nc_lines)
  )
  nc_raster$values <- rep(0, ncell(nc_raster)) # 
  nc_raster$id <- 1:ncell(nc_raster)
  
  # 3.) Make Polygons from raster
  # (I would love to get rid of this step)
  nc_poly <- raster::rasterToPolygons(nc_raster)
  
  # 4.) Calc intersections between lines and polys
  # (This is the computationally most expensive step.)
  nc_lines_intersec <- raster::intersect(nc_lines, nc_poly) %>%
    st_as_sf()
  
  # 5.) Calc lengths of line segments using sf and sum up lengths
  nc_lines_intersec$length <- sf::st_length(nc_lines_intersec)
  nc_line_lengths <- tapply(nc_lines_intersec$length,
                            nc_lines_intersec$id,
                            sum)
  
  # 6.) Assign lengths to raster cells by id
  # This will override "id" and "values" layer
  nc_raster[as.integer(names(nc_line_lengths))] <- nc_line_lengths
  # drop "id" layer, it's unnecessary
  nc_raster <- dropLayer(nc_raster, "id")
  
  
  # plot(nc_raster, main = "Line length per raster cell") +
  #  plot(nc_lines, add =T)
  
  return(nc_raster)
}

plot

system.time(calc_lengths_per_cell())
#   user  system elapsed 
#  0.176   0.002   0.177 
system.time(calc_lengths_per_cell(100, 200))
#   user  system elapsed 
#  5.800   0.001   5.797 

2 Answers 2

3

Here is a concise approach with 'terra'.

Example data:

library(terra)
v <- vect(system.file("shape/nc.shp", package = "sf")) |> as.lines()
r <- rast(v, ncol=20)

Solution:

x <- rasterizeGeom(v, r, "length")    
plot(x); lines(v)

I do not know how well it scales, but that can be fixed if need be.

2
  • Oh wow I didn't see terra::rasterizeGeom() because I was still on terra terra_1.5-31. What a game changer!
    – gosz
    Jun 29, 2022 at 11:14
  • This solution worked like a charm for me - on one core this solution took ~30h for my 40 000 lines on 700 000 raster cells. Memory use was suprisingly low, below 20GB...
    – gosz
    Jul 5, 2022 at 20:48
1

I have a solution that is approximately twice faster using the same idea with terra and grouping with data.table. It is also likely less memory greedy.

calc_lengths_per_cell2 <-  function(n_cell_rows = 10, n_cell_cols = 20) {
  
  nc_lines <- vect(nc_lines)
  nc_raster <- rast(ext(nc_lines), nrows = n_cell_rows, ncols = n_cell_cols, crs = crs(nc_lines))
  nc_raster <- init(nc_raster, "cell")
  nc_poly <- as.polygons(nc_raster)
  
  nc_lines_intersec <- intersect(nc_poly, nc_lines)
  nc_lines_intersec = st_as_sf(nc_lines_intersec)
  nc_lines_intersec$length = st_length(nc_lines_intersec)
  nc_lines_intersec = st_drop_geometry(nc_lines_intersec)

  data.table::setDT(nc_lines_intersec)
  lengths = nc_lines_intersec[, .(length = sum(length)), by = lyr.1]

  nc_raster = init(nc_raster, 0)
  nc_raster[lengths$lyr.1] = lengths$length
  return(nc_raster)
}

Then if it takes too much memory you can crop your lines in chunks of e.g. 100 km and loop through the chunks

1
  • Great solution, I’ll try it out on my real data and let u know…
    – gosz
    Jun 28, 2022 at 16:00

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.