library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
我有一些数据,想要计算每个点(站点)沿着定义好的路径之间的距离。
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
)
我想计算34-35-...-39号站之间的距离,但需要沿着路径计算(站点编号决定顺序)。首先遇到的问题是路线(路径)不相连,站点也没有与路线相连。我首先尝试提取路径和站点:
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
)
我卡在了下一步。我首先尝试合并线条,然后使用st_snap()
将点捕捉到最近的线条,但没有成功。希望能得到任何帮助。
此内容由 reprex package (v2.0.1) 于2021-12-01创建