Reconstructing way–node incidence from `oe_get()` output — is there a better way?

[gustavo-tm]

Hi — first, thanks for the package, it's been really useful.

I'm working on a road safety project where I need the way-node incidence relation.

As far as I can tell this relation isn't available directly from oe_get(), since the OSM driver reconstructs geometries and doesn't keep node references. So I've been rebuilding it by exploding the linestrings into vertices and matching those coordinates against the points layer.

The part I'm unsure about is that this only works if I round the coordinates to 7 decimal places before joining — on raw doubles the match fails for a lot of nodes. So in practice I'm snapping and rebuilding the geometries before the join, which feels like something I probably shouldn't need to do. Rounding to 7 dp seems to line up with OSM's native storage precision, so maybe that's the right thing, but I'm not confident.

Is coordinate matching the intended approach here, or is there a supported way to get node references that I've missed?

[agila5]

Hi @gustavo-tm! Unfortunately, this is a known limitation of the approach employed in our package and, as far as I know, it's not possible to reconstruct the node/way relationship from GDAL OSM drivers (unless you are using a spatial match operation as you describe).

Did you already check the data returned by osmdata? As far as I know, they are derived from the Overpass API and they preserve the OSM IDs that can be used to match nodes and ways. For example:

Load packages

library(sf)
#> Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

Get OSM roads for a small abandoned town in northern Italy

cn <- opq("Consonno, Italy", timeout = 75) |> add_osm_feature(key = "highway") |> osmdata_sf()

Data are structured as follows.

cn
#> Object of class 'osmdata' with:
#>                  $bbox : 45.7642766,9.3735149,45.8042766,9.4135149
#>         $overpass_call : The call submitted to the overpass API
#>                  $meta : metadata including timestamp and version numbers
#>            $osm_points : 'sf' Simple Features Collection with 7437 points
#>             $osm_lines : 'sf' Simple Features Collection with 828 linestrings
#>          $osm_polygons : 'sf' Simple Features Collection with 5 polygons
#>        $osm_multilines : NULL
#>     $osm_multipolygons : NULL

We can extract the ways/lines:

head(cn$osm_lines)
#> Simple feature collection with 6 features and 76 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 9.400817 ymin: 45.80063 xmax: 9.413871 ymax: 45.80382
#> Geodetic CRS:  WGS 84
#>            osm_id                          name abandoned abandoned:highway
#> 24179481 24179481         Via per Calolziocorte      <NA>              <NA>
#> 24179638 24179638 Via Professore Mario Redaelli      <NA>              <NA>
#> 46940506 46940506 Via Professore Mario Redaelli      <NA>              <NA>
#> 52047343 52047343         Via Don Carlo Gnocchi      <NA>              <NA>
#> 54369532 54369532              Strada del Ronco      <NA>              <NA>
#> 54373564 54373564              Via Valcavellino      <NA>              <NA>
#> (Omitted...)
#>          width                       geometry
#> 24179481  <NA> LINESTRING (9.41281 45.8020...
#> 24179638  <NA> LINESTRING (9.41281 45.8020...
#> 46940506  <NA> LINESTRING (9.408931 45.803...
#> 52047343  <NA> LINESTRING (9.413447 45.800...
#> 54369532  <NA> LINESTRING (9.400994 45.802...
#> 54373564  <NA> LINESTRING (9.402688 45.802...

and may plot them:

plot(st_geometry(cn$osm_lines))

We can see that, for example, the first way has ID 24179481. We can extract all nodes that form such way as follows

