CODEWITHSUNDEEP
rspatial
jesspi
jesspi

Reputation: 39

Calculating distance between spatial point and spatial line

I have a simple features object with spatial points (coordinates) in it, and I have a large SpatialLinesDataframe (representing roads). I have about 9000 points, and I am wanting to calculate the shortest distance between each point and the nearest road in the SpatialLinesDataframe. As rgdal and rgeos are now deprecated, I am unsure how I can best calculate these distances efficiently (ie gDistance won't work). Any help is appreciated!

Upvotes: 0

Views: 585

Answers (2)

margusl
margusl

Reputation: 17514

If it's OK for you to switch to sf package, you can combine or union your road network (presumably LINESTRING geometries) into a single MULTILINESTRING and use st_distance():

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
# as an example, Roxel road network from sfnetworks package
roxel_ <- sfnetworks::roxel
roxel_
#> Simple feature collection with 851 features and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 7.522594 ymin: 51.94151 xmax: 7.546705 ymax: 51.9612
#> Geodetic CRS:  WGS 84
#> # A tibble: 851 × 3
#>    name                  type                                           geometry
#>  * <chr>                 <fct>                                  <LINESTRING [°]>
#>  1 Havixbecker Strasse   residential      (7.533722 51.95556, 7.533461 51.95576)
#>  2 Pienersallee          secondary   (7.532442 51.95422, 7.53236 51.95377, 7.53…
#>  3 Schulte-Bernd-Strasse residential (7.532709 51.95209, 7.532823 51.95239, 7.5…
#>  4 <NA>                  path        (7.540063 51.94468, 7.539696 51.94479, 7.5…
#>  5 Welsingheide          residential       (7.537673 51.9475, 7.537614 51.94562)
#>  6 <NA>                  footway     (7.543791 51.94733, 7.54369 51.94686, 7.54…
#>  7 <NA>                  footway           (7.54012 51.94478, 7.539931 51.94514)
#>  8 <NA>                  path        (7.53822 51.94546, 7.538131 51.94549, 7.53…
#>  9 <NA>                  track       (7.540063 51.94468, 7.540338 51.94468, 7.5…
#> 10 <NA>                  track       (7.5424 51.94599, 7.54205 51.94629, 7.5419…
#> # ℹ 841 more rows

# sample points 
roxel_bbox <- 
  roxel_ |>
  st_bbox() |>
  st_as_sfc()
sample_points <- st_sf(id = 1:20, geometry = st_sample(roxel_bbox, 20))

plot(st_geometry(roxel_))
plot(st_geometry(sample_points), add = TRUE, pch = 19, col = "red")



# convert linestrings to a single multilinestring
roxel_u <- st_union(roxel_)
roxel_u
#> Geometry set for 1 feature 
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 7.522594 ymin: 51.94151 xmax: 7.546705 ymax: 51.9612
#> Geodetic CRS:  WGS 84
#> MULTILINESTRING ((7.533722 51.95556, 7.533461 5...

# add distance to a multilinestring for each point
sample_points$dist_to_road <- st_distance(sample_points, roxel_u)
sample_points
#> Simple feature collection with 20 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 7.522862 ymin: 51.94221 xmax: 7.545885 ymax: 51.96073
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>    id                  geometry   dist_to_road
#> 1   1   POINT (7.52293 51.9458) 191.023154 [m]
#> 2   2  POINT (7.537142 51.9531)  20.392243 [m]
#> 3   3 POINT (7.544911 51.94645)  76.427547 [m]
#> 4   4 POINT (7.523751 51.95759) 110.614399 [m]
#> 5   5 POINT (7.545885 51.96069) 173.888472 [m]
#> 6   6 POINT (7.523605 51.94263) 234.182639 [m]
#> 7   7 POINT (7.522862 51.95603) 101.840416 [m]
#> 8   8 POINT (7.525905 51.94221)  94.636355 [m]
#> 9   9  POINT (7.53068 51.95891)   3.990801 [m]
#> 10 10 POINT (7.525584 51.94721)   4.575924 [m]

Created on 2024-01-18 with reprex v2.0.2

Upvotes: 1

P. Luchs
P. Luchs

Reputation: 147

Given that the points are already an sf table, you can convert your SpatialLinesDataframe to an sf object and then use st_nearest_feature from the sf package to return the indexes of the nearest road for each point.

lines_sf <- sf::st_as_sf(lines_df)

# Creates a vector of indexes 
nearest_lines <- sf::st_nearest_feature(points, lines_sf)

If you wanted to add attributes (i.e., join the two) from the lines to the points df or vice versa you can use a spatial join like so:

# This will return the points with any attributes from the nearest lines attached 
lines_with_nearest_point <- st_join(points, lines_sf, join = st_nearest_feature)

Upvotes: 0

Related Questions