By JPD


2014-10-28 10:13:01 8 Comments

I have a line shapefile representing a road network. I wish to rasterize this data, with the resulting values in the raster showing the total length of the lines that fall within the raster cell.

The data is in the British National Grid projection so the units will be meters.

Ideally I would like to perform this operation using R, and i'm guessing that the rasterize function from the raster package would play a part in achieving this, I just can't work out what the function applied should be.

6 comments

@Sergio 2018-07-10 20:03:27

Let me present you the package vein with several functions to work spatial lines and imports sf and data.table

library(vein)
library(sf)
library(cptcity)
data(net)
netsf <- st_as_sf(net) #Convert Spatial to sf
netsf <- st_transform(netsf, 31983) # Project data
netsf$length_m  <- st_length(netsf)
netsf <- netsf[, "length_m"]
g <- make_grid(netsf, width = 1000) #Creat grid of 1000m spacing with columns id for each feature
# Number of lon points: 12
# Number of lat points: 11

gnet <- emis_grid(netsf, g)
plot(gnet["length_m"])

sf_to_raster <- function(x, column, ncol, nrow){
  x <- sf::as_Spatial(x)
  r <- raster::raster(ncol = ncol, nrow = nrow)
  raster::extent(r) <- raster::extent(x)
  r <- raster::rasterize(x, r, column)
  return(r)
}

rr <- sf_to_raster(gnet, "length_m", 12, 11)
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(pal = 5176), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = lucky(), scales = list(draw = T))
# Colour gradient: neota_flor_apple_green, number: 6165

enter image description here

enter image description here

enter image description here

@derekB 2016-03-02 17:13:15

This may sound a bit naive but if it is a road system then select the roads and save to a clip board then look find a tool that allows you to add a buffer to the clipboard set it to the legal width of the road ie 3 meters +/- remember that the buffer is from center line to edge * 2 i for each side so a 3 meter buffer is actually a 6 meter road from side to side.

@alphabetasoup 2016-03-02 18:42:34

What does road width have to do with road length. This answer does not address the question.

@Matt SM 2015-03-19 17:20:05

Here's yet another approach. It deviates from those already given by using the spatstat package. As far as I can tell, this package has its own version of spatial objects (e.g. im vs. raster objects), but the maptools package allows conversion back and forth between spatstat objects and standard spatial objects.

This approach is taken from this R-sig-Geo post.

require(sp)
require(raster)
require(rgdal)
require(spatstat)
require(maptools)
require(RColorBrewer)

# Load data and transform to UTM
roads <- shapefile('data/TZA_roads.shp')
roadsUTM <- spTransform(roads, CRS("+init=epsg:21037"))

# Need to convert to a line segment pattern object with maptools
roadsPSP <- as.psp(as(roadsUTM, 'SpatialLines'))

# Calculate lengths per cell
roadLengthIM <- pixellate.psp(roadsUTM, dimyx=10)

# Convert pixel image to raster in km
roadLength <- raster(dtanz / 1000, crs=projection(roadsUTM))

# Plot
spplot(rtanz, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", roadsUTM), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))

Imgur

The slowest bit is converting the roads from SpatialLines to a Line Segment Pattern (i.e. spatstat::psp). Once that's done the actual length calculation parts is quite quick, even for much higher resolutions. For example, on my old 2009 MacBook:

system.time(pixellate(tanzpsp, dimyx=10))
#    user  system elapsed 
#   0.007   0.001   0.007

system.time(pixellate(tanzpsp, dimyx=1000))
#    user  system elapsed 
#   0.146   0.032   0.178 

@Matt SM 2015-03-19 17:31:59

Hmmmm...I sure wish those axes weren't in scientific notation. Anyone have any idea how to fix that?

@Jeffrey Evans 2016-03-02 18:25:06

You can modify the global R settings and turn off notation using: options(scipen=999)) however, I do not know if lattice will honor the global environment settings.

@Robert Hijmans 2015-03-18 01:10:14

The below is modified from Jeffrey Evans' solution. This solution is much faster as it does not use rasterize

library(raster)
library(rgdal)
library(rgeos)