osm_points(cn, "24179481")
#> Simple feature collection with 7 features and 56 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 9.41281 ymin: 45.80202 xmax: 9.413871 ymax: 45.80247
#> Geodetic CRS:  WGS 84
#>                  osm_id name abandoned access advertising amenity barrier bench
#> 261853435     261853435 <NA>      <NA>   <NA>        <NA>    <NA>    <NA>  <NA>
#> 444538884     444538884 <NA>      <NA>   <NA>        <NA>    <NA>    <NA>  <NA>
#> 1981147433   1981147433 <NA>      <NA>   <NA>        <NA>    <NA>    <NA>  <NA>
#> 1981147464   1981147464 <NA>      <NA>   <NA>        <NA>    <NA>    <NA>  <NA>
#> 5520942304   5520942304 <NA>      <NA>   <NA>        <NA>    <NA>    <NA>  <NA>
#> 12169517036 12169517036 <NA>      <NA>   <NA>        <NA>    <NA>    <NA>  <NA>
#> 14110607089 14110607089 <NA>      <NA>   <NA>        <NA>    <NA>    <NA>  <NA>
#> (Omitted...)
#>             wheelchair wikimedia_commons                  geometry
#> 261853435         <NA>              <NA>  POINT (9.41281 45.80202)
#> 444538884         <NA>              <NA> POINT (9.413871 45.80247)
#> 1981147433        <NA>              <NA> POINT (9.413417 45.80229)
#> 1981147464        <NA>              <NA> POINT (9.413747 45.80242)
#> 5520942304        <NA>              <NA> POINT (9.413835 45.80245)
#> 12169517036       <NA>              <NA> POINT (9.412873 45.80205)
#> 14110607089       <NA>              <NA> POINT (9.413294 45.80223)

Graphical check

plot(st_geometry(cn$osm_lines)[1], reset = FALSE)
plot(osm_points(cn, "24179481")[, "geometry"], pch = 16, col = "red", add = TRUE)

^{Created on 2026-08-26 with reprex v2.1.1.9000}

[gustavo-tm]

Hello, @agila5!

Thanks for the fast response!

Unfortunately I need the geometries from São Paulo, Brazil, which is a huge city thus I can't get all the data from osmdata. Currently I have the following workflow for building incidence. Hopefully this procedure is robust for correct matching.

Would something along these lines be in scope for the package either as a vignette on recovering the node/way relation, or as a helper that returns the incidence table (and optionally node degree) directly from an extract?

tidy_osm_ways <- function(municipality_geometry){
  
  if (!file.exists("data/osm_ways.gpkg")){
    # ways database
    ways <- "Sao Paulo" |> 
      
      # Get OSM data
      oe_get(
        place = _, 
        layer = "lines",
        provider = "openstreetmap_fr",
        extra_tags = c("maxspeed", "oneway", "lanes", "alt_name", "ref", "junction"),
        boundary = municipality_geometry, 
        boundary_type = "spat",
        force_vectortranslate = TRUE) |> 
      
      # Select and rename columns
      select(
        id_way = osm_id,
        name,
        name_alt = alt_name,
        name_ref = ref,
        type = highway,
        lanes,
        speed_limit = maxspeed,
        oneway,
        junction,
        geom = geometry) |> 
      
      # Filter only roads 
      # discard footways, aeroways, waterways, railways, etc.
      filter(type %in% c(
        "motorway", "trunk", "primary", "secondary", "tertiary", "residential", 
        # "unclassified", "service",
        "motorway_link", "trunk_link", "primary_link", "secondary_link")) |> 
      
      # Only the municipality region
      st_join(municipality_geometry, left = FALSE) |> 
      
      # Correcting column formats
      mutate(lanes = as.numeric(lanes),
             speed_limit = as.numeric(speed_limit),
             oneway = case_when(oneway == "yes" ~ TRUE,
                                oneway == "no" ~ FALSE,
                                .default = NA)) |> 
      
      # Identifying roundabouts
      mutate(roundabout = junction %in% c("roundabout", "circular")) |> 
      select(-junction)
    
    # Save data for reproducibility; OSM servers overwrite old data
    ways |> st_write("data/osm_ways.gpkg", append = FALSE)
  } 
  else {
    ways <- st_read("data/osm_ways.gpkg")
  }
  return(ways)
}

