I have a number of routes (5 in this example data) stored in a SpatialLinesDataFrame called r. These have different origins but all converge on the same destination (Manchester Airport):

lines

Each contains the variable All, representing data on the number of people regularly using each route.

The question is this: how can I merge the lines such that the sections where they overlap are allocated a value of sum(All) for all the lines passing through there?

The red and yellow lines, for example, have All values of 210 and 395, so the orange line going north-south should be a segment with a value of 605:

overlap

> r_overlap$All + r[4,]$All
[1] 605

But I need all the merged lines in one layer, with each new Lines object representing a segment with a single sum(All) value, not just 1 as above.

The data above is stored at raw.githubusercontent.com/npct/pct-data/master/test-data/airport.geojson so the solution can be demonstrable and reproducible. An R-based solution would be preferable to me.

Please find below code used to load and view the data using the new geojsonio and leaflet packages:

pkgs <- c("geojsonio", "leaflet", "sp")
lapply(pkgs, library, character.only = T)
download.file("https://raw.githubusercontent.com/npct/pct-data/master/test-data/airport.geojson", destfile = "l.geojson", method = "wget")
r <- geojsonio::geojson_read("l.geojson")

plot(r) # check the data is there
r@data # look at the data

r_overlap <- gIntersection(r[1,], r[4,])
plot(r)
plot(r_overlap, col = "red", add = T, lwd = 5)

r_overlap$All + r[4,]$All

leaflet() %>%
  addTiles() %>%
  addPolylines(data = r, color = c("red", "blue", "green", "yellow", "orange"))
  • geojsonio is non-CRAN, yes? This: devtools::install_github("ropensci/geojsonio") – Spacedman Mar 20 '15 at 12:49
  • Yes that's right. Equally r <- readOGR("l.geojson", layer = "OGRGeoJSON") works - just playing with new packages - no idea if they're any different in terms of performance - geojson_read is shorter is only difference I see. – RobinLovelace Mar 20 '15 at 13:18
  • So your output is a set of simple linear features corresponding to the sections of road in that network with an attribute that is the total for all the input simple linear features that go over it? I think rgeos::gOverlaps and gIntersection will be the key here... – Spacedman Mar 20 '15 at 14:29
  • Yes I think that's a fair summary of the problem. Any example code? plot(gIntersection(r[1, ], r[2, ])) shows it's pulling out the joint geometry - major step forward. But still cannot see how to generalise the solution. My approach would be to build a double or triple nested for loop that starts at gIntersection(r[1, ], r[2, ]) and ends and gIntersection(r[4, ], r[5, ]) at sums at each stage until there are no more overlaps. Sound like a reasonable approach? If so I'll try it. @Spacedman – RobinLovelace Mar 20 '15 at 14:47
  • 1
    Seems to me this approach is going to require something that gives the opposite of gIntersection, like a gDifference. It's easy enough to get all the overlapping segments with the summed attribute, but then need to merge in all the pieces that didn't overlap (i.e. the result of my hypothetical gDifference function). – Matt SM Mar 20 '15 at 19:27
up vote 4 down vote accepted

Here you go. A couple of utility functions and then the meat in one function (and no for loops :))

islines <- function(g1, g2){
    ## return TRUE if geometries intersect as lines, not points
    inherits(gIntersection(g1,g2),"SpatialLines")
}

sections <- function(sl){
    ## union and merge and disaggregate to make a
    ## set of non-overlapping line segments
    disaggregate(gLineMerge(gUnion(sl,sl)))
}

aggit <- function(sldf, attr, fun=sum){
    ## simplify down to SpatialLines
    sl = as(sldf, "SpatialLines")
    ## get the line sections that make the network
    slu = sections(sl)
    ## overlay network with routes
    overs = over(slu, sl, returnList=TRUE)
    ## overlay is true if end points overlay, so filter them out:
    overs = lapply(1:length(overs), function(islu){
        Filter(function(isl){
            islines(sl[isl,],slu[islu,])
        }, overs[[islu]])
    })
    ## now aggregate the required attribute using fun():
    aggs = sapply(overs, function(os){fun(sldf[[attr]][os])})

    ## make a SLDF with the named attribute:
    sldf = SpatialLinesDataFrame(slu, data.frame(Z=aggs))
    names(sldf)=attr
    sldf
}

lineLabels <- function(sldf, attr){
    text(coordinates(gCentroid(sldf,byid=TRUE)),labels=sldf[[attr]])
}

Usage:

> r <- readOGR("airport.geojson", layer = "OGRGeoJSON")
> ag = aggit(r,"All")
> plot(ag)
> lineLabels(ag,"All")

ag is now a spatial lines data frame with the aggregated variable. The plot lets you check it all adds up. Because your routes are a bit like a river network, you can see how the "flows" add up at the junctions:

enter image description here

Seems right to me...

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.