Philippe Massicotte
Reputation: 1529
How to calculate distances between points along a path
library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
I have some data where I would like to calculate the distances between each point (station) along defined paths.
dat <-
structure(
list(
name = c(
"Untitled Path",
"St34B",
"St35N",
"St36F",
"St37N",
"St38B",
"Untitled Path",
"St39N"
),
description = c(
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_
),
timestamp = structure(
c(
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_
),
class = c("POSIXct", "POSIXt"),
tzone = ""
),
begin = structure(
c(
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_
),
class = c("POSIXct", "POSIXt"),
tzone = ""
),
end = structure(
c(
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_,
NA_real_
),
class = c("POSIXct", "POSIXt"),
tzone = ""
),
altitude_mode = c(
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_
),
tessellate = c(
1L, -1L, -1L, -1L,
-1L, -1L, 1L, -1L
),
extrude = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),
visibility = c(-1L, -1L, -1L, -1L, -1L, -1L, -1L, -1L),
draw_order = c(
NA_integer_,
NA_integer_,
NA_integer_,
NA_integer_,
NA_integer_,
NA_integer_,
NA_integer_,
NA_integer_
),
icon = c(
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_,
NA_character_
),
geometry = structure(
list(
structure(
c(
-213231.809501996,
-205487.607705256,
-784028.913066238,
-708301.049327739
),
.Dim = c(
2L,
2L
),
class = c("XY", "LINESTRING", "sfg")
),
structure(
c(
-213529.323058115,
-785232.982945769
),
class = c("XY", "POINT", "sfg")
),
structure(
c(
-212176.423266777,
-773238.391709674
),
class = c("XY", "POINT", "sfg")
),
structure(
c(
-210268.431741568,
-756818.73172344
),
class = c("XY", "POINT", "sfg")
),
structure(
c(
-208050.517190725,
-737973.862632309
),
class = c("XY", "POINT", "sfg")
),
structure(
c(
-206040.836893304,
-709783.744787448
),
class = c("XY", "POINT", "sfg")
),
structure(
c(
-204426.676405507,
-160265.400475699,
-708310.127055397,
-727750.877479657
),
.Dim = c(
2L,
2L
),
class = c("XY", "LINESTRING", "sfg")
),
structure(
c(
-179260.597288432,
-718361.477655825
),
class = c("XY", "POINT", "sfg")
)
),
n_empty = 0L,
crs = structure(
list(input = "EPSG:3411", wkt = "PROJCRS[\"NSIDC Sea Ice Polar Stereographic North\",\n BASEGEOGCRS[\"Unspecified datum based upon the Hughes 1980 ellipsoid\",\n DATUM[\"Not specified (based on Hughes 1980 ellipsoid)\",\n ELLIPSOID[\"Hughes 1980\",6378273,298.279411123064,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4054]],\n CONVERSION[\"US NSIDC Sea Ice polar stereographic north\",\n METHOD[\"Polar Stereographic (variant B)\",\n ID[\"EPSG\",9829]],\n PARAMETER[\"Latitude of standard parallel\",70,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8832]],\n PARAMETER[\"Longitude of origin\",-45,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8833]],\n PARAMETER[\"False easting\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting (X)\",south,\n MERIDIAN[45,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",south,\n MERIDIAN[135,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"World - N hemisphere - north of 60°N\"],\n BBOX[60,-180,90,180]],\n ID[\"EPSG\",3411]]"),
class = "crs"
),
class = c(
"sfc_GEOMETRY",
"sfc"
),
precision = 0,
bbox = structure(
c(
xmin = -213529.323058115,
ymin = -785232.982945769,
xmax = -160265.400475699,
ymax = -708301.049327739
),
class = "bbox"
),
classes = c(
"LINESTRING",
"POINT",
"POINT",
"POINT",
"POINT",
"POINT",
"LINESTRING",
"POINT"
)
)
),
row.names = c(
NA,
8L
),
sf_column = "geometry",
agr = structure(
c(
name = NA_integer_,
description = NA_integer_,
timestamp = NA_integer_,
begin = NA_integer_,
end = NA_integer_,
altitude_mode = NA_integer_,
tessellate = NA_integer_,
extrude = NA_integer_,
visibility = NA_integer_,
draw_order = NA_integer_,
icon = NA_integer_
),
class = "factor",
.Label = c(
"constant",
"aggregate", "identity"
)
),
class = c("sf", "data.frame")
)
dat
#> Simple feature collection with 8 features and 11 fields
#> Geometry type: GEOMETRY
#> Dimension: XY
#> Bounding box: xmin: -213529.3 ymin: -785233 xmax: -160265.4 ymax: -708301
#> Projected CRS: NSIDC Sea Ice Polar Stereographic North
#> name description timestamp begin end altitude_mode tessellate
#> 1 Untitled Path <NA> <NA> <NA> <NA> <NA> 1
#> 2 St34B <NA> <NA> <NA> <NA> <NA> -1
#> 3 St35N <NA> <NA> <NA> <NA> <NA> -1
#> 4 St36F <NA> <NA> <NA> <NA> <NA> -1
#> 5 St37N <NA> <NA> <NA> <NA> <NA> -1
#> 6 St38B <NA> <NA> <NA> <NA> <NA> -1
#> 7 Untitled Path <NA> <NA> <NA> <NA> <NA> 1
#> 8 St39N <NA> <NA> <NA> <NA> <NA> -1
#> extrude visibility draw_order icon geometry
#> 1 0 -1 NA <NA> LINESTRING (-213231.8 -7840...
#> 2 0 -1 NA <NA> POINT (-213529.3 -785233)
#> 3 0 -1 NA <NA> POINT (-212176.4 -773238.4)
#> 4 0 -1 NA <NA> POINT (-210268.4 -756818.7)
#> 5 0 -1 NA <NA> POINT (-208050.5 -737973.9)
#> 6 0 -1 NA <NA> POINT (-206040.8 -709783.7)
#> 7 0 -1 NA <NA> LINESTRING (-204426.7 -7083...
#> 8 0 -1 NA <NA> POINT (-179260.6 -718361.5)
ggplot() +
geom_sf(data = dat) +
geom_sf_text(
data = dat,
aes(label = name),
size = 3,
hjust = 0
)
I would like to calculate the distance between stations 34 - 35 - … - 39 but along the path (station numbers determine the order).The first problems I see is that the lines (paths) are not connected and the stations are not connected to the lines. I first tried to extract the paths and the stations:
stations <- dat %>%
filter(str_starts(name, "St"))
paths <- dat %>%
filter(str_starts(name, "Untitled"))
ggplot() +
geom_sf(data = paths, color = "red") +
geom_sf(data = stations, color = "blue") +
geom_sf_text(
data = stations,
aes(label = name),
color = "blue",
size = 3,
hjust = 0
)
I am stuck on the next steps. I first tried to merge the lines and then
snap the points to the closest line using st_snap() without success. Any
help is appreciated.
Created on 2021-12-01 by the reprex package (v2.0.1)
Upvotes: 3
Views: 1245
Answers (1)
lovalery
Reputation: 4652
Please find a detailed reprex that provides a solution to your request using the sf, sfnetworks, units, dplyr and ggplot2 libraries.
Reprex
- STEP 1: Create a 'sfnetworks' object only based 'on connected lines(i.e.edges)
library(sf)
library(units)
library(sfnetworks)
options(sfn_max_print_active = 15, sfn_max_print_inactive = 15)
library(dplyr)
library(ggplot2)
network <- dat %>%
filter(st_geometry_type(.) == "LINESTRING") %>% # selects only the lines from 'sf' object 'dat'
st_snap(.,., tolerance = 10000) %>% # coerces the snapping using a big tolerance value!
as_sfnetwork() # creates the network
autoplot(network)
- STEP 2: Create a 'sf' object with only points (i.e. nodes)
nodes <- dat %>%
filter(st_geometry_type(.) == "POINT")
nodes
#> Simple feature collection with 6 features and 11 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -213529.3 ymin: -785233 xmax: -179260.6 ymax: -709783.7
#> Projected CRS: NSIDC Sea Ice Polar Stereographic North
#> name description timestamp begin end altitude_mode tessellate extrude
#> 1 St34B <NA> <NA> <NA> <NA> <NA> -1 0
#> 2 St35N <NA> <NA> <NA> <NA> <NA> -1 0
#> 3 St36F <NA> <NA> <NA> <NA> <NA> -1 0
#> 4 St37N <NA> <NA> <NA> <NA> <NA> -1 0
#> 5 St38B <NA> <NA> <NA> <NA> <NA> -1 0
#> 6 St39N <NA> <NA> <NA> <NA> <NA> -1 0
#> visibility draw_order icon geometry
#> 1 -1 NA <NA> POINT (-213529.3 -785233)
#> 2 -1 NA <NA> POINT (-212176.4 -773238.4)
#> 3 -1 NA <NA> POINT (-210268.4 -756818.7)
#> 4 -1 NA <NA> POINT (-208050.5 -737973.9)
#> 5 -1 NA <NA> POINT (-206040.8 -709783.7)
#> 6 -1 NA <NA> POINT (-179260.6 -718361.5)
STEP 3: Add the nodes of the 'sf' object into the 'network'
1. Code
new_network <- network %>%
st_network_blend(., nodes, tolerance = 10000) %>% # snap the nodes on the network based on the given tolerance
filter(.,!is.na(name)) %>% # keeps only the nodes from the 'sf' object 'nodes'
st_as_sf %>% # convert into sf object (mandatory step for the next one to work properly)
as_sfnetwork(., edges_as_lines = TRUE) # reconstructs the network only with the nodes from the 'sf' object 'nodes'
#> Warning: st_network_blend assumes attributes are constant over geometries
2. Specifications of the network
new_network
#> # A sfnetwork with 6 nodes and 5 edges
#> #
#> # CRS: EPSG:3411
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Node Data: 6 x 12 (active)
#> # Geometry type: POINT
#> # Dimension: XY
#> # Bounding box: xmin: -213231.8 ymin: -784028.9 xmax: -179639.4 ymax:
#> # -709824.4
#> name description timestamp begin end
#> <chr> <chr> <dttm> <dttm> <dttm>
#> 1 St34B <NA> NA NA NA
#> 2 St35N <NA> NA NA NA
#> 3 St36F <NA> NA NA NA
#> 4 St37N <NA> NA NA NA
#> 5 St38B <NA> NA NA NA
#> 6 St39N <NA> NA NA NA
#> # ... with 7 more variables: altitude_mode <chr>, tessellate <int>,
#> # extrude <int>, visibility <int>, draw_order <int>, icon <chr>,
#> # geometry <POINT [m]>
#> #
#> # Edge Data: 5 x 3
#> # Geometry type: LINESTRING
#> # Dimension: XY
#> # Bounding box: xmin: -213231.8 ymin: -784028.9 xmax: -179639.4 ymax:
#> # -709824.4
#> from to geometry
#> <int> <int> <LINESTRING [m]>
#> 1 1 2 (-213231.8 -784028.9, -212128.8 -773243.3)
#> 2 2 3 (-212128.8 -773243.3, -210447.3 -756800.4)
#> 3 3 4 (-210447.3 -756800.4, -208517.2 -737926.1)
#> 4 4 5 (-208517.2 -737926.1, -205643.4 -709824.4)
#> 5 5 6 (-205643.4 -709824.4, -179639.4 -719222)
3. Visualization of the network
# option 1 with autoplot:
autoplot(new_network) +
geom_sf_text(
data = st_as_sf(new_network),
aes(label = name),
size = 3,
hjust = 0
)
# if you prefer, option 2 with only ggplot:
ggplot() +
geom_sf(data = st_as_sf(new_network, "edges"), col = "grey50") +
geom_sf(data = st_as_sf(new_network, "nodes")) +
geom_sf_text(
data = st_as_sf(new_network),
aes(label = name),
size = 3,
hjust = 0
)
- STEP 4: Computes the length of edges between each node along the network and creates the dataframe
distances(i.e. tibble class)
distances <- new_network %>%
activate("edges") %>%
mutate(length = set_units(edge_length(),km)) %>%
st_as_sf() %>%
st_drop_geometry
distances
#> # A tibble: 5 x 3
#> from to length
#> * <int> <int> [km]
#> 1 1 2 10.8
#> 2 2 3 16.5
#> 3 3 4 19.0
#> 4 4 5 28.2
#> 5 5 6 27.6
- STEP 5: Replace ids of columns "from" and "to" of the
distancesdataframe by the names of nodes
1. Extract names of nodes and map them to the id's of distances dataframe
names_id <- new_network %>%
activate("nodes") %>%
st_as_sf() %>%
mutate(ID = seq(name)) %>%
select(., c("ID", "name")) %>%
st_drop_geometry
names_id
#> # A tibble: 6 x 2
#> ID name
#> * <int> <chr>
#> 1 1 St34B
#> 2 2 St35N
#> 3 3 St36F
#> 4 4 St37N
#> 5 5 St38B
#> 6 6 St39N
2. Modify the dataframe distances to get the names of nodes in 'from' and 'to' columns using two left_join()
distances <- left_join(distances, names_id, by = c("from" = "ID")) %>%
mutate(from = name) %>%
select(-name) %>%
left_join(., names_id, by = c("to" = "ID")) %>%
mutate(to = name) %>%
select(-name)
3. Final output
distances
#> # A tibble: 5 x 3
#> from to length
#> <chr> <chr> [km]
#> 1 St34B St35N 10.8
#> 2 St35N St36F 16.5
#> 3 St36F St37N 19.0
#> 4 St37N St38B 28.2
#> 5 St38B St39N 27.6
Created on 2021-12-06 by the reprex package (v2.0.1)
Upvotes: 5
Related Questions
- How to measure distance between points in separate data frames?
- Using R, how to calculate the distance from one point to a line?
- How to Calculate Distance between points from a reference point using R?
- R - How to calculate distance between two sets of coordinate points?
- Measuring the distance of all points along a line in R (linestring in sf)
- How to calculate geographic distance between two points along a line in R?
- R - Calculate distance between two points along a polyline
- Calculating the distance between points in R
- Find Distance Between Locations in R Data Frame
- How to find whole distance between two points in a curved line in R?