R - 查找给定半径内最近的相邻点和相邻点的数量，坐标经纬度答案

【问题标题】：R - Finding closest neighboring point and number of neighbors within a given radius, coordinates lat-longR - 查找给定半径内最近的相邻点和相邻点的数量，坐标经纬度
【发布时间】：2014-03-25 12:33:00
【问题描述】：

我试图弄清楚我的数据集中某些点的孤立程度。我使用两种方法来确定隔离，最近邻居的距离和给定半径内的相邻站点的数量。我所有的坐标都是经纬度

这是我的数据的样子：

    pond            lat         long        area    canopy  avg.depth   neighbor    n.lat   n.long  n.distance  n.area  n.canopy    n.depth n.avg.depth radius1500
    A10             41.95928    -72.14605   1500    66      60.61538462                                 
    AA006           41.96431    -72.121     250     0       57.77777778                                 
    Blacksmith      41.95508    -72.123803  361     77      71.3125                                 
    Borrow.Pit.1    41.95601    -72.15419   0       0       41.44444444                                 
    Borrow.Pit.2    41.95571    -72.15413   0       0       37.7                                    
    Borrow.Pit.3    41.95546    -72.15375   0       0       29.22222222                                 
    Boulder         41.918223   -72.14978   1392    98      43.53333333

我想把最近的相邻池塘的名字放在neighbor这个列中，它的lat和long在n.lat和n.long中，两个池塘之间的距离在n.distance中，以及面积，canopy和avg .depth 在每个适当的列中。

其次，我想把目标池塘1500m内的池塘数量放到radius1500中。

有谁知道可以帮助我计算我想要的距离/数字的函数或包？如果有问题，输入我需要的其他数据并不难，但最近邻居的名称和距离，加上1500m内的池塘数量是我真正需要帮助的。

谢谢。

【问题讨论】：

标签： r distance latitude-longitude

【解决方案1】：

最好的选择是使用库 sp 和 rgeos，它们使您能够构建空间类并执行地理处理。

library(sp)
library(rgeos)

读取数据并将其转换为空间对象：

mydata <- read.delim('d:/temp/testfile.txt', header=T)

sp.mydata <- mydata
coordinates(sp.mydata) <- ~long+lat

class(sp.mydata)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

现在计算点之间的成对距离

d <- gDistance(sp.mydata, byid=T)

求第二最短的距离（最近的距离是指向自身的，因此使用第二最短的）

min.d <- apply(d, 1, function(x) order(x, decreasing=F)[2])

用所需的变量构造新的数据框

newdata <- cbind(mydata, mydata[min.d,], apply(d, 1, function(x) sort(x, decreasing=F)[2]))

colnames(newdata) <- c(colnames(mydata), 'neighbor', 'n.lat', 'n.long', 'n.area', 'n.canopy', 'n.avg.depth', 'distance')

newdata
            pond      lat      long area canopy avg.depth     neighbor    n.lat    n.long n.area n.canopy n.avg.depth
6            A10 41.95928 -72.14605 1500     66  60.61538 Borrow.Pit.3 41.95546 -72.15375      0        0    29.22222
3          AA006 41.96431 -72.12100  250      0  57.77778   Blacksmith 41.95508 -72.12380    361       77    71.31250
2     Blacksmith 41.95508 -72.12380  361     77  71.31250        AA006 41.96431 -72.12100    250        0    57.77778
5   Borrow.Pit.1 41.95601 -72.15419    0      0  41.44444 Borrow.Pit.2 41.95571 -72.15413      0        0    37.70000
4   Borrow.Pit.2 41.95571 -72.15413    0      0  37.70000 Borrow.Pit.1 41.95601 -72.15419      0        0    41.44444
5.1 Borrow.Pit.3 41.95546 -72.15375    0      0  29.22222 Borrow.Pit.2 41.95571 -72.15413      0        0    37.70000
6.1      Boulder 41.91822 -72.14978 1392     98  43.53333 Borrow.Pit.3 41.95546 -72.15375      0        0    29.22222
        distance
6   0.0085954872
3   0.0096462277
2   0.0096462277
5   0.0003059412
4   0.0003059412
5.1 0.0004548626
6.1 0.0374480316

