The following example converts a line vector layer (road network) into a density raster. It seems to be working for others, but it looks like it's creating data gaps at my resolution (50m).
library(spatstat) library(raster) library(rgdal) library(magrittr) # Set common CRS crs <- "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs " # Read in our vector & raster data study_area <- readRDS(gzcon(url('http://web.pdx.edu/~porlando/study_area.RDS'))) %>% spTransform(CRSobj = crs) %>% raster(res = 50, crs = crs) roads <- readRDS(gzcon(url('http://web.pdx.edu/~porlando/fwy.RDS'))) %>% spTransform(CRSobj = crs) %>% crop(study_area) # Use spatstat to create density raster roadsPSP <- as.psp(as(roads, "SpatialLines")) roadLengthIM <- pixellate.psp(roadsPSP, eps = 50) # meters roadLengthRaster <- raster(roadLengthIM, crs = projection(roads)) # The resolution isn't exactly 50m x 50m so I resample again roadLengthRaster50m <- resample(roadLengthRaster, study_area) # Recode values to exaggerate presence/absence values(roadLengthRaster50m)[values(roadLengthRaster50m)<=0] <- NA values(roadLengthRaster50m)[values(roadLengthRaster50m)>0] <- 100 # Plot raster layer on top of vector layer plot(roads, col="black") plot(roadLengthRaster50m, add = T)
There are many sections where the black vector layer is poking through the rasterized version. At this point, I'm not sure whether it's a visualization issue, or a limitation of this approach entirely!
Edit: I've followed all of the various solutions found here and I'm getting similar results, however, the non
spatstat methods take significantly longer! I'm open to trying any library that gets the job done though.