使用SF将点捕捉到线段上的最近点
我想使用将点捕捉到道路段上最接近的点sf::st_snap。但是，该函数似乎返回了错误的结果，它将我的观点捕捉到了路段起点处的一点。有想法该怎么解决这个吗？
下面提供了可重现的示例，包括对使用sf::st_snapvs 时得到的结果的比较maptools::snapPointsToLines
使用 sf::st_snap
# Max distance
cut_dist = 200 # meters
# snap points to closest road
new_point <- sf::st_snap(point1, roads, tolerance = cut_dist)
# view points on the map
mapView(point1, color="red") + mapView( st_buffer(point1, dist = cut_dist)) + mapView(new_point) + mapView(roads)
# Distance between pont1 and new_point
st_distance( point1, new_point)
> 1591 meters # note this is above the set maximun distance
使用maptools::snapPointsToLines（我期望的结果）
# convert sf to sp
point_sp <- as_Spatial(point1)
roads_sp <- as_Spatial(roads)
# snap points
new_point_sp <- snapPointsToLines(point_sp, roads_sp, maxDist = cut_dist)
# view points on the map
mapView(point1, color="red") + mapView( st_buffer(point1, dist = cut_dist)) + mapView(new_point_sp) + mapView(roads)
# Distance between pont1 and new_point
spDistsN1( point_sp, new_point_sp)
> 116 meters
资料库
library(sf)
library(mapview)
library(maptools)
library(sp)
point1 <- structure(list(idhex = 9L, geometry = structure(list(structure(c(665606.970079183,
6525003.41418009), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(xmin = 665606.970079183,
ymin = 6525003.41418009, xmax = 665606.970079183, ymax = 6525003.41418009
), class = "bbox"), crs = structure(list(epsg = 32633L, proj4string = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(idhex = NA_integer_), .Label = c("constant",
"aggregate", "identity"), class = "factor"), row.names = 2L, class = c("sf",
"data.table", "data.frame"))
roads <- structure(list(id = 139885, osm_id = 250886593, geometry = structure(list(
structure(c(665387.589147313, 665367.867159783, 665363.008143169,
665363.051860059, 665363.308104069, 665366.519781353, 665368.635421323,
665370.846894641, 665370.829724196, 665367.910645335, 665361.777524054,
665355.967776345, 665351.649946698, 665343.44353825, 665334.917779131,
665313.306069501, 665309.001351385, 665310.66019677, 665313.528620709,
665341.742306731, 665351.854389331, 665354.981775569, 665360.254611229,
665365.006104512, 665379.034588542, 665394.435039616, 665409.282519288,
665410.676785182, 665448.890797438, 665458.917562631, 665471.042094353,
665485.485001236, 665495.899212422, 665504.535684257, 665509.674854913,
665506.145837246, 665483.727146874, 665481.426949686, 665462.311063365,
665445.215460573, 665450.424049663, 665450.837897892, 665491.036360788,
665491.419140717, 665469.507518623, 665458.677850808, 665455.926197775,
665462.873809047, 665460.283684337, 665426.046702616, 665396.279686035,
665368.373253059, 665357.878521323, 665304.347529357, 665221.04051616,
665170.777462125, 665144.670345016, 665106.030568334, 665073.2789218,
665018.208956171, 664947.693178271, 664921.708297412, 664861.659061389,
664797.900403384, 664745.001666066, 664730.200174759, 664717.892651619,
664706.473711845, 664697.750102392, 664688.215719591, 664681.544531593,
664672.960647368, 664665.064067202, 664636.446517023, 664622.930521655,
664518.065243846, 664442.725560545, 664423.048166559, 664411.132259582,
664407.05972929, 664398.364646172, 664391.348502443, 664382.558239303,
664372.012526058, 664354.354954718, 664332.995014599, 664311.609706282,
664271.102641808, 664228.816287751, 664150.088321471, 664069.895400484,
6526138.02793883, 6526135.40749336, 6526130.11578605, 6526111.34403368,
6526087.4978365, 6526054.13445288, 6526022.49962268, 6525982.74830288,
6525959.40435839, 6525944.55197219, 6525918.33886077, 6525894.18611795,
6525874.55473851, 6525840.53410542, 6525813.96628006, 6525767.42907088,
6525745.21917638, 6525733.51582599, 6525713.24841331, 6525627.57847652,
6525608.06984863, 6525568.30170735, 6525550.71644271, 6525539.76231607,
6525491.25651378, 6525446.12690364, 6525433.36256694, 6525431.23562504,
6525372.98235432, 6525354.13376808, 6525331.3288195, 6525309.59511696,
6525293.92174422, 6525270.21980161, 6525256.11455612, 6525228.35885783,
6525217.10943051, 6525215.95489587, 6525195.91355696, 6525158.79257025,
6525134.01851773, 6525131.70940566, 6525050.96446632, 6524950.68358502,
6524851.23226232, 6524806.24052727, 6524749.34394609, 6524714.63617193,
6524660.07336072, 6524612.21010524, 6524583.84484865, 6524562.03540982,
6524557.38094998, 6524533.67136837, 6524510.74454804, 6524495.56823805,
6524486.9387399, 6524475.63373441, 6524465.4404841, 6524468.04929815,
6524475.95178632, 6524478.86036788, 6524470.76472937, 6524447.96214429,
6524448.06967557, 6524443.4855897, 6524435.86812114, 6524425.93373791,
6524417.67487537, 6524409.79262886, 6524399.64960133, 6524378.79085156,
6524360.33496349, 6524303.24355601, 6524302.70486651, 6524293.01335665,
6524290.81442892, 6524298.30279414, 6524309.46697681, 6524313.27442914,
6524337.22831533, 6524364.43083297, 6524376.27944935, 6524382.92319852,
6524389.6474774, 6524406.74565716, 6524430.82326744, 6524462.46041311,
6524492.20009833, 6524544.74318075, 6524591.10483188), .Dim = c(91L,
2L), class = c("XY", "LINESTRING", "sfg"))), class = c("sfc_LINESTRING",
"sfc"), precision = 0, bbox = structure(c(xmin = 664069.895400484,
ymin = 6524290.81442892, xmax = 665509.674854913, ymax = 6526138.02793883
), class = "bbox"), crs = structure(list(epsg = 32633L, proj4string = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"), class = "crs"), n_empty = 0L)), row.names = 139885L, class = c("sf",
"data.frame"), sf_column = "geometry", agr = structure(c(id = NA_integer_,
osm_id = NA_integer_), .Label = c("constant", "aggregate", "identity"
), class = "factor"))
鉴于这个问题几乎在四年前就得到了回答，我想我应该添加一种替代的、更快的方法，这在当时可能是不可能的。@TimSalabim 的st_snap_points()功能很棒，但如果您有sf collection一些要点，直接st_nearest_points()在sf collection.
正如 Tim 和其他人在讨论中提到的，st_nearest_points()返回一个LINESTRING对象 - 一开始有点令人困惑，但可以很容易地变成POINT. 听起来 中点的顺序LINESTRING总是从第一个几何图形到第二个几何图形，因此您可以选择 的起点或终点LINESTRING，具体取决于您将几何图形放入的顺序st_nearest_points()。
看?st_nearest_points：
价值
an sfc object with all two-point LINESTRING geometries of point pairs from the first to the second geometry, of length x * y, with y cycling fastest. See examples for ideas how to convert these to POINT geometries.
Using your example, we can find the closest point on the road by making the road the second object, and retrieving the endpoint (second point) of the LINESTRING as follows:
# find the LINESTRING object which is the shortest line between the two geometries
st_nearest_points(point1, roads) %>%
{. ->> my_linestring}
my_linestring
# Geometry set for 1 feature
# Geometry type: LINESTRING
# Dimension: XY
# Bounding box: xmin: 665491.2 ymin: 6525003 xmax: 665607 ymax: 6525003
# Projected CRS: WGS 84 / UTM zone 33N
# LINESTRING (665607 6525003, 665491.2 6525003)
# then, convert the LINESTRING to POINTs
# and, pull out the second point, because we want the point on the 'roads' object,
# which was supplied second to st_nearest_points()
my_linestring %>%
st_cast('POINT') %>%
.[2] %>%
{. ->> closest_point}
closest_point
# Geometry set for 1 feature
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 665491.2 ymin: 6525003 xmax: 665491.2 ymax: 6525003
# Projected CRS: WGS 84 / UTM zone 33N
# POINT (665491.2 6525003)
And we can look at it on a map to confirm:
library(tidyverse)
ggplot()+
geom_sf(data = roads, col = 'red')+
geom_sf(data = point1, col = 'black')+
geom_sf(data = closest_point, col = 'blue')
To find the distance between the points we use st_distance():
# work out the distance
st_distance(point1, closest_point)[1]
# 115.7514 [m]
To incorporate a distance threshold into the "snapping", we can make a function which includes an if statement:
# we can combine the code and snap the point conditionally
f1 <- function(x, y, max_dist){
st_nearest_points(x, y) %>%
{. ->> my_linestring} %>%
st_cast('POINT') %>%
.[2] %>%
{. ->> closest_point} %>%
st_distance(., x) %>%
as.numeric %>%
{. ->> my_distance}
if(my_distance <= max_dist){
return(closest_point)
} else {
return(st_geometry(point1))
}
}
# run the function with threshold of 400 m (it should snap)
f1(point1, roads, max_dist = 400) %>%
{. ->> p1_400}
# run the function with a threshold of 50 m (it should NOT snap)
f1(point1, roads, max_dist = 50) %>%
{. ->> p1_50}
# plot it to see the results
ggplot()+
geom_sf(data = roads, col = 'red')+
geom_sf(data = p1_400, col = 'blue')+
geom_sf(data = p1_50)
# the blue point is on the road because it was within the threshold
# the black point is in the original spot because it was outside the threshold
Working on an sf collection rather than a single point
Now, I'll assume that many instances will require this to be done on multiple points rather than a single point, and we can modify it slightly to use this approach on an sf collection with multiple points. To demonstrate this, first make up 100 random points within the bbox of roads:
set.seed(69)
# let's make some random points around the road
tibble(
x = runif(n = 100, min = st_bbox(roads)$xmin, max = st_bbox(roads)$xmax),
y = runif(n = 100, min = st_bbox(roads)$ymin, max = st_bbox(roads)$ymax),
id = 1:100
) %>%
st_as_sf(
x = .,
coords = c('x', 'y'),
crs = st_crs(roads)
) %>%
{. ->> my_points}
my_points
# Simple feature collection with 100 features and 1 field
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 664090.7 ymin: 6524310 xmax: 665500.6 ymax: 6526107
# Projected CRS: WGS 84 / UTM zone 33N
# # A tibble: 100 x 2
# id geometry
# * <int> <POINT [m]>
# 1 1 (664834.1 6525742)
# 2 2 (665176.8 6524370)
# 3 3 (664999.9 6525528)
# 4 4 (665315.7 6525242)
# 5 5 (664601 6525464)
# 6 6 (665320.7 6525600)
# 7 7 (664316.2 6525065)
# 8 8 (665204 6526025)
# 9 9 (664319.8 6524335)
# 10 10 (664101.7 6525830)
# # ... with 90 more rows
# now show the points on a figure
ggplot()+
geom_sf(data = roads, col = 'red')+
geom_sf(data = my_points, shape = 1)
To find the closest point on roads to each of the random points, we use the same functions as the single point example above, but we save the LINESTRING as a column for each point. Then, instead of retrieving just the second point of the LINESTRING, we pull out every second point of the collection (column) of LINETRINGs. You have the option of pulling out the second point of every LINESTRING on a row-by-row basis, by including rowwise() before mutate(), but this makes for much slower computation. Instead, by working on the sf collection as a whole and pulling out every second point (the second point from each row), you get the result much much faster. This is how they demonstrate it in the example of ?st_nearest_points.
So basically, we make new columns in the sf collection using mutate() which includes the LINESTRING and the closest POINT.
# now, let's find the closest points
my_points %>%
mutate(
my_linestring = st_nearest_points(geometry, roads),
closest_point = st_cast(my_linestring, 'POINT')[seq(2, nrow(.)*2, 2)],
) %>%
{. ->> closest_points}
closest_points
# Simple feature collection with 100 features and 1 field
# Active geometry column: geometry
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 664090.7 ymin: 6524310 xmax: 665500.6 ymax: 6526107
# Projected CRS: WGS 84 / UTM zone 33N
# # A tibble: 100 x 4
# id geometry my_linestring closest_point
# * <int> <POINT [m]> <LINESTRING [m]> <POINT [m]>
# 1 1 (664834.1 6525742) (664834.1 6525742, 665309 6525745) (665309 6525745)
# 2 2 (665176.8 6524370) (665176.8 6524370, 665142.8 6524486) (665142.8 6524486)
# 3 3 (664999.9 6525528) (664999.9 6525528, 665337.9 6525639) (665337.9 6525639)
# 4 4 (665315.7 6525242) (665315.7 6525242, 665454.2 6525178) (665454.2 6525178)
# 5 5 (664601 6525464) (664601 6525464, 665317.8 6525700) (665317.8 6525700)
# 6 6 (665320.7 6525600) (665320.7 6525600, 665348.4 6525615) (665348.4 6525615)
# 7 7 (664316.2 6525065) (664316.2 6525065, 664069.9 6524591) (664069.9 6524591)
# 8 8 (665204 6526025) (665204 6526025, 665367.7 6526036) (665367.7 6526036)
# 9 9 (664319.8 6524335) (664319.8 6524335, 664354.4 6524390) (664354.4 6524390)
# 10 10 (664101.7 6525830) (664101.7 6525830, 665309 6525745) (665309 6525745)
# # ... with 90 more rows
# and have a look at a figure
ggplot()+
geom_sf(data = roads, col = 'red')+
geom_sf(data = my_points, shape = 1)+
geom_sf(data = closest_points$my_linestring, linetype = 'dotted')+
geom_sf(data = closest_points$closest_point, shape = 1, col = 'blue')
如果您想在分析中包含距离阈值，我们会计算点之间的距离（或 的长度my_linestring），并创建一个新的几何列，有条件地包含原始几何或新的捕捉点，具体取决于distance（请参阅此处） 。我们使用ifelse()400 m 的阈值。
# now, use the closest_point if the distance is within a threshold, otherwise keep the original point
closest_points %>%
mutate(
distance = st_distance(geometry, roads)[,1],
distance2 = st_length(my_linestring), # not used, demo purposes only
snapped_point_cond = st_sfc(ifelse(as.numeric(distance) <= 400, st_geometry(closest_point), geometry), crs = st_crs(roads))
) %>%
{. ->> closest_points_threshold}
closest_points_threshold
# Simple feature collection with 100 features and 3 fields
# Active geometry column: geometry
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 664090.7 ymin: 6524310 xmax: 665500.6 ymax: 6526107
# Projected CRS: WGS 84 / UTM zone 33N
# # A tibble: 100 x 7
# id geometry my_linestring closest_point distance distance2 snapped_point_cond
# * <int> <POINT [m]> <LINESTRING [m]> <POINT [m]> [m] [m] <POINT [m]>
# 1 1 (664834.1 6525742) (664834.1 6525742, 665309 6525745) (665309 6525745) 475. 475. (664834.1 6525742)
# 2 2 (665176.8 6524370) (665176.8 6524370, 665142.8 6524486) (665142.8 6524486) 121. 121. (665142.8 6524486)
# 3 3 (664999.9 6525528) (664999.9 6525528, 665337.9 6525639) (665337.9 6525639) 356. 356. (665337.9 6525639)
# 4 4 (665315.7 6525242) (665315.7 6525242, 665454.2 6525178) (665454.2 6525178) 153. 153. (665454.2 6525178)
# 5 5 (664601 6525464) (664601 6525464, 665317.8 6525700) (665317.8 6525700) 755. 755. (664601 6525464)
# 6 6 (665320.7 6525600) (665320.7 6525600, 665348.4 6525615) (665348.4 6525615) 31.2 31.2 (665348.4 6525615)
# 7 7 (664316.2 6525065) (664316.2 6525065, 664069.9 6524591) (664069.9 6524591) 534. 534. (664316.2 6525065)
# 8 8 (665204 6526025) (665204 6526025, 665367.7 6526036) (665367.7 6526036) 164. 164. (665367.7 6526036)
# 9 9 (664319.8 6524335) (664319.8 6524335, 664354.4 6524390) (664354.4 6524390) 64.5 64.5 (664354.4 6524390)
# 10 10 (664101.7 6525830) (664101.7 6525830, 665309 6525745) (665309 6525745) 1210. 1210. (664101.7 6525830)
# # ... with 90 more rows
# and plot it
ggplot()+
geom_sf(data = roads %>% st_buffer(400), col = 'red', fill = NA)+
geom_sf(data = roads, col = 'red')+
geom_sf(data = closest_points_threshold$snapped_point_cond, shape = 1, col = 'blue')+
geom_sf(data = closest_points_threshold %>% filter(as.numeric(distance) <= 400) %>% .$my_linestring, linetype = 'dotted')+
geom_sf(data = closest_points_threshold %>% filter(as.numeric(distance) <= 400), shape = 1)
因此，蓝点是最终的阈值捕捉点。正如我们所看到的，400 m 缓冲区之外的点（红色多边形）保留在其原始位置，而内部的任何内容（黑点和虚线）都会捕捉到道路。
st_nearest_points()与@TimSalabim相比st_snap_points()
从时间上看，该st_nearest_points()方法似乎比st_snap_points(). 据我所知（这不是我的强项），这是因为st_snap_points()当在 上使用时，它是在逐点的基础上或逐行的基础上工作sf collection，而我们能够在整个一次使用collections 。LINESTRINGst_nearest_points()
对于 100 点的示例，我们讨论的是 <0.2 秒与 ~11 秒。
此外，如果我们尝试对 进行逐点 ( rowwise()) 方法st_nearest_points()，它是所有三种方法中最慢的，大约需要 15 秒。
#~~~~~ st_nearest_points() ~~~~~
system.time(
my_points %>%
mutate(
my_linestring = st_nearest_points(geometry, roads),
closest_point = st_cast(my_linestring, 'POINT')[seq(2, nrow(.)*2, 2)],
distance = st_distance(geometry, roads)[,1],
snapped_point_cond = st_sfc(ifelse(as.numeric(distance) <= 400, st_geometry(closest_point), geometry), crs = st_crs(roads))
) %>%
{. ->> closest_points_f1}
)
# <0.2 secs
#~~~~~ st_snap_points() ~~~~~
st_snap_points = function(x, y, max_dist = 1000) {
if (inherits(x, "sf")) n = nrow(x)
if (inherits(x, "sfc")) n = length(x)
out = do.call(c,
lapply(seq(n), function(i) {
nrst = st_nearest_points(st_geometry(x)[i], y)
nrst_len = st_length(nrst)
nrst_mn = which.min(nrst_len)
if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
return(st_cast(nrst[nrst_mn], "POINT")[2])
})
)
return(out)
}
system.time(
st_snap_points(my_points, roads, max_dist = 400) %>%
{. ->> closest_points_f2}
)
# 11 secs
#~~~~~ st_nearest_points() rowwise ~~~~~
system.time(
my_points %>%
rowwise %>%
mutate(
my_linestring = st_nearest_points(geometry, roads),
closest_point = st_cast(my_linestring, 'POINT')[2], # pull out the second point of each row
distance = st_distance(geometry, roads)[,1],
snapped_point_cond = st_sfc(ifelse(as.numeric(distance) <= 400, st_geometry(closest_point), geometry), crs = st_crs(roads))
) %>%
{. ->> closest_points_f1_rowwise}
)
# 15 secs
所有三种方法都会产生相同的最终数字：
ggplot()+
geom_sf(data = roads, col = 'red')+
geom_sf(data = closest_points_f1$snapped_point_cond, shape = 1, col = 'blue')
但是，该st_nearest_points()方法返回sf collection包含所有计算列的 ，而st_snap_points()仅生成最终的捕捉点。我认为拥有所有原始和工作列非常有用，特别是对于故障排除。另外，这种方法速度明显更快。
closest_points_f1
# Simple feature collection with 100 features and 2 fields
# Active geometry column: geometry
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 664090.7 ymin: 6524310 xmax: 665500.6 ymax: 6526107
# Projected CRS: WGS 84 / UTM zone 33N
# # A tibble: 100 x 6
# id geometry my_linestring closest_point distance snapped_point_cond
# * <int> <POINT [m]> <LINESTRING [m]> <POINT [m]> [m] <POINT [m]>
# 1 1 (664834.1 6525742) (664834.1 6525742, 665309 6525745) (665309 6525745) 475. (664834.1 6525742)
# 2 2 (665176.8 6524370) (665176.8 6524370, 665142.8 6524486) (665142.8 6524486) 121. (665142.8 6524486)
# 3 3 (664999.9 6525528) (664999.9 6525528, 665337.9 6525639) (665337.9 6525639) 356. (665337.9 6525639)
# 4 4 (665315.7 6525242) (665315.7 6525242, 665454.2 6525178) (665454.2 6525178) 153. (665454.2 6525178)
# 5 5 (664601 6525464) (664601 6525464, 665317.8 6525700) (665317.8 6525700) 755. (664601 6525464)
# 6 6 (665320.7 6525600) (665320.7 6525600, 665348.4 6525615) (665348.4 6525615) 31.2 (665348.4 6525615)
# 7 7 (664316.2 6525065) (664316.2 6525065, 664069.9 6524591) (664069.9 6524591) 534. (664316.2 6525065)
# 8 8 (665204 6526025) (665204 6526025, 665367.7 6526036) (665367.7 6526036) 164. (665367.7 6526036)
# 9 9 (664319.8 6524335) (664319.8 6524335, 664354.4 6524390) (664354.4 6524390) 64.5 (664354.4 6524390)
# 10 10 (664101.7 6525830) (664101.7 6525830, 665309 6525745) (665309 6525745) 1210. (664101.7 6525830)
# # ... with 90 more rows
closest_points_f2
# Geometry set for 100 features
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 664069.9 ymin: 6524380 xmax: 665480.3 ymax: 6526107
# Projected CRS: WGS 84 / UTM zone 33N
# First 5 geometries:
# POINT (664834.1 6525742)
# POINT (665142.8 6524486)
# POINT (665337.9 6525639)
# POINT (665454.2 6525178)
# POINT (664601 6525464)
我以 1000 点再次运行该示例，st_neareast_points()花费了约 0.31 秒（长了 1.5 倍），并st_snap_points()花费了 120 秒（长了约 10 倍）。我还跑st_nearest_points()了 10,000 点，花费了约 2.5 秒（比 100 点长 12 倍；比 1000 点长 8 倍）。我没有运行st_snap_points()10,000 行。
我不知道为什么st_snap返回线串的起点/终点。也许这是https://github.com/r-spatial/sf上的一个问题。
这是一种似乎可以确定正确点的解决方法。请注意，这st_nearest_points只是最近才引入的，因此您可能需要从github 安装sf。
nrst = st_nearest_points(point1, roads)
new_point2 = st_cast(nrst, "POINT")[2]
mapView(point1, color="red") + st_buffer(point1, dist = cut_dist) + new_point2 + roads
包装这是一个函数，用于返回新的几何集，其中的捕捉点低于某个特定值max_dist：
library(sf)
library(mapview)
st_snap_points = function(x, y, max_dist = 1000) {
if (inherits(x, "sf")) n = nrow(x)
if (inherits(x, "sfc")) n = length(x)
out = do.call(c,
lapply(seq(n), function(i) {
nrst = st_nearest_points(st_geometry(x)[i], y)
nrst_len = st_length(nrst)
nrst_mn = which.min(nrst_len)
if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
return(st_cast(nrst[nrst_mn], "POINT")[2])
})
)
return(out)
}
brew = st_transform(breweries, st_crs(trails))
tst = st_snap_points(brew, trails, 500)
mapview(breweries) + mapview(tst, col.regions = "red") + trails
- 这看起来好多了。谢谢。也许您想建议将此函数合并到 `sf` 库中 https://github.com/r-spatial/sf/issues/792 (2认同)