编辑：如果坐标以度为单位，并且您想以公里为单位计算距离，请使用包geosphere

library(geosphere)

d <- distm(sp.mydata)

# rest is the same

如果点分散在全球各地并且坐标以度为单位，这应该会提供更好的结果

【讨论】：

非常感谢。你推荐的库正是我需要的！
使用以下函数：sort(x[x>0], decreasing=F)[1]
@ike; sort 将行/列中的值从小到大排序，反之亦然。但由于您想省略零距离，首先您必须过滤数据 - x[x>0]。然后对它们进行排序，最后只取排序数组 ([1]) 中的第一个值。都清楚了吗？
距离以公里计算？
@NicoCoallier 它取决于坐标系的单位（对于 WGS，例如，它只是一个数字，它不使用半正弦距离）

【解决方案2】：

@Zbynek 提出的解决方案非常好，但是如果您像我一样正在寻找以公里为单位的两个邻居之间的距离，我提出了这个解决方案。

   earth.dist<-function(lat1,long1,lat2,long2){

           rad <- pi/180
           a1 <- lat1 * rad
           a2 <- long1 * rad
           b1 <- lat2 * rad
           b2 <- long2 * rad
           dlat <- b1-a1
           dlon<- b2-a2
           a <- (sin(dlat/2))^2 +cos(a1)*cos(b1)*(sin(dlon/2))^2
           c <- 2*atan2(sqrt(a),sqrt(1-a))
           R <- 6378.145
           dist <- R *c
           return(dist)
           }


    Dist <- matrix(0,ncol=length(mydata),nrow=length(mydata.sp))

  for (i in 1:length(mydata)){
      for(j in 1:length(mydata.sp)){
          Dist[i,j] <- earth.dist(mydata$lat[i],mydata$long[i],mydata.sp$lat[j],mydata.sp$long[j])
 }}



     DDD <- matrix(0, ncol=5,nrow=ncol(Dist))   ### RECTIFY the nb of col by the number of variable you want

   for(i in 1:ncol(Dist)){
       sub<- sort(Dist[,i])[2]
       DDD[i,1] <- names(sub) 
       DDD[i,2] <- sub
       DDD[i,3] <- rownames(Dist)[i]
       sub_neig_atr <- Coord[Coord$ID==names(sub),]
       DDD[i,4] <- sub_neig_atr$area
       DDD[i,5] <- sub_neig_atr$canopy
       ### Your can add any variable you want here 

   }

    DDD <- as.data.frame(DDD)

    names(DDD)<-c("neigboor_ID","distance","pond","n.area","n.canopy")
   data <- merge(mydata,DDD, by="pond")

如果您的坐标又长又宽，那么您最终会得到以公里为单位的距离。

有什么改进的建议吗？

【讨论】：

不用自己装，已经有geosphere包-cran.r-project.org/web/packages/geosphere/geosphere.pdf
该包中的哪个函数可以计算以公里为单位的欧几里得距离？
我认为distm 你可以选择精确的公式 - 默认是Haversine，但还有更多选项（参见手册）

【解决方案3】：

我在下面添加了一个使用 spatialrisk 包的解决方案。此包中的关键函数是用 C++ (Rcpp) 编写的，因此速度非常快。

首先，加载数据：

df <- data.frame(pond = c("A10", "AA006", "Blacksmith", "Borrow.Pit.1", 
                          "Borrow.Pit.2", "Borrow.Pit.3", "Boulder"), 
                 lat = c(41.95928, 41.96431, 41.95508, 41.95601, 
                         41.95571, 41.95546, 41.918223), 
                 long = c(-72.14605, -72.121, -72.123803, -72.15419, 
                          -72.15413, -72.15375, -72.14978), 
                 area = c(1500, 250, 361, 0, 0, 0, 1392), 
                 canopy = c(66, 0, 77, 0, 0, 0, 98), 
                 avg.depth = c(60.61538462, 57.77777778, 71.3125, 41.44444444,
                               37.7, 29.22222222, 43.53333333))

