Extracting GPS coordinates from SF network build from OSM Data, provides those in a wrong order

[Andreas4242]

I encountered a strange issue while working with pathfinding within SF Network. My goal is to find the shortest path and display the resulting GPS trace. However, while extracting the GPS coordinates, I sometimes get them in the wrong order.. meaning backward.

This is how I get the Path and the coordinates.

path <- shortest_paths(my_sfn, from = start_point, to = dest_point,mode = "out", weights = weight, output = "both")
path<-as.integer(path$epath[[1]])
coordinates_of_the_path<-as.data.table(st_coordinates(get.edge.attribute(graph = my_sfn, name = "geometry", index = path)))

This creates the following picture (dont mind the blue background):

Or this:

The data Looks similar to this:

Although the position of the GPS points is correct, the connection between the id blocks (from the path) is sometimes wrong i.e backwards. For example the data with ID 5102 is correct from position but needs to connect to the data with ID 5105 not with the point NR 21 but with the point NR 23 ie the whole "sub data frame" should be reversed. (just for example to explain the issue).

Is it something I do wrong or is it something I missed or a bug? I apologies beforehand if its something I missed!

SF Network is created by using OSM data like this:

my_sfn <- st_as_sf(sf_data)
my_sfn <- as_sfnetwork(sf_data[, c("osm_id", "highway","bridge","tunnel","ref","surface")], directed = FALSE)

 my_sfn <- convert(my_sfn, to_spatial_subdivision, .clean = TRUE)
 my_sfn <- my_sfn %E>% mutate(weight = edge_length())

 my_sfn <- my_sfn %>% activate("edges") %>%
   filter(!edge_is_multiple()) %>%
   filter(!edge_is_loop())

To provide another example for a better understanding:

I get a path with edges 5105,5102 ,5101, 5162,5163 Those transferred to gps points create this table:

now if i calculate the distances between points in a row Nr1 and Nr8 and point Nr7 and Nr8:

  #point 1  to point  8 = 35 m
  geosphere::distGeo(c(11.8534718 ,48.1514218),c(11.8539547 ,48.1514243))
  #point 7 to  point  8 = 217 m
  geosphere::distGeo(c(11.8535345, 48.1494923),c(11.8539547 ,48.1514243))

That means either all the points with ref St 2081 or EBE 5 need to be reversed.

[agila5]

Hi @andreas-ullrich! Can you please paste the output of dput(sf_data) or upload sf_data after saving the data as a GeoPackage?

[Andreas4242]

Hi @agila5 here you go. I saved my_sfn and path that you can try as rds objects. my_sfn_and_path.rds.zip

[agila5]

Can you please provide a full and reproducible example showing your problem? I cannot replicate the issue that you mentioned. For example: What is start_point? dest_point? weight? How did you create those plots?

[Andreas4242]

Oh sorry. Sure here you go. This will generate the gpx file. andrea_example.zip

[agila5]

Hi @andreas-ullrich! The problem is that, for some reason, some segments are defined in the opposite order with respect to the others. For example

# First, load packages
library(sfnetworks)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(igraph, warn.conflicts = FALSE)
library(tidygraph, warn.conflicts = FALSE)

# and data
load("data.RData")

# Define two points for a toy example
start_point <- 5329
dest_point<- 4447

# Compute the shortest path among them and extract the edges, then subset only
# the third and fourth ones
path <- shortest_paths(my_sfn, start_point, dest_point, weights = weight, output = "both")
path <- as.integer(path$epath[[1]])
subset_segments <- get.edge.attribute(graph = my_sfn, name = "geometry", index = path)
subset_segments <- subset_segments[3:4]

# If we extract their coordinates, we can see that, for some reason, they are
# defined in the opposite order (see the first an last row)
st_coordinates(subset_segments)
#>             X        Y L1
#> [1,] 11.83332 48.12696  1
#> [2,] 11.83309 48.12680  1
#> [3,] 11.83384 48.12715  2
#> [4,] 11.83358 48.12706  2
#> [5,] 11.83332 48.12696  2

# Therefore
plot(subset_segments)

# while, if we plot the raw coordinates and connect them... 
plot(st_coordinates(subset_segments), type = "l")