tidy_osm_points <- function(municipality_geometry){
  
  if (!file.exists("data/osm_points.gpkg")){
    
    # Points database
    points <- "Sao Paulo" |> 
      
      # Get OSM data
      oe_get(
        place = _, 
        layer = "points",
        provider = "openstreetmap_fr",
        extra_tags = c("crossing", "traffic_signals", "traffic_calming"),
        boundary = municipality_geometry, 
        boundary_type = "spat",
        force_vectortranslate = TRUE) |> 
      
      # Create dummies for information of interest
      mutate(
        traffic_signal = highway %in% "traffic_signals" | crossing %in% "traffic_signals",
        crossing = highway %in% "crossing" | (!is.na(crossing) & !crossing %in% "no"),
        traffic_calming = !is.na(traffic_calming) & !traffic_calming %in% "no",
        roundabout = highway %in% "mini_roundabout",
        speed_camera = highway %in% "speed_camera"
      ) |> 
      
      # Select wanted columns
      select(id_point = osm_id, 
             traffic_signal, crossing, traffic_calming, speed_camera, roundabout,
             geom = geometry) |> 
      
      # Filter only points with information of interest
      filter(traffic_signal + crossing + traffic_calming + speed_camera + roundabout > 0) |> 
      
      # Only the municipality region
      st_join(municipality_geometry, left = FALSE)
    
    # Save data for reproducibility; OSM servers overwrite old data
    points |> st_write("data/osm_points.gpkg", append = FALSE)
  } 
  else {
    points <- st_read("data/osm_points.gpkg")
  }
  return(points)
}

tidy_osm_nodes_preprocess <- function(ways){
  
  # Turn the linestring into multiple points, while keeping the order (position)
  ways_to_nodes <- ways |> 
    st_coordinates() |> 
    as.data.frame() |> 
    mutate(id_way = ways$id_way[L1]) |>
    mutate(.by = L1,
           position = row_number()) |> 
    
    # Rounding coordinates is important for exact node join
    # Original data is in 7 digit precision
    mutate(across(c(X, Y), ~ round(.x, 7)))
  
  # Create a unique id for each node
  nodes_ids <- ways_to_nodes |> 
    distinct(X, Y) |> 
    mutate(id_node = row_number())
  
  return(list(ways_to_nodes, nodes_ids))
}

tidy_osm_incidence <- function(nodes_preprocess){
  
  ways_to_nodes <- nodes_preprocess[[1]]
  nodes_ids <- nodes_preprocess[[2]]
  
  # Each node now has an unique ID
  incidence <- ways_to_nodes |> 
    left_join(nodes_ids) |> 
    select(id_way, id_node, position)
  
  return(incidence)
  
}

tidy_osm_nodes <- function(points, incidence, nodes_preprocess){
  
  nodes_ids <- nodes_preprocess[[2]]
  
  # Extract coordinates for exact join with nodes
  points_prejoin <- points |> 
    st_coordinates() |> 
    as.data.frame() |> 
    mutate(id_point = points$id_point) |> 
    
    # Rounding coordinates is important for exact node join
    # Original data is in 7 digit precision
    mutate(across(c(X, Y), ~ round(.x, 7)))
  
  # Calculate node degree
  # A node is terminal if it is the first or last of a linestring
  degree <- incidence |> 
    mutate(.by = id_way,
           terminal = position == 1 | position == max(position)) |>
    summarise(.by = id_node,
              degree = sum(if_else(terminal, 1L, 2L)),
              n_ways = n_distinct(id_way))
  
  # Combine unique nodes with 
  # (i) osm id for them if they exist i.e. the node has information
  # (ii) node degree
  nodes <- nodes_ids |> 
    left_join(points_prejoin) |> 
    left_join(degree) |> 
    select(id_node, id_point, degree, X, Y)
  
  return(nodes)
}

[agila5]

Hi @gustavo-tm! Ok, I see you point.