roads <- shapefile("TZA_roads.shp")
roads <- spTransform(roads, CRS("+proj=utm +zone=37 +south +datum=WGS84"))
rs <- raster(extent(roads), crs=projection(roads))
rs[] <- 1:ncell(rs)

# Intersect lines with raster "polygons" and add length to new lines segments
rsp <- rasterToPolygons(rs)
rp <- intersect(roads, rsp)
rp$length <- gLength(rp, byid=TRUE) / 1000
x <- tapply(rp$length, rp$layer, sum)
r <- raster(rs)
r[as.integer(names(x))] <- x

@Matt SM 2015-03-19 16:53:07

Seems the most efficient and elegant method of those listed. Also, hadn't seen raster::intersect() before, I like that it combines the attributes of the intersected features, unlike rgeos::gIntersection().

@Jeffrey Evans 2016-03-02 18:19:07

+1 always nice seeing more efficient solutions!

@Doon_Bogan 2016-05-23 11:42:04

@RobertH, I tried using this solution for another problem, where I want to do the same as asked in this thread, but with a very large shapefile of roads (for an entire continent). It seems to have worked, but when I try to do the figure as done by @ fdetsch, I do not get a contiguous grid but just a few colored squares in the picture (see in tinypic.com/r/20hu87k/9).

@user3386170 2018-12-21 19:10:55

And by most efficient... with my sample dataset: run time 0.6sec for this solution vs. 8.25sec for the most upvotes solution.

@Jeffrey Evans 2014-11-04 21:40:48

You do not need a for loop. Just intersect everything at once and then add line lengths to the new line segments using the "SpatialLinesLengths" function in sp. Then, using the raster package rasterize function with the fun=sum argument you can create a raster with the sum of the line length(s) intersecting each cell. Using the above answer and associated data here is code that will generate the same results.

require(rgdal)
require(raster)
require(sp)
require(rgeos)

setwd("D:/TEST/RDSUM")
roads <- readOGR(getwd(), "TZA_roads")
  roads <- spTransform(roads, CRS("+init=epsg:21037"))
    rrst <- raster(extent(roads), crs=projection(roads))

# Intersect lines with raster "polygons" and add length to new lines segments
rrst.poly <- rasterToPolygons(rrst)
  rp <- gIntersection(roads, rrst.poly, byid=TRUE)
    rp <- SpatialLinesDataFrame(rp, data.frame(row.names=sapply(slot(rp, "lines"), 
                               function(x) slot(x, "ID")), ID=1:length(rp), 
                               length=SpatialLinesLengths(rp)/1000) ) 

# Rasterize using sum of intersected lines                            
rd.rst <- rasterize(rp, rrst, field="length", fun="sum")

# Plot results
require(RColorBrewer)
spplot(rd.rst, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", rp), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))

@fdetsch 2014-11-05 08:27:09

First time I see SpatialLinesLengths. Guess it is never too late to learn, thank you (: rasterize takes quite long, though (7 times longer than the upper approach on my machine).

@Jeffrey Evans 2014-11-05 15:27:15

I did notice that rasterize was being slow. However, for large problems a for loop is going to really slow things down and I think that your will see a much more optimized solution with rasterize. Besides the developer of the raster package is striving to make every release more optimized and faster.

@Matt SM 2015-03-14 19:09:32

One potential issue I found with this technique is that the rasterize() function includes all lines that touch a given cell. This results in line segment lengths being counted twice in some cases: once in the cell they're supposed to and once in the adjacent cell that the line endpoint just touches.

@Jeffrey Evans 2015-03-14 21:22:25

Yes, but "rp" is the object that is being rasterized which it the intersection of polygons and points.

@fdetsch 2014-10-29 09:46:53

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))

lines2raster

@philiporlando 2018-11-17 00:18:47

This is awesome! I changed sapply() to pbsapply() and used the cluster argument cl = detectCores()-1. Now I'm able to run this example in parallel!

Related Questions

Sponsored Content

1 Answered Questions

[SOLVED] Rasterize point layer with weighted distance

1 Answered Questions

2 Answered Questions

[SOLVED] Rasterizing vector layer with PostGIS?

1 Answered Questions

[SOLVED] Rasterizing getCover not working in R?

1 Answered Questions

Sponsored Content