Following a recent question, you may want to make use of the functionalities offered by the rgeos package to solve your problem. For reasons of reproducibility, I downloaded a shapefile of Tanzanian roads from DIVA-GIS and put it in my current working directory. For the upcoming tasks, you will need three packages:
- rgdal for general spatial data handling
- raster for rasterization of the shapefile data
- rgeos to check intersection of roads with raster template and calculate road lengths
Consequently, your first lines of could should look like this:
library(rgdal)
library(raster)
library(rgeos)
After that, you need to import the shapefile data. Note that DIVA-GIS shapefiles are distributed in EPSG:4326, so I will project the shapefile to EPSG:21037 (UTM 37S) to deal with meters rather than degrees.
roads <- readOGR(dsn = ".", layer = "TZA_roads")
roads_utm <- spTransform(roads, CRS("+init=epsg:21037"))
For subsequent rasterization, you will need a raster template that covers the spatial extent of your shapefile. The raster template consists of 10 rows and 10 columns by default, thus avoiding too extensive computation times.
roads_utm_rst <- raster(extent(roads_utm), crs = projection(roads_utm))
Now that the template is set up, loop through all cells of the raster (that currently consists of NA values only). By assigning a value of '1' to the current cell and subsequently executing rasterToPolygons
, the resulting shapefile 'tmp_shp' automatically holds the extent of the currently processed pixel. gIntersects
detects whether this extent overlaps with roads. If not, the function will return a value of '0'. Otherwise, the road shapefile is cropped by the current cell and the total length of 'SpatialLines' within that cell is being calculated using gLength
.
lengths <- sapply(1:ncell(roads_utm_rst), function(i) {
tmp_rst <- roads_utm_rst
tmp_rst[i] <- 1
tmp_shp <- rasterToPolygons(tmp_rst)
if (gIntersects(roads_utm, tmp_shp)) {
roads_utm_crp <- crop(roads_utm, tmp_shp)
roads_utm_crp_length <- gLength(roads_utm_crp)
return(roads_utm_crp_length)
} else {
return(0)
}
})
Finally, you can insert the calculated lengths (which are converted to kilometers) into the raster template and visually verify your results.
roads_utm_rst[] <- lengths / 1000
library(RColorBrewer)
spplot(roads_utm_rst, scales = list(draw = TRUE), xlab = "x", ylab = "y",
col.regions = colorRampPalette(brewer.pal(9, "YlOrRd")),
sp.layout = list("sp.lines", roads_utm),
par.settings = list(fontsize = list(text = 15)), at = seq(0, 1800, 200))
vignette('over', package = 'sp')
might help.