0

I'm looking to turn a shapefile with roads (which includes a column of length per road) in the Eastern half of the USA into a raster of 1x1km of road density, using R.

I can't find a straightforward way in Arcmap (Line density works with a radius from the cell center instead of just the cell).

TylerH
  • 20,799
  • 66
  • 75
  • 101
Beardedant
  • 134
  • 8
  • Please update your question with a reproducible example (see [here](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) for tips). If you don't update your question this should be closed for lacking detail – Mako212 Sep 02 '21 at 18:47
  • I can't really do that as I have no code currently, I'm just looking for directions on what package/function to use – Beardedant Sep 02 '21 at 18:58
  • See https://geocompr.robinlovelace.net/ – Phil Sep 03 '21 at 04:46

1 Answers1

1

Here is a solution that creates polygons from the raster cells (adapted from my answer here). You may need to to this for subsets of your dataset and then combine.

Example data

library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))
roads <- as.lines(v)
rs <- rast(v)

Solution

values(rs) <- 1:ncell(rs)
names(rs) <- "rast"    
rsp <- as.polygons(rs)

rp <- intersect(roads, rsp)

rp$length <- perim(rp) / 1000 #km
x <- tapply(rp$length, rp$rast, sum)

r <- rast(rs)
r[as.integer(names(x))] <- as.vector(x)

plot(r)
lines(roads)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thank you, @Robert Hijmans , the example works but because I'm not familiar with Terra I'm unsure how to conver my shapefile (SpatiallinesDataFrame / SP) into a Spatvector. I've been searching the web but don't find an immediate solution... – Beardedant Sep 03 '21 at 15:40
  • 1
    The 2nd line of my code shows the function to use: `vect`. Either with a file or a sp object. – Robert Hijmans Sep 03 '21 at 18:15
  • Thank you @Robert Hijmans ! Definitely making progress, but something's still off in the resolution, I need a very fine (relative to the study area) resolution of 1 sqkm, but res(r) gives me 345410.4x226742.4 (meters, supposedly). What part of the code you gave me controls the resolution, and how could I change the output raster to 1sqkm? – Beardedant Sep 07 '21 at 12:40
  • See ?rast or do something like `res(rs) <- 1000` – Robert Hijmans Sep 10 '21 at 07:26