Line density function in R equivalent to Line density tool in ArcMap (arcpy)

Question

I need to calculate the magnitude-per-unit area of polylines that fall within a radius around each cell. Essentially I need to calculate a km/km2 road density within a 500m pixel search radius. ArcMap has a quick and easy tool that handles this, but I need a pure R solution.

Here is a link on how line density works: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-line-density-works.htm

And this is how to use it in a python (arcpy) script: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/line-density.htm

I currently execute a backwards approach using raster::focal function, calculating a density of burned in road features. I then convert the km2/km2 output to km/km2.

#Import libraries
 library(raster)
 library(rgdal)
 library(gdalUtils)

#Read-in an already created raster mask (cells are all set to 0) 
mask <- raster("x://path to raster mask...")
#Make a copy of the mask to burn features in, keeping the original untouched
roads_mask <- file.copy(mask, "x://output path ...//roads.tif")

#Read-in road features (shapefile format)
roads_sldf <- readOGR("x://path to shapefile" , "roads")

#Rasterize spatial lines data frame ie. burn road features into mask
#Where road features get a value of 1, mask extent gets a value of 0
roads_raster <- gdalUtils::gdal_rasterize(src_datasource = roads_sldf, 
                dst_filename = "x://output path ...//roads.tif", b = 1,      
                burn = 1, l = "roads", output_Raster = TRUE)

#Run a 1km circular radius density function (be mindful of edge effects)
weight <- raster::focalWeight(roads_raster,1000,type = "circle")
1km_rdDensity <- raster::focal(roads_raster, weight, fun=sum, filename = '', 
                  na.rm=TRUE, pad=TRUE, NAonly=FALSE, overwrite=TRUE)

#Convert km2/km2 road density to km/km2
#Set up the moving window
weight <- raster::focalWeight(roads_raster,1000,type = "circle")
#Count how many records in each column of the moving window are > 0
columnCount <- apply(weight,2,function(x) sum(x > 0))
#Get the sum of the column count
number_of_cells <- sum(columnCount)
#multiply km2/km2 density by number of cells in the moving window
step1 <- roads_raster * number_of_cells
#Rescale step1 output with respect to cell size(30m) and radius of a circle
final_rdDensity <- (step1*0.03)/3.14159265

#Write out final km/km2 road density raster
writeRaster(final_rdDensity,"X://path to output...", datatype = 'FLT4S', overwrite = TRUE)

After some more research I think I may be able to use a kernel function, however I don't want to apply the smoothing algorithm... As well the output is an 'im' object which I would need to write to as a 'tif'

#Import libraries
library(spatstat)
library(rgdal)

#Read-in road features (shapefile format)
roads_sldf <- readOGR("x://path to shapefile" , "roads")

#Convert roads spatial lines data frame to psp object
psp_roads <- as.psp(roads_sldf)
#Apply kernel density, however this is where I am unsure of the arguments
road_density <- spatstat::density.psp(psp_roads, sigma = 0.01, eps = 500)

Cheers.

I like the compactness of the density.psp approach, but I am still trying to understand the output generated — dwiz, Feb 22 '18 at 16:07

score 0 · Answer 1 · answered Jun 02 '22 at 19:22

See this question https://gis.stackexchange.com/questions/138861/calculating-road-density-in-r-using-kernel-density Tried to mark as a duplicate but doesn't work because the other Q is on gis stack exchange

Short answer is use spatstat.geom::pixellate()

I also needed spatstat.geom::as.psp(sf::st_geometry(x)) to convert an sf lines object to the correct format and maptools::as.im.RasterLayer(r) to convert a raster. I was able to convert the result to RasterLayer with raster::raster(pix_res)

Robert Hijmans · Answer 2 · 2022-06-03T04:50:51.620

Perhaps you can use terra::rasterizeGeom which is available in the development version that you can install with install.packages('terra', repos='https://rspatial.r-universe.dev')

Example data

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f) |> as.lines()
r <- rast(v, res=.1)

Solution

x <- rasterizeGeom(v, r, fun="length", "km")

And then use focal sum, but you would not have a perfect circle.

What you could do instead, if your dataset is not too large, is create a circle for each grid cell and use intersect. Something like this:

p <- xyFromCell(r, 1:ncell(r)) |> vect(crs="+proj=longlat")
p$id <- 1:ncell(r)
b <- buffer(p, 10000)
values(v) <- NULL
i <- intersect(v, b)
x <- aggregate(perim(i), list(id=i$id), sum)
r[x$id] <- x[,2]

Line density function in R equivalent to Line density tool in ArcMap (arcpy)

2 Answers2