Hi,

Following Barry's code and using spatstat::project2segment, I implemented
the integration of extra nodes to the graph

I tried to give different attributes to vertexes depending on whether they
were on the lines or whether they're extra nodes, in order to facilitate
queries but I don't find my code very clean.

It's :

  E(graph)$weight=lengths
  nVtx <- length(V(graph))

  # TYPE ATTRIBUTE
  V(graph)$type <- "Line"
  V(graph)[(nVtx-nrow(pts)+1):nVtx]$type <- "Point"


Would you see any improvement to my code?

Here is the "new" function based on Barry's one

buildTopo <- function(lines, pts){
  require(rgeos)
  require(igraph)
  require(spatstat)

  # (1) LINEAR NETWORK
  g = gIntersection(f,f)
  edges1 = do.call(rbind,lapply(g@lines
[[1]]@Lines,function(ls){as.vector(t(ls@coords))}))

  # (2) POINT LOCATIONS
  froms2 <- as(pts,"ppp")
  tos2 <- project2segment(froms2, as(g, "psp"))$Xproj
  edges2 <- cbind(froms2$x, froms2$y, tos2$x, tos2$y)

  # MERGING BOTH ((1)+(2))
  edges <- rbind(edges1, edges2)

  # LENGTH CALCULATION
  lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)

  # GRAPH
  froms = paste(edges[,1],edges[,2])
  tos = paste(edges[,3],edges[,4])

  graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
  E(graph)$weight=lengths
  nVtx <- length(V(graph))

  # TYPE ATTRIBUTE
  V(graph)$type <- "Line"
  V(graph)[(nVtx-nrow(pts)+1):nVtx]$type <- "Point"

  return(graph)
}

m = readOGR(".","Line")
pts <- readOGR(".", "Point")
gg = buildTopo(m, pts)
E(gg)[get.shortest.paths(gg,V(gg)[1],V(gg)[14])[[1]]-1]

I'm glad to have, now, a totally R topology function!

Mathieu

2012/8/31 Mathieu Rajerison <mathieu.rajerison@gmail.com>

> Wow, I'm really impressed by your geohacking skills. Thanks a lot.
>
> Yes, there is a GRASS::v.clean function that cuts lines.
> Also, GRASS already includes allocation algorithms.
>
> We could also imagine performing shortest path calculations for extra
> vertexes that aren't initally on the graph: shop locations, particular
> events, houses...
>
> For this, spatstat::project2segment might be the appropriate function.
>
> Based on your initial code, I think we could write a function that would
> combine a linear network with extra vertexes into a single graph. The
> vertex labels could discriminate the intersection nodes from the thematic
> ones.
>
>
>
> 2012/8/31 Barry Rowlingson <b.rowlingson@lancaster.ac.uk>
>
>> On Fri, Aug 31, 2012 at 11:58 AM, Rowlingson, Barry
>> <b.rowlingson@lancaster.ac.uk> wrote:
>>
>> >  If I didn't have to be on a train in three hours I'd code this up...
>>
>>  Oh who am I kidding:
>>
>> buildTopo <- function(lines){
>>   require(rgeos)
>>   require(igraph)
>>
>>   g = gIntersection(lines,lines)
>>   edges = do.call(rbind,lapply(g@lines
>> [[1]]@Lines,function(ls){as.vector(t(ls@coords))}))
>>   lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)
>>
>>   froms = paste(edges[,1],edges[,2])
>>   tos = paste(edges[,3],edges[,4])
>>
>>   graph = graph.edgelist(cbind(froms,tos),directed=FALSE)
>>   E(graph)$weight=lengths
>>   return(graph)
>>
>> }
>>
>>  > m = readOGR(".","Line")
>>  > gg = buildTopo(m)
>>
>> compute shortest distance between two vertices (note the -1 is because
>> I'm still on old igraph with 0-indexing):
>>
>>  > E(gg)[get.shortest.paths(gg,V(gg)[1],V(gg)[14])[[1]]-1]
>>
>> Edge sequence:
>>
>> [0]  0.81711193 47.57711634 -- -0.28546154 48.4287593
>> [1]  2.2151139 46.49728033  -- 0.81711193 47.57711634
>> [9]  1.57160852 45.36348513 -- 2.2151139 46.49728033
>> [10] 1.39973834 45.06066625 -- 1.57160852 45.36348513
>> [11] 1.29755246 44.88062446 -- 1.39973834 45.06066625
>> [12] 0.91544563 43.67971729 -- 1.29755246 44.88062446
>> [13] 0.88815229 43.43407719 -- 0.91544563 43.67971729
>>
>> I think this could be improved by adding the coordinates to the nodes.
>> But anyway, 90% done.
>>
>> Barry
>>
>
>

	[[alternative HTML version deleted]]

