R - 如何围绕一个点绘制半径并使用该结果过滤其他点？答案

【问题标题】：R - How do I draw a radius around a point and use that result to filter other points?R - 如何围绕一个点绘制半径并使用该结果过滤其他点？
【发布时间】：2020-10-29 21:24:39
【问题描述】：

我希望围绕一个经纬度点绘制一个半径，然后使用该缓冲区过滤适合其中的其他点。例如：

#stores datasets
stores = data.frame(store_id = 1:3,
                    lat = c("40.7505","40.7502","40.6045"),
                    long = c("-73.8456","-73.8453","-73.8012")
                    )

#my location
me  = data.frame(lat = "40.7504", long = "-73.8456")

#draw a 100 meter radius around me 


#use the above result to check which points in dataset stores are within that buffer

不知道如何解决这个问题。我之前曾与over 合作过点和多边形相交，但不知道如何在孤点上运行类似的场景。

【问题讨论】：

测量所有点到给定点的距离（使用毕达哥拉斯定理），然后选择距离小于所需阈值的点。
有各种提供圆的包，但实际上可以简单地生成角度 (seq(0,2*pi,len=51)) 并使用 sin 和 cos 将其转换为 x 的向量和ys。但是，您可能会遇到问题，因为其中大部分是基于有数字的，而您的数据中有字符串。另一个复杂因素：在纬度/经度坐标中，固定半径的圆在技术上更加困难（如方位/距离计算），因为毕达哥拉斯距离仅在小范围内近似正确；否则，事情就会开始崩溃（尤其是在 40 度和更多的北方）。
@Rodrigo，并不是那么简单，因为经纬度坐标不是平面坐标。您不能只将勾股定理应用于角度并期望它全部解决。
@clbieganek 是的，你是对的，它只是一个近似值。更好的选择是以米为单位计算距离，将地球在每个纬度的周长视为赤道周长乘以该纬度的余弦。然后使用每个距离的平均纬度作为测量该距离（以米为单位）的基础。这将提供更好的近似值。
你应该试试 geosphere 包中的 destPoint，它的设计目的是在 lon/lat 坐标上绘制圆圈 stackoverflow.com/a/54638040/7877917

标签： r spatial sp sf

【解决方案1】：

您可以假设尝试在球体或椭圆体的表面上进行几何计算，但通常在执行几何地图操作时，会使用地图投影，将经纬度坐标投影到平面上。

以下是使用 sf 包的方法。首先，在经纬度坐标中创建您的点：

library(sf)

lat <- c(40.7505, 40.7502, 40.6045)
lon <- c(-73.8456, -73.8453, -73.8012)

stores <- st_sfc(st_multipoint(cbind(lon, lat)), crs = 4326)

me <- st_sfc(st_point(c(-73.8456, 40.7504)), crs = 4326)

crs = 4326 参数指定经纬度坐标系的 EPSG 代码。接下来我们需要选择一个地图投影。在本例中，我将使用 UTM 区域 18，其中包含以上几点：

stores_utm <- st_transform(stores, "+proj=utm +zone=18")
me_utm     <- st_transform(me, "+proj=utm +zone=18")

现在我们可以将代表自己的点缓冲 100 米，生成一个半径为 100 米的圆：

circle <- st_buffer(me_utm, 100)

现在，我们几乎可以使用几何谓词来测试圆中的点了。但是，stores_utm 当前是MULTIPOINT，因此几何谓词会将其视为一个几何实体。我们可以通过将stores_utm 转换为POINT 来解决此问题，这将为我们提供三个不同点的集合：

stores_utm_column <- st_cast(stores_utm, "POINT")
stores_utm_column
# Geometry set for 3 features 
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 597453 ymin: 4495545 xmax: 601422.3 ymax: 4511702
# epsg (SRID):    32618
# proj4string:    +proj=utm +zone=18 +ellps=WGS84 +units=m +no_defs
# POINT (597453 4511702)
# POINT (597478.7 4511669)
# POINT (601422.3 4495545)

现在我们可以测试圆圈中有哪些点：

> st_contains(circle, stores_utm_column, sparse = FALSE)
#      [,1] [,2]  [,3]
# [1,] TRUE TRUE FALSE

表示前两点在圆内，第三点不在。

当然，每个地图投影都会引入一些失真。您选择的投影取决于您的问题的性质。

【讨论】：

【解决方案2】：

spatialrisk 包中的 points_in_circle() 函数处理这个问题。

例如，使用您的数据：

library(spatialrisk)

# Stores 
stores <- data.frame(store_id = 1:3,
                     lat = c(40.7505, 40.7502, 40.6045),
                     long = c(-73.8456, -73.8453, -73.8012))

# My location
me <- data.frame(lat = 40.7504, long = -73.8456)

> spatialrisk::points_in_circle(stores, me$long[1], me$lat[1], radius = 100, lon = long)
# store_id     lat     long distance_m
#        1 40.7505 -73.8456   11.13195
#        2 40.7502 -73.8453   33.70076

【讨论】：