如何计算沿路径的点之间的距离
library(tidyverse)\nlibrary(sf)\n#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1\n我有一些数据，我想计算沿定义路径的每个点（站）之间的距离。
\ndat <-\n structure(\n list(\n name = c(\n "Untitled Path",\n "St34B",\n "St35N",\n "St36F",\n "St37N",\n "St38B",\n "Untitled Path",\n "St39N"\n ),\n description = c(\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_\n ),\n timestamp = structure(\n c(\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_\n ),\n class = c("POSIXct", "POSIXt"),\n tzone = ""\n ),\n begin = structure(\n c(\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_\n ),\n class = c("POSIXct", "POSIXt"),\n tzone = ""\n ),\n end = structure(\n c(\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_,\n NA_real_\n ),\n class = c("POSIXct", "POSIXt"),\n tzone = ""\n ),\n altitude_mode = c(\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_\n ),\n tessellate = c(\n 1L, -1L, -1L, -1L,\n -1L, -1L, 1L, -1L\n ),\n extrude = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),\n visibility = c(-1L, -1L, -1L, -1L, -1L, -1L, -1L, -1L),\n draw_order = c(\n NA_integer_,\n NA_integer_,\n NA_integer_,\n NA_integer_,\n NA_integer_,\n NA_integer_,\n NA_integer_,\n NA_integer_\n ),\n icon = c(\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_,\n NA_character_\n ),\n geometry = structure(\n list(\n structure(\n c(\n -213231.809501996,\n -205487.607705256,\n -784028.913066238,\n -708301.049327739\n ),\n .Dim = c(\n 2L,\n 2L\n ),\n class = c("XY", "LINESTRING", "sfg")\n ),\n structure(\n c(\n -213529.323058115,\n -785232.982945769\n ),\n class = c("XY", "POINT", "sfg")\n ),\n structure(\n c(\n -212176.423266777,\n -773238.391709674\n ),\n class = c("XY", "POINT", "sfg")\n ),\n structure(\n c(\n -210268.431741568,\n -756818.73172344\n ),\n class = c("XY", "POINT", "sfg")\n ),\n structure(\n c(\n -208050.517190725,\n -737973.862632309\n ),\n class = c("XY", "POINT", "sfg")\n ),\n structure(\n c(\n -206040.836893304,\n -709783.744787448\n ),\n class = c("XY", "POINT", "sfg")\n ),\n structure(\n c(\n -204426.676405507,\n -160265.400475699,\n -708310.127055397,\n -727750.877479657\n ),\n .Dim = c(\n 2L,\n 2L\n ),\n class = c("XY", "LINESTRING", "sfg")\n ),\n structure(\n c(\n -179260.597288432,\n -718361.477655825\n ),\n class = c("XY", "POINT", "sfg")\n )\n ),\n n_empty = 0L,\n crs = structure(\n 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\xc2\xb0N\\"],\\n BBOX[60,-180,90,180]],\\n ID[\\"EPSG\\",3411]]"),\n class = "crs"\n ),\n class = c(\n "sfc_GEOMETRY",\n "sfc"\n ),\n precision = 0,\n bbox = structure(\n c(\n xmin = -213529.323058115,\n ymin = -785232.982945769,\n xmax = -160265.400475699,\n ymax = -708301.049327739\n ),\n class = "bbox"\n ),\n classes = c(\n "LINESTRING",\n "POINT",\n "POINT",\n "POINT",\n "POINT",\n "POINT",\n "LINESTRING",\n "POINT"\n )\n )\n ),\n row.names = c(\n NA,\n 8L\n ),\n sf_column = "geometry",\n agr = structure(\n c(\n name = NA_integer_,\n description = NA_integer_,\n timestamp = NA_integer_,\n begin = NA_integer_,\n end = NA_integer_,\n altitude_mode = NA_integer_,\n tessellate = NA_integer_,\n extrude = NA_integer_,\n visibility = NA_integer_,\n draw_order = NA_integer_,\n icon = NA_integer_\n ),\n class = "factor",\n .Label = c(\n "constant",\n "aggregate", "identity"\n )\n ),\n class = c("sf", "data.frame")\n )\n\ndat\n#> Simple feature collection with 8 features and 11 fields\n#> Geometry type: GEOMETRY\n#> Dimension: XY\n#> Bounding box: xmin: -213529.3 ymin: -785233 xmax: -160265.4 ymax: -708301\n#> Projected CRS: NSIDC Sea Ice Polar Stereographic North\n#> name description timestamp begin end altitude_mode tessellate\n#> 1 Untitled Path <NA> <NA> <NA> <NA> <NA> 1\n#> 2 St34B <NA> <NA> <NA> <NA> <NA> -1\n#> 3 St35N <NA> <NA> <NA> <NA> <NA> -1\n#> 4 St36F <NA> <NA> <NA> <NA> <NA> -1\n#> 5 St37N <NA> <NA> <NA> <NA> <NA> -1\n#> 6 St38B <NA> <NA> <NA> <NA> <NA> -1\n#> 7 Untitled Path <NA> <NA> <NA> <NA> <NA> 1\n#> 8 St39N <NA> <NA> <NA> <NA> <NA> -1\n#> extrude visibility draw_order icon geometry\n#> 1 0 -1 NA <NA> LINESTRING (-213231.8 -7840...\n#> 2 0 -1 NA <NA> POINT (-213529.3 -785233)\n#> 3 0 -1 NA <NA> POINT (-212176.4 -773238.4)\n#> 4 0 -1 NA <NA> POINT (-210268.4 -756818.7)\n#> 5 0 -1 NA <NA> POINT (-208050.5 -737973.9)\n#> 6 0 -1 NA <NA> POINT (-206040.8 -709783.7)\n#> 7 0 -1 NA <NA> LINESTRING (-204426.7 -7083...\n#> 8 0 -1 NA <NA> POINT (-179260.6 -718361.5)\n\nggplot() +\n geom_sf(data = dat) +\n geom_sf_text(\n data = dat,\n aes(label = name),\n size = 3,\n hjust = 0\n )\n
我想计算车站 34 - 35 - \xe2\x80\xa6 - 39\n之间的距离，但沿着路径（车站编号决定顺序）。第一个问题\n我看到的是线路（路径）未连接并且车站未\n连接到线路。\n我首先尝试提取路径和车站：
\nstations <- dat %>%\n filter(str_starts(name, "St"))\n\npaths <- dat %>%\n filter(str_starts(name, "Untitled"))\n\nggplot() +\n geom_sf(data = paths, color = "red") +\n geom_sf(data = stations, color = "blue") +\n geom_sf_text(\n data = stations,\n aes(label = name),\n color = "blue",\n size = 3,\n hjust = 0\n )\n
我被困在接下来的步骤中。我首先尝试合并这些线，然后将点捕捉到最近的线，但st_snap()没有成功。任何帮助表示赞赏。
由reprex 包于 2021 年 12 月 1 日创建(v2.0.1)
\n请查找详细的表示，它使用sf、sfnetworks、units和库为您的请求提供解决方案。dplyrggplot2
雷普莱克斯
- 步骤 1：仅基于“连接线（ieedges）”创建“sfnetworks”对象
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)
- 步骤 2：创建仅包含点（即节点）的“sf”对象
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)
步骤 3：将“sf”对象的节点添加到“network”中
1. 代码
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. 网络规格
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. 网络可视化
# 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
)
- 步骤 4：计算网络中每个节点之间的边的长度并创建数据帧
distances（即 tibble 类）
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
distances步骤 5：用节点名称替换数据帧的“from”和“to”列的 ID
1. 提取节点名称并将它们映射到距离数据帧的 id
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. 修改数据框distances以使用两个函数获取“from”和“to”列中的节点名称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.最终输出
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
由reprex 包于 2021 年 12 月 6 日创建(v2.0.1)