从R中的已知中心点查找半径内的纬度/经度

Posted 2023-02-23

【中文标题】从R中的已知中心点查找半径内的纬度/经度

【英文标题】：Find lat/lon within a radius from a known centered point in R

【发布时间】：2017-08-19 19:38:47

【问题描述】：

我有一个已知纬度和经度的已知点（中心），并且我在 df 中有一些坐标（以纬度/经度为单位），我想看看哪些坐标在距离中心点 5 公里或更小的半径内。

中心

mylon <- c(-2.106472)
mylat <- c(57.14455)

坐标

longitude latitude
 [1,] -2.141177 57.16278
 [2,] -2.090960 57.18079
 [3,] -2.118894 57.12292
 [4,] -2.140090 57.13763
 [5,] -2.113988 57.13855
 [6,] -2.123892 57.13585
 [7,] -2.144685 57.17207
 [8,] -2.220046 57.19150
 [9,] -2.114343 57.09301
[10,] -2.285314 57.20138
[11,] -2.092354 57.14279
[12,] -2.149571 57.15334
[13,] -2.126605 57.15615
[14,] -2.097045 57.10443
[15,] -2.183441 57.15051
[16,] -2.166915 57.15089
[17,] -2.133863 57.15201
[18,] -2.100909 57.18968
[19,] -2.106770 57.15670
[20,] -2.251495 57.19315
[21,] -2.118894 57.12292
[22,] -2.140090 57.13763
[23,] -2.123201 57.12686
[24,] -2.114343 57.09301
[25,] -2.140327 57.15676
[26,] -2.148826 57.17355
[27,] -2.120553 57.12507
[28,] -2.133902 57.16279
[29,] -2.094246 57.17180
[30,] -2.113170 57.14125
[31,] -2.251495 57.19315
[32,] -2.090960 57.18079
[33,] -2.212955 57.10941
[34,] -2.118894 57.12292
[35,] -2.183501 57.19596
[36,] -2.140090 57.13763
[37,] -2.249217 57.10063
[38,] -2.123201 57.12686
[39,] -2.114343 57.09301
[40,] -2.092354 57.14279
[41,] -2.148826 57.17355
[42,] -2.120553 57.12507
[43,] -2.117338 57.15301
[44,] -2.116486 57.14484
[45,] -2.094981 57.13614
[46,] -2.232998 57.14629
[47,] -2.118894 57.12292
[48,] -2.140090 57.13763
[49,] -2.123201 57.12686
[50,] -2.104092 57.14485
[51,] -2.114343 57.09301
[52,] -2.148826 57.17355
[53,] -2.175179 57.15079
[54,] -2.090713 57.14755
[55,] -2.090960 57.18079
[56,] -2.118894 57.12292
[57,] -2.140090 57.13763

如果有任何帮助，我将不胜感激。非常感谢

【问题讨论】：

geosphere::distHaversine(coord, c(-2.106472, 57.14455)) / 1000 < 5

我认为这些点是以度为单位的，Haversine 公式以弧度为单位的点作为参数。那么我们不需要先将度数转换为弧度数吗？

@RichPauloo 不，它需要经度和纬度，它们几乎总是以度数表示。地图投影是一个问题，但这超出了问题的范围。

这两个答案都以自己的方式很好。我将使用最短的一行代码。干杯

【参考方案1】：

我们可以使用sf 包中的st_distance 函数来完成这个任务。 sf 是 R 中用于处理空间数据的下一代数据类型。下面的例子展示了如何创建两个sf对象，point_sf：参考，target_sf：检查目标点。之后，很容易应用st_distance 函数。 target_sf2 为最终输出，包含 44 条参考点 5​​ 公里范围内的记录。

library(tidyverse)
library(sf)

## Create reference point data
point <- data_frame(mylon = -2.106472, mylat = 57.14455) 
# Specify the source of X and Y coordinates
point_sf <- st_as_sf(point, coords = c("mylon", "mylat"))
# Set the projection to EPSG 4326 (long-lat)
st_crs(point_sf) <- 4326

## Create target point data
target <- read.table(text = "longitude latitude
-2.141177 57.16278
-2.090960 57.18079
-2.118894 57.12292
-2.140090 57.13763
-2.113988 57.13855
-2.123892 57.13585
-2.144685 57.17207
-2.220046 57.19150
-2.114343 57.09301
-2.285314 57.20138
-2.092354 57.14279
-2.149571 57.15334
-2.126605 57.15615
-2.097045 57.10443
-2.183441 57.15051
-2.166915 57.15089
-2.133863 57.15201
-2.100909 57.18968
-2.106770 57.15670
-2.251495 57.19315
-2.118894 57.12292
-2.140090 57.13763
-2.123201 57.12686
-2.114343 57.09301
-2.140327 57.15676
-2.148826 57.17355
-2.120553 57.12507
-2.133902 57.16279
-2.094246 57.17180
-2.113170 57.14125
-2.251495 57.19315
-2.090960 57.18079
-2.212955 57.10941
-2.118894 57.12292
-2.183501 57.19596
-2.140090 57.13763
-2.249217 57.10063
-2.123201 57.12686
-2.114343 57.09301
-2.092354 57.14279
-2.148826 57.17355
-2.120553 57.12507
-2.117338 57.15301
-2.116486 57.14484
-2.094981 57.13614
-2.232998 57.14629
-2.118894 57.12292
-2.140090 57.13763
-2.123201 57.12686
-2.104092 57.14485
-2.114343 57.09301
-2.148826 57.17355
-2.175179 57.15079
-2.090713 57.14755
-2.090960 57.18079
-2.118894 57.12292
-2.140090 57.13763",
                     header = TRUE, stringsAsFactors = FALSE)