函数spatialrisk::points_in_circle 计算中心点半径内的观测值。请注意，距离是使用 Haversine 公式计算的。由于输出的每个元素都是一个数据框，purrr::map_dfr 用于将它们行绑定在一起：

ans1 <- purrr::map2_dfr(df$long, df$lat, 
                        ~spatialrisk::points_in_circle(df, .x, .y, 
                                                       lon = long, 
                                                       radius = 100000)[2,])

colnames(ans1) <- c("neighbor", "n.lat", "n.long", "n.area", 
                    "n.canopy", "n.avg.depth", "distance_m")

      neighbor    n.lat    n.long n.area n.canopy n.avg.depth distance_m
1 Borrow.Pit.1 41.95601 -72.15419      0        0    41.44444  765.87823
2   Blacksmith 41.95508 -72.12380    361       77    71.31250 1053.35200
3        AA006 41.96431 -72.12100    250        0    57.77778 1053.35200
4 Borrow.Pit.2 41.95571 -72.15413      0        0    37.70000   33.76321
5 Borrow.Pit.1 41.95601 -72.15419      0        0    41.44444   33.76321
6 Borrow.Pit.2 41.95571 -72.15413      0        0    37.70000   42.00128
7 Borrow.Pit.3 41.95546 -72.15375      0        0    29.22222 4158.21978

现在计算目标池塘 1500m 范围内的池塘数量。函数spatialrisk::concentration 将中心点半径内的观测数相加。从池塘数量中减去 1 以排除池塘本身。

df$npond <- 1  
radius1500 <- spatialrisk::concentration(df, df, npond, lon_sub = long, 
                                         lon_full = long, radius = 1500, 
                                         display_progress = FALSE)$concentration - 1

将数据框列绑定在一起：

cbind(df, ans1, radius1500)

          pond      lat      long area canopy avg.depth     neighbor    n.lat    n.long n.area n.canopy n.avg.depth distance_m radius1500
1          A10 41.95928 -72.14605 1500     66  60.61538 Borrow.Pit.1 41.95601 -72.15419      0        0    41.44444  765.87823          3
2        AA006 41.96431 -72.12100  250      0  57.77778   Blacksmith 41.95508 -72.12380    361       77    71.31250 1053.35200          1
3   Blacksmith 41.95508 -72.12380  361     77  71.31250        AA006 41.96431 -72.12100    250        0    57.77778 1053.35200          1
4 Borrow.Pit.1 41.95601 -72.15419    0      0  41.44444 Borrow.Pit.2 41.95571 -72.15413      0        0    37.70000   33.76321          3
5 Borrow.Pit.2 41.95571 -72.15413    0      0  37.70000 Borrow.Pit.1 41.95601 -72.15419      0        0    41.44444   33.76321          3
6 Borrow.Pit.3 41.95546 -72.15375    0      0  29.22222 Borrow.Pit.2 41.95571 -72.15413      0        0    37.70000   42.00128          3
7      Boulder 41.91822 -72.14978 1392     98  43.53333 Borrow.Pit.3 41.95546 -72.15375      0        0    29.22222 4158.21978          0

【讨论】：

【解决方案4】：

我在下面添加了一个使用较新的sf 包的替代解决方案，供那些感兴趣并现在访问此页面的人（就像我一样）。

首先，加载数据并创建sf 对象。

# Using sf
mydata <- structure(
  list(pond = c("A10", "AA006", "Blacksmith", "Borrow.Pit.1", 
                "Borrow.Pit.2", "Borrow.Pit.3", "Boulder"), 
       lat = c(41.95928, 41.96431, 41.95508, 41.95601, 41.95571, 41.95546, 
               41.918223), 
       long = c(-72.14605, -72.121, -72.123803, -72.15419, -72.15413, 
                -72.15375, -72.14978), 
       area = c(1500L, 250L, 361L, 0L, 0L, 0L, 1392L), 
       canopy = c(66L, 0L, 77L, 0L, 0L, 0L, 98L), 
       avg.depth = c(60.61538462, 57.77777778, 71.3125, 41.44444444, 
                     37.7, 29.22222222, 43.53333333)), 
  class = "data.frame", row.names = c(NA, -7L))


