How to identify raster cells that intersect with lines in R?-CodePudding

Problem

I have Spatial Lines (as an terra::vect() object) that I would like to intersect with a Raster (terra::rast()) in order to identify the cells by id that each line intersects.

As an output I want a long form data.frame with columns line_id and cell_id (see example below).

Possible solutions

I would prefer a terra solution, but I'm also fine with a solution using sf, stars, raster, etc. packages.

Vectorising the raster would be an option, but my real life data is pretty big - I might easily run into performance issues. So I would like to get away without vectorising the raster.

I thought about terra::intersect(), but it does not take SpatRasters as an argument. And terra::rasterizeGeom()` does not return ids, but calculates length, etc...

Example

library("terra")

# Create example lines
v <- vect(system.file("ex/lux.shp", package="terra"))
lns <- as.lines(v)[c(1,7), "ID_1"] # filter for two lines only & only keep ID
names(lns) <- "line_id"
lns
# class       : SpatVector 
# geometry    : lines 
# dimensions  : 2, 1  (geometries, attributes)
# extent      : 5.826232, 6.381701, 49.46498, 50.18162  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       : line_id
# type        :   <num>
# values      :       1
#                     2

# Create example raster
r <- rast(v, res=.26)
r <- init(r, fun=1:ncell(r))
names(r) <- "cell_id"
r
# class       : SpatRaster 
# dimensions  : 3, 3, 1  (nrow, ncol, nlyr)
# resolution  : 0.26, 0.26  (x, y)
# extent      : 5.74414, 6.52414, 49.44781, 50.22781  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# source      : memory 
# name        : cell_id 
# min value   :       1 
# max value   :       9 

# Plot example
plot(r)
text(r) # add cell_ids as labels to plot
plot(lns, add = T)

For this example I want my output data.frame to look like this:

# data.frame(line_id = c(1,1,1,2,2),
#             cell_id = c(1,2,5,8,9))
#
#     line_id cell_id
# 1       1       1
# 2       1       2
# 3       1       5
# 4       2       8
# 5       2       9

CodePudding user response：

Am I missing something simply making use of terra::extract()?

extract(r, lns)
#>   ID cell_id
#> 1  1       1
#> 2  1       2
#> 3  1       5
#> 4  2       8
#> 5  2       9