# Specify the source of X and Y coordinates
target_sf <- st_as_sf(target, coords = c("longitude", "latitude"))
# Set the projection to EPSG 4326 (long-lat)
st_crs(target_sf) <- 4326

# Calculate the distance
target_sf2 <- target_sf %>%
  mutate(Dist = as.numeric(st_distance(point_sf, target_sf, by_element = TRUE))) %>%
  # Filter the records with Dist <= 5000 (m)
  filter(Dist <= 5000)

【参考方案2】：

这可以通过距离公式和转换来完成，但是一旦距离开始变大，精度就会变得很差。可以编写自己的更严格的方法，但使用 geosphere 包要容易得多，它提供了各种功能来确定地理空间距离，具有不同程度的准确度和所需的计算能力。 distHaversine 是一个很好的起点：

coord <- cbind("longitude" = c(-2.141177, -2.09096, -2.118894, -2.14009, -2.113988, -2.123892, -2.144685, -2.220046, -2.114343, -2.285314, -2.092354, -2.149571, -2.126605, -2.097045, -2.183441, -2.166915, -2.133863, -2.100909, -2.10677, -2.251495, -2.118894, -2.14009, -2.123201, -2.114343, -2.140327, -2.148826, -2.120553, -2.133902, -2.094246, -2.11317, -2.251495, -2.09096, -2.212955, -2.118894, -2.183501, -2.14009, -2.249217, -2.123201, -2.114343, -2.092354, -2.148826, -2.120553, -2.117338, -2.116486, -2.094981, -2.232998, -2.118894, -2.14009, -2.123201, -2.104092, -2.114343, -2.148826, -2.175179, -2.090713, -2.09096, -2.118894, -2.14009), 
               "latitude" = c(57.16278, 57.18079, 57.12292, 57.13763, 57.13855, 57.13585, 57.17207, 57.1915, 57.09301, 57.20138, 57.14279, 57.15334, 57.15615, 57.10443, 57.15051, 57.15089, 57.15201, 57.18968, 57.1567, 57.19315, 57.12292, 57.13763, 57.12686, 57.09301, 57.15676, 57.17355, 57.12507, 57.16279, 57.1718, 57.14125, 57.19315, 57.18079, 57.10941, 57.12292, 57.19596, 57.13763, 57.10063, 57.12686, 57.09301, 57.14279, 57.17355, 57.12507, 57.15301, 57.14484, 57.13614, 57.14629, 57.12292, 57.13763, 57.12686, 57.14485, 57.09301, 57.17355, 57.15079, 57.14755, 57.18079, 57.12292, 57.13763))

str(coord)    # like data above, a matrix, not a data.frame
#>  num [1:57, 1:2] -2.14 -2.09 -2.12 -2.14 -2.11 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : NULL
#>   ..$ : chr [1:2] "longitude" "latitude"

# make a data.frame to hold both numeric and logical values
coord_df <- data.frame(coord, 
                       within_5km = geosphere::distHaversine(
                           coord, 
                           c(-2.106472, 57.14455)
                       ) / 1000 < 5)    # convert m to km, check < 5
str(coord_df)
#> 'data.frame':    57 obs. of  3 variables:
#>  $ longitude : num  -2.14 -2.09 -2.12 -2.14 -2.11 ...
#>  $ latitude  : num  57.2 57.2 57.1 57.1 57.1 ...
#>  $ within_5km: logi  TRUE TRUE TRUE TRUE TRUE TRUE ...

table(coord_df$within_5km)
#> 
#> FALSE  TRUE 
#>    13    44

【参考方案3】：

在 Google Maps API 上的this link 上，您可以找到一些解释和一些 SQL 代码来查找 n 英里范围内的内容。

要在 5 公里内使用它，您只需将公里转换为英里。所以你最终会得到：

SELECT id, 
    ( 3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * sin( radians( lat ) ) ) ) AS distance 
FROM markers 
HAVING distance < (5 /*km*/ / 0.621371192 /*convert to miles*/) 
ORDER BY distance

【讨论】：

这被标记为r。

你能告诉我这些数字代表什么，3959, 3, ...等