[R-sig-Geo] Generating Random Planar Graph with 1m Edges
Krenz, Kimon-Vincent
kimon-vincent.krenz.12 at ucl.ac.uk
Tue Oct 18 20:47:52 CEST 2016
Dear Roger,
Apologies for the late reply.
Your suggestion was a life saver, the computation time by which RANN produces nearest neighbour indexes and distances is incredible. Even with more than 1.5million points it took only a minute to compute 30 neighbours and their distance for each pair.
The only problem I faced with the RANN approach is that one arrives with the nearest points first. Meaning that if a point is surrounded by lets say 30 other points within a distance of 500 meter, one will never arrive with a pair that is in distance of 5000 meter. To overcome this one can only increase the number of k neighbours, which in my case would mean minimum k=100, leading to 100million combinations of which I then sampled randomly the amount necessary to proceed.
Still the best and fastest option. Thanks for your help!
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
On 06 Oct 2016, at 13:32, Roger Bivand <Roger.Bivand at nhh.no<mailto:Roger.Bivand at nhh.no>> wrote:
Just a partial suggestion: the RANN package will let you index the points byk-neighbourness, avoiding some of the burden of searching for points within your 100-1000m band.
Hope this helps,
Roger
On Wed, 5 Oct 2016, Krenz, Kimon-Vincent wrote:
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><mailto:ucftkr3 at ucl.ac.uk>
phone. 0044 7784 329089
web. www.kimonkrenz.com<http://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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org<mailto:R-sig-Geo at r-project.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no<mailto:Roger.Bivand at nhh.no>
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list