^{Created on 2023-04-08 with reprex v2.0.2}

Unfortunately, I cannot help you more than this unless you add the code you used to produce my_sfn.

[Andreas4242]

Hi Andrea,

this is correct. I got a remedy for this problem suggested by someone on stack overflow but it dies not solve the problem about why this is happening.

https://stackoverflow.com/questions/75938726/reverse-elements-of-list-of-data-frames-based-on-the-distance-of-points-within-t

I wrote to you in the previous mail how I create sf sf netowrk:

OSM data is in a postgres database. Here Im connecting to it and getting an ST object and create sf network. Im not doing anything to it to make this "wrong" result.

sf_data = st_read(con1, query = query) dbDisconnect(con1)

my_sfn <- st_as_sf(sf_data)

my_sfn <- as_sfnetwork(sf_data[, c("osm_id", "highway","bridge","tunnel","ref","surface")], directed = FALSE)

my_sfn <- convert(my_sfn, to_spatial_subdivision, .clean = TRUE) my_sfn <- my_sfn %E>% mutate(weight = edge_length())

my_sfn <- my_sfn %>% activate("edges") %>% filter(!edge_is_multiple()) %>% filter(!edge_is_loop())

What is strange: if I use

lines_merged <- st_cast(st_line_merge( st_union(st_cast(get.edge.attribute(graph = my_sfn, name = "geometry", index = path ), "MULTILINESTRING"))), "LINESTRING")

This connects the points in the correct direction.

Best regards Andreas

Am 08.04.2023 um 09:10 schrieb Andrea Gilardi @.***>:

Hi @andreas-ullrich! The problem is that, for some reason, some segments are defined in the opposite order with respect to the others. For example

First, load packages

library(sfnetworks) library(sf)

> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(igraph, warn.conflicts = FALSE) library(tidygraph, warn.conflicts = FALSE)

and data

load("data.RData")

Define two points for a toy example

start_point <- 5329 dest_point<- 4447

Compute the shortest path among them and extract the edges, then subset only

the third and fourth ones

path <- shortest_paths(my_sfn, start_point, dest_point, weights = weight, output = "both") path <- as.integer(path$epath[[1]]) subset_segments <- get.edge.attribute(graph = my_sfn, name = "geometry", index = path) subset_segments <- subset_segments[3:4]

If we extract their coordinates, we can see that, for some reason, they are

defined in the opposite order (see the first an last row)

st_coordinates(subset_segments)

> X Y L1

> [1,] 11.83332 48.12696 1

> [2,] 11.83309 48.12680 1

> [3,] 11.83384 48.12715 2

> [4,] 11.83358 48.12706 2

> [5,] 11.83332 48.12696 2

Therefore

plot(subset_segments)

while, if we plot the raw coordinates and connect them...

plot(st_coordinates(subset_segments), type = "l") Created on 2023-04-08 with reprex v2.0.2 Unfortunately, I cannot help you more than this unless you add the code you used to produce my_sfn. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>

[agila5]

OSM data is in a postgres database.

And that is the OSM data for the whole planet? Moreover, can you share also which query do you use?

[Andreas4242]

Its this file:

https://download.geofabrik.de/europe/dach-latest.osm.pbf

Sure the query is like this:

query<-paste( "SELECT osm_id,highway,bridge,tunnel,ref,surface,ST_Transform(way,4326) AS geometry FROM planet_osm_line WHERE highway IN ('primary','secondary','tertiary','primary_link','secondary_link','tertiary_link') AND planet_osm_line.way && ST_Transform( ST_MakeEnvelope(",bbox_dimensions[[1]],",",bbox_dimensions[[2]],",",bbox_dimensions[[3]],",",bbox_dimensions[[4]],",4326),3857 );") sf_data = st_read(con1, query = query)

with the bbox dimensions for a small part of the file:

xmin ymin xmax ymax 

11.71541 47.77093 12.32288 48.17883

Am 08.04.2023 um 09:26 schrieb Andrea Gilardi @.***>:

OSM data is in a postgres database. And that is the OSM data for the whole planet? Moreover, can you share also which query do you use? — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>