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

Question

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

#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 来使点和多边形相交，但不确定如何在孤立点上运行类似的场景。

Answer 1

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

以下是如何使用 sf 包完成此操作。首先，在 lon-lat 坐标中创建点:

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 代码。接下来我们需要选择一个 map 投影。对于这个例子，我将使用 UTM zone 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

这表明前两个点在圆圈内，而第三个不在圆圈内。

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