1

I'm looking for a way to get the xy-coordinates of all the intersections within one SpatialLines object or SpatialLinesDataFrame. I have found the function gIntersect of rgeos but that only looks at the intersection between two datasets. Since I am working with a dataset of over half a million lines it would take too much time to make a separate file of every line and check whether any line intersects with another. In ArcMap there is the Intersect function that is able to do it in a couple of seconds and I was wondering whether there was also such a function in R. Thanks!

jbaums
  • 27,115
  • 5
  • 79
  • 119
jonas van duijvenbode
  • 61
  • 10

1 Answers1

1

If you convert your SpatialLines object into a psp object from spatstat you can use the spatstat function selfcrossing.psp. However, I'm not sure how it will cope with half a million lines since the number of crossings potentially could be enormous. The code below generates a random segment pattern and finds the self crossings.

BEWARE that this code potentially can take up a lot of memory and kill R, so try with progressively increasing examples before processing a half million lines. The code below used quite a bit of memory on my 5 year old laptop and took 5 seconds to run.

set.seed(42)
N <- 1e4
x <- psp(runif(N), runif(N), runif(N), runif(N), owin(), check=FALSE)
y <- selfcrossing.psp(x)
Ege Rubak
  • 4,347
  • 1
  • 10
  • 18
  • Many thanks Ege. I ran your script but indeed as you said it takes a lot of processing power. when I tried 100.000 lines my 8gb memory got filled so half a million is not going to work. I know that the function that both ArcMap and QGIS use can run the same dataset in a couple of seconds but I would really like to integrate it in R – jonas van duijvenbode Oct 29 '14 at 09:23
  • Maybe you can split the problem into several smaller ones by splitting the dataset into e.g. 50 line patterns with 10,000 lines each and the make a loop to find the crossings between all of them (there is a function called `crossing.psp` for crossings between patterns). It probably won't be very fast, but at least you will have enough memory, and get the job done in R. – Ege Rubak Oct 29 '14 at 12:39
  • I got around 100.000 points out of the crossings in ArcMap. I am now trying to incorporate R into QGIS, which also seems like an impossible task to do but that will at least make the process chain shorter. – jonas van duijvenbode Oct 29 '14 at 17:24
  • OK. The algorithm in `spatstat` is quite generic C code which doesn't take advantage of the fact that few lines intersect. It is probably only competitive in terms of speed and memory when the line pattern is dense with a large proportion of the lines crossing each other. We might implement a more modular version, but it probably wont be in the very near future. Sorry about that. – Ege Rubak Oct 30 '14 at 19:50
  • I tried the selfcrossing.psp function and since there were not as many intersects as in your function the method worked okay. It seems like the long execution time is caused by the great number of intersects as a result of your sample script rather than the actual number of line elements. thank you very much for the tips. – jonas van duijvenbode Oct 30 '14 at 21:20
  • OK. That is good news. How do the computation times compare between ArcMap and `spatstat`? If the answer solves you problem you should consider accepting it, so other people can see the problem is solved. – Ege Rubak Oct 31 '14 at 07:10