For the time being, I believe that we can showcase these functions into the new "tips and tricks" vignette (https://github.com/ropensci/osmextract/blob/master/vignettes/tips-and-tricks.Rmd) using a small OSM extract. Would you like to start working on the first draft? Please present a self-contained example using an OSM extract which is a small as possible (possibly the ones already included in the package so the functions do not require internet connection).

Moreover, i briefly checked the OSM wiki (https://wiki.openstreetmap.org/wiki/Node) and it looks like 7 decimals of precision are ok.

[gustavo-tm]

@agila5, I might have found a solution for the problem.

I assume that for most users correctly combining points and lines is enough. In that case, it seems like adding sf::st_set_precision(1e7) to the geometries fixes the line-point join issue. However, for building incidence (which is my case), the precision is not extended to sf::st_coordinates sadly.

This probably needs further investigation, but it looks like adding sf::st_set_precision(1e7) to imported data should be a package default on osmextract::oe_get and others. On the other hand, this entire discussion about decimal precision could be somewhere to help people who want to recreate node relations, as exploding lines into nodes with sf::st_coordinates does cause this issue and manual snapping is required.

Image of naive join failing:

library(osmextract)
library(sf)
library(dplyr)
library(ggplot2)

lines  <- oe_get("ITS Leeds", layer = "lines",  quiet = TRUE) |> select(osm_id)
points <- oe_get("ITS Leeds", layer = "points", quiet = TRUE) |> select(osm_id)

nrow(lines)   #> 189
nrow(points)  #> 186


# The problem -------------------------------------------------------------

joined_geoms <- lines |>
  st_join(points)

orphan_points <- points |>
  anti_join(joined_geoms |> st_drop_geometry() |> select(osm_id = osm_id.y),
            by = "osm_id")

n_distinct(na.omit(joined_geoms$osm_id.y))  #> 14  of 186 points
nrow(orphan_points)                         #> 172

# 172 orphans out of 186. Nothing errors, nothing warns.


# The reason --------------------------------------------------------------

vertices <- lines |>
  st_coordinates() |>
  as.data.frame() |>
  distinct(X, Y) |>
  st_as_sf(coords = c("X", "Y"), crs = 4326)

d <- st_distance(orphan_points,
                 vertices[st_nearest_feature(orphan_points, vertices), ],
                 by_element = TRUE) |> as.numeric()

id <- orphan_points$osm_id[d < 1e-6][1]   #> "21069418"

p <- points |> filter(osm_id == id) |> st_coordinates()
v <- lines |> st_coordinates() |> as.data.frame() |>
  filter(round(X, 7) == round(p[1], 7), round(Y, 7) == round(p[2], 7)) |> slice(1)

sprintf("%.17g", c(p[1], v$X))   #> "-1.5525082000000001"  "-1.5525081999999999"
sprintf("%.17g", c(p[2], v$Y))   #> "53.808604800000005"   "53.808604799999998"



# The fix -----------------------------------------------------------------

joined_snapped <- st_set_precision(lines, 1e7) |>
  st_join(st_set_precision(points, 1e7))

n_distinct(na.omit(joined_snapped$osm_id.y))  #> 40



# Figure ------------------------------------------------------------------

m0 <- unique(na.omit(joined_geoms$osm_id.y))
m7 <- unique(na.omit(joined_snapped$osm_id.y))

points <- points |> mutate(status = case_when(
  osm_id %in% m0 ~ "matched by st_join (14)",
  osm_id %in% m7 ~ "orphan, but lying ON a line (26)",
  .default       = "orphan, genuinely off-network (146)"))

ggplot() +
  geom_sf(data = lines, colour = "grey72", linewidth = 0.35) +
  geom_sf(data = points |> arrange(desc(status)), aes(colour = status, size = status)) +
  scale_colour_manual(values = c(
    "matched by st_join (14)"            = "#2ca02c",
    "orphan, but lying ON a line (26)"   = "#d62728",
    "orphan, genuinely off-network (146)" = "grey82"), name = NULL) +
  scale_size_manual(values = c(2.2, 2.4, 1.0), guide = "none") +
  theme_void() +
  theme(legend.position = "bottom")

[agila5]

Hi @gustavo-tm and sorry for my late reply. I checked your response and I have a few comments:

• If possible, could you please explain more precisely why you need / want to combine points and lines? I'm asking this question since I'm not really sure I understand the underlying objectives of these analyses.
• I would also like to point out that the points included into the points layer are only the OSM nodes with a significat tag attached (e.g., crossing = traffic_signal). If you really need to extract all OSM nodes you need to tune the OSM config file (check the vignette; I'm happy to provide an example if you want).
• I'm not able to run your code since, the following line errors:

id <- orphan_points$osm_id[d < 1e-6][1] #> "21069418"

It says that object d is not found.

• I don't understand what you mean by

However, for building incidence (which is my case), the precision is not extended to sf::st_coordinates sadly.

Can you please add more details or an example? Please notice that, as reported in the sf docs, the effects of st_precision are applied immediately before the GEOS routines are run inside the sf package and has no visible effect when printing the coordinates.

A few more comments based on your reply:

This probably needs further investigation, but it looks like adding sf::st_set_precision(1e7) to imported data should be a package default on osmextract::oe_get and others.

I'm not sure I like such ideas since the exact precision is really analysis-dependent and I'm not sure we should force any particular threshold. @Robinlovelace, what do you think?

On the other hand, this entire discussion about decimal precision could be somewhere to help people who want to recreate node relations, as exploding lines into nodes with sf::st_coordinates does cause this issue and manual snapping is required.

I like to idea and, as already mentioned, I'm happy to work jointly on updating the "Tips and tricks" vignette showcasing these issues.

[gustavo-tm]

Hi, @agila5,

You're right, I did not include the d variable in the code I sent, my mistake! I edited my previous answer to include the following

vertices <- lines |>
  st_coordinates() |>
  as.data.frame() |>
  distinct(X, Y) |>
  st_as_sf(coords = c("X", "Y"), crs = 4326)

d <- st_distance(orphan_points,
                 vertices[st_nearest_feature(orphan_points, vertices), ],
                 by_element = TRUE) |> as.numeric()

About the context

If possible, could you please explain more precisely why you need / want to combine points and lines? I'm asking this question since I'm not really sure I understand the underlying objectives of these analyses.

I am currently working on a research about traffic safety and OSM data is my main source of information about the road infrastructure. However, the data "as it comes" imposes many challenges for my study strategy.

• First, the length of the segments and they way they are split whenever there is a change in some tag is problematic for me. I need segments with relatively equal length and for them to break in relatively exogenous locations, not whenever there is a feature change. Thus, I need to manually reconstruct the lines geometries.
• Second, I need to calculate many variables that don't come by default. For example, I need to measure the amount and types of intersections: if they have crossings, signals, what road types participate in this intersection (primary x secondary or primary x residential). And for all those calculations it is imperative that I can perfectly join points and lines.

For this reason I am building a relational database from OSM source data that looks something like this

erDiagram
    ways {
        string id_way PK
        string name
        string name_alt
        string name_ref
        string type
        numeric lanes
        numeric speed_limit
        boolean oneway
        boolean roundabout
        LINESTRING geom
    }
 
 
    incidence {
        string id_way PK, FK
        integer position PK
        integer id_node FK
    }
 
    nodes {
        integer id_node PK
        string id_point FK "nullable"
        integer degree
        numeric X
        numeric Y
    }
 
    points {
        string id_point PK
        boolean traffic_signal
        boolean crossing
        boolean traffic_calming
        boolean speed_camera
        boolean roundabout
        POINT geom
    }
 
    junction_nodes {
        integer id_node PK, FK
        integer id_junction FK
        string role
    }
 
    junctions {
        integer id_junction PK
        integer n
        string types
        boolean roundabout
        integer traffic_signal
        integer crossing
        integer traffic_calming
        numeric degree
    }
 
    segments {
        integer id_chain PK "see segmentation diagram"
        integer id_segment PK
    }
 
    ways ||--o{ incidence : "contains"
    nodes ||--o{ incidence : "belongs to"
    points ||--o| nodes : "tags"
    nodes ||--o| junction_nodes : "member of"
    junctions ||--|{ junction_nodes : "groups"
    ways }o--o{ segments : "overlaps"
    nodes }o--o{ segments : "lies within"

Loading

And segments are built according to these parts

erDiagram
    ways {
        string id_way PK
    }
 
    nodes {
        integer id_node PK
    }
 
    roads {
        string id_way PK, FK
        integer id_road
    }
 
    chains {
        integer id_chain PK
        integer id_road
        numeric length
        LINESTRING geom
    }
 
    chains_nodes {
        integer id_chain PK, FK
        integer position_chain PK
        string id_way PK, FK
        integer id_road
        integer id_node FK
        numeric length_cumulative
    }
 
    chains_ways {
        integer id_chain PK, FK
        string id_way PK, FK
        integer id_road
        numeric way_start
        numeric way_end
    }
 
    segments {
        integer id_chain PK, FK
        integer id_segment PK
        integer id_road
        numeric length
        numeric speed_limit
        numeric lanes
        boolean lanes_change
        boolean oneway
        integer junctions
        integer traffic_signals
        integer crossings
        integer traffic_calming
        LINESTRING geom
    }
 
    segments_ways {
        integer id_chain PK, FK
        integer id_segment PK, FK
        string id_way PK, FK
        numeric overlap
    }
 
    segments_nodes {
        integer id_chain PK, FK
        integer id_segment PK, FK
        integer id_node PK, FK
    }
 
    ways ||--|| roads : "grouped into"
    roads ||--o{ chains : "split into"
    chains ||--|{ chains_nodes : "traversed by"
    nodes ||--o{ chains_nodes : "visited at"
    ways ||--o{ chains_nodes : "carries"
    chains ||--|{ chains_ways : "composed of"
    ways ||--o{ chains_ways : "spans"
    chains ||--|{ segments : "cut into"
    segments ||--|{ segments_ways : "composed of"
    ways ||--o{ segments_ways : "spans"
    segments ||--o{ segments_nodes : "contains"
    nodes ||--o{ segments_nodes : "assigned to"

Loading

And essentially, when I take the segments, I need to be able to merge it with information regarding points, junctions and ways.

About the following comment

I'm not sure I like such ideas since the exact precision is really analysis-dependent and I'm not sure we should force any particular threshold.

I don't think 1e7 is an analysis-dependent threshold, because it isn't a threshold at all. From what I understand, it's the resolution the source data already has, but R does not know about it and mistakenly assumes that is higher.

• a precision coarser than 1e7 (1e6, 1e5, ...) genuinely is an analysis choice: it discards real distinctions, and I agree the package shouldn't pick one for the user;
• a precision finer than 1e7, including the default of 0, doesn't buy any accuracy, because the accuracy ceiling is upstream. What it does buy is room for the sub-grid noise, which is exactly the failure in the reprex above, where many points silently fail to join.

[gustavo-tm]

There's also a practical reason this matters at import time rather than later: once you reproject, the grid is gone.

The 7-decimal grid only exists in the source CRS. If you set the precision first, sf snaps the coordinates onto that grid and then reprojects, so
everything stays consistent. If you reproject first, the coordinates become arbitrary metres and there is nothing left to round to, making the exact join impossible.

oe_get() |> st_transform(new_crs) is a completely natural first line to write, and it silently breaks exact join if precision was not set beforehand.

Besides, supporting the idea that this is not an arbitrary threshold, after st_transform the precision is set back to the default (0)

lines <- oe_get("ITS Leeds", layer = "lines", quiet = TRUE) |> st_geometry()

st_precision(lines) # 0

lines |> 
  st_set_precision(1e7) |> 
  st_precision() # 1e+07

lines |> 
  st_set_precision(1e7) |> 
  st_transform("EPSG:31983") |> 
  st_precision() # 0