library(sf)
data_sf <- st_as_sf(mydata, coords = c("long", "lat"),
                    # Change to your CRS
                    crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
st_is_longlat(data_sf)

sf::st_distance 在使用纬度/经度数据时使用大圆距离计算以米为单位的距离矩阵。

dist.mat <- st_distance(data_sf) # Great Circle distance since in lat/lon
# Number within 1.5km: Subtract 1 to exclude the point itself
num.1500 <- apply(dist.mat, 1, function(x) {
  sum(x < 1500) - 1
})

# Calculate nearest distance
nn.dist <- apply(dist.mat, 1, function(x) {
  return(sort(x, partial = 2)[2])
})
# Get index for nearest distance
nn.index <- apply(dist.mat, 1, function(x) { order(x, decreasing=F)[2] })

n.data <- mydata
colnames(n.data)[1] <- "neighbor"
colnames(n.data)[2:ncol(n.data)] <- 
  paste0("n.", colnames(n.data)[2:ncol(n.data)])
mydata2 <- data.frame(mydata,
                      n.data[nn.index, ],
                      n.distance = nn.dist,
                      radius1500 = num.1500)
rownames(mydata2) <- seq(nrow(mydata2))

mydata2
          pond      lat      long area canopy avg.depth     neighbor    n.lat    n.long n.area n.canopy
1          A10 41.95928 -72.14605 1500     66  60.61538 Borrow.Pit.1 41.95601 -72.15419      0        0
2        AA006 41.96431 -72.12100  250      0  57.77778   Blacksmith 41.95508 -72.12380    361       77
3   Blacksmith 41.95508 -72.12380  361     77  71.31250        AA006 41.96431 -72.12100    250        0
4 Borrow.Pit.1 41.95601 -72.15419    0      0  41.44444 Borrow.Pit.2 41.95571 -72.15413      0        0
5 Borrow.Pit.2 41.95571 -72.15413    0      0  37.70000 Borrow.Pit.1 41.95601 -72.15419      0        0
6 Borrow.Pit.3 41.95546 -72.15375    0      0  29.22222 Borrow.Pit.2 41.95571 -72.15413      0        0
7      Boulder 41.91822 -72.14978 1392     98  43.53333 Borrow.Pit.3 41.95546 -72.15375      0        0
  n.avg.depth n.distance radius1500
1    41.44444  766.38426          3
2    71.31250 1051.20527          1
3    57.77778 1051.20527          1
4    37.70000   33.69099          3
5    41.44444   33.69099          3
6    37.70000   41.99576          3
7    29.22222 4149.07406          0

为了在计算距离后获得最近的邻居，您可以使用sort() 和partial = 2 参数。根据数据量的不同，这可能比使用order 快得多，如在之前的解决方案中。 Rfast 包可能更快，但我避免在此处包含其他包。有关各种解决方案的讨论和基准测试，请参阅此相关帖子：https://stackoverflow.com/a/53144760/12265198

【讨论】：

【解决方案5】：

在 Rfast 中， 有一个名为 “dista” 的函数，仅计算欧几里得或曼哈顿距离（目前）。它提供了计算 k 最小距离的选项。或者，它可以返回具有最小距离的观测值的索引。余弦距离与欧几里德距离基本相同（不包括常数，我认为是 2）。

【讨论】：

这方面的例子有机会吗？

【解决方案6】：

另一个可能较慢的答案可能对 dplyr 上瘾者具有直观的吸引力。

您创建一个包含纬度/经度的每个可能组合的巨型网格，然后您可以使用 geosphere 找到距离最短的那个。

该示例是您有两个具有不同点的数据集要比较 - 但您可以通过复制第一个数据集轻松调整它。

    library(tidyverse)
    library(geosphere)
    library(data.table)

    #This function creates a big dataframe with every possible combination
    expand.grid.df <- function(...) Reduce(function(...) merge(..., by=NULL), list(...))

shortest_distance <- expand.grid.df(df1,df2) %>%
      mutate(distance = distHaversine(p1 = cbind(lon_2,lat_2),
                                      p2 = cbind(lon,lat))) %>% 
      group_by(ACCIDENT_NO) %>% 
      slice(which.min(distance))

【讨论】：