0

1. Problem

I have the coordinates for over 5000 points (cities) that lie over a road shape file, as partially illustrated below:

enter image description here

I'm working on R, and I need to find the shortest distance between all pairs of points along the roads.

2. What I've already tried, without success:

  • stplanr package: had a problem converting the full road shape file to a SpatialLinesNetwork object ("non finite line lengths");
  • osrm package: can't compute all distance pairs (over 30 million) due to server restrictions (up to 10k pairs).

3. Sample data

I provide a sample data here: https://github.com/cassianord/example_data/blob/main/example_data.7z?raw=true

8
  • sl = rgdal::readOGR(".","roads_ex") ; sln = SpatialLinesNetwork(sl) works for me with no warnings or errors about non finite line lengths. Version 0.8.1 of stplanr.
    – Spacedman
    Feb 21, 2021 at 16:58
  • Yes, in the sample data that might work, but with the full shape it doesn't =/
    – cassianord
    Feb 21, 2021 at 17:10
  • Can you supply a subset where it doesn't work? If we knew where your "non-finite line length" lines were.... Is it simply lines that are length zero and could be dropped anyway?
    – Spacedman
    Feb 21, 2021 at 17:53
  • 1
    Sure, sorry for not doing it before. Here follows a larger sample where the issue happens: drive.google.com/file/d/1Nt0WWSb0pOgRcS4pGOHrR1S3sO0lkgSF/…
    – cassianord
    Feb 21, 2021 at 18:47
  • Hi @cassianord. What's the CRS of the data? 4326?
    – agila
    Feb 21, 2021 at 20:25

1 Answer 1

1

At least one row of your roads is causing this error when SpatialLinesNetwork tries to get its length:

> SpatialLinesLengths(l[65403,])
Error in FUN(X[[i]], ...) : non-finite line lengths

This row appears to have zero length, which is messing up the length calculation when it is using lat-long coordinates:

> l[65403,]
class       : SpatialLinesDataFrame 
features    : 1 
extent      : -41.98302, -41.98302, -15.25437, -15.25437  (xmin, xmax, ymin, ymax)

I suspect a divide-by-zero somewhere when working out the distance on a sphere by lat-long.

The adjacent segments are fine:

> SpatialLinesLengths(l[65402,])
[1] 1.813751
> SpatialLinesLengths(l[65404,])
[1] 1.802447

You can trick it into not using lat-long coordinates - there is an option but I think having a CRS on the lines overrides this anyway. So set the CRS to nothing:

> proj4string(l) = ""

And then I can build a spatial lines network without errors:

> ln = SpatialLinesNetwork(l)
> 

But that will give you a spatial network based on euclidean lat-long distance which might not be very valid.

Or you could project to an appropriate cartesian grid system for your data and use that.

Alternatively you could try an eliminate these rows from your data, but in some cases it seems to be caused by zero-length segments within road segments, which is hard to track down.

1
  • Thanks, @Spacedman, it solved the issue with the SpatialLinesNetwork. Now I'm facing other problems in attributing nodes to the network... also, I think that, after all, maybe stplanr will not help me to find the distances between all pairs. Nevertheless, I have upvoted your answer and accepted it. Thanks again!
    – cassianord
    Feb 22, 2021 at 17:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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