[R-sig-Geo] Generating Random Planar Graph with 1m Edges

[Krenz, Kimon-Vincent]

Dear All,

I started a week ago to use R to solve a problem I am currently facing in my PhD.
Apologies in advance for the long-winded explanation of my problem.

I am trying to generate a random planar graph with approximately 1 million edges, where:

A) nodes (points) feature spatial coordinates
B) the network has a given boundary

C) edges (lines) are created if two points fall within a given distance (e.g. 100 - 1000 metres)
D) degree (connectivity) ranges between given k max and min (e.g. k ≤ 5)
E) edges do not intersect

This will result in something one might want to compare to a random street network.

I am following a method proposed here: Masucci, A. P., Smith, D., Crooks, A., & Batty, M. (2009). Random planar graphs and the London street network. European Physical Journal B, 71(2), 259–271. http://doi.org/10.1140/epjb/e2009-00290-4) link to paper: https://goo.gl/6XWW8P

Masucci et al.: "We first introduce a random model for a static planar graph. ... To build an ERPG we start with a Poisson distribution of N points in a plane and we choose a distance r. To build the first segment, we randomly pick up two points of the distribution that have a distance less then r and we connect them. Then we continue to randomly pick up pairs of points P and Q in the given points distribution that have a distance less then r. If the segment PQ does not intersect any other line of the graph, we add it to the graph. The process ends when we add the desired number of edges E.”

I hence, started with generating random points using the Poisson distribution in a given spatial box (A and A):

require(spatstat)
require(maps)
library(sp)
w <- as.owin(list(xrange=c(32275175,32475175),yrange=c(5611910,5811910)))
Y <- rpoispp(0.00001, win=w)
Ydata <- data.frame(Y)
list(Y)

That seems to be quite straightforward in R.
I then followed the proposed method and wrote a simple loop, that selects two points from Ydata
based on a random sample and adds the projection of the coordinate system.

#399717 is = N from the rpoispp

repeat {
  l1 <- sample(1:399717, 2, replace=F)

 Ydata.sp<-Ydata[l1, c('x','y')]
  coordinates(Ydata.sp)<-~x+y
  crs.geo <- CRS("+init=epsg:4647 +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=32500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
  proj4string(Ydata.sp) <- crs.geo

  L = SpatialLines(list(Lines(list(Line(coordinates(Ydata.sp))),"X")))
  if (SpatialLinesLengths(L)<=4000){
    break
  }
}
plot(L)
SpatialLinesLengths(L)

This fulfils the requirement of C) and I have not yet reached the point to deal with D) and E).
My problem here is that the process takes way too long for a single pair (approx 20 seconds, resulting in 250 days computational time).

Is there a quicker way solve this in a more efficient manner?

I have thought about selecting the first point by random and the second one randomly based on an evaluation
of all points within a given radius to the first point, rather than two complete random points.

This would at least cut down the test of several thousand meaningless combinations. However, I couldn't find a way to do this.
Another option might be gDistance from the (rgeos), but with 400k points, the result seems to be not computable.

I am also happy for any suggestions regarding requirements D and E, or help with the task in general.

Best,
Kimon

Kimon Krenz

PhD Researcher
MSc SDAC Course Coordinator<http://www.bartlett.ucl.ac.uk/space-syntax/programmes/mres-msc/msc-spatial-design>
Space Syntax Laboratory<http://www.bartlett.ucl.ac.uk/space-syntax>

mail.       ucftkr3 at ucl.ac.uk<mailto:ucftkr3 at ucl.ac.uk>
phone.    0044 7784 329089
web.       www.kimonkrenz.com<http://www.kimonkrenz.com/>

Bartlett School of Architecture<http://www.bartlett.ucl.ac.uk/>
UCL
140 Hampstead Road
London     NW1 2BX     UK


	[[alternative HTML version deleted]]