循环并查找子组的最近数据点

Question

我正在尝试复制在真正的 cool nearest neighbor question 中完成的工作，但要针对我数据框中的每个区域而不是整个组进行复制。

我的数据ncbaby（别问）是这样的：

     id      printid     areaname latitude longitude
1  7912048 233502729     073    36.06241 -80.44229
2   735253 171241999 Area 12-06 35.54452 -78.75388
3  4325564  85564887 Area 12-04 35.49328 -78.73756
4  4222241  85461255 Area 12-06 35.53621 -78.75553
5 11997754 356053648 Area 12-04 35.49328 -78.73756
6 13444458 536073775 Area 12-06 35.53987 -78.74922

我想为每个区域名称运行该函数。我尝试调用 split 但距离函数不会调用列表。

splitfile <- split(ncbaby, ncbaby$precinctname)

c <- gDistance(splitfile, byid=TRUE)

Error in (function (classes, fdef, mtable)  : 
unable to find an inherited method for function ‘is.projected’ for signature ‘"list"’

我得到的最接近的是：

library(sp)
library(rgeos)

uniq <- unique(unlist(ncbaby$areaname))
for (i in 1:length(uniq)){
    data_1 <- subset(ncbaby, areaname == uniq[i])
    sp.mydata <- data_1
    coordinates(sp.mydata) <- ~longitude+latitude
    d <- gDistance(sp.mydata, byid=TRUE)
    min.d <- apply(d, 1, function(x) order(x, decreasing=F)[2])
    newdata <- cbind(data_1, data_1[min.d,], apply(d, 1, function(x) sort(x, decreasing=F)[2]))
    colnames(newdata) <- c(colnames(data_1), 'n.ID',
                          'n.printid', '.Areaname', 'n.latitude', 'n.longitude',
                          'distance')
}

这里的问题是它最终只踢出最后一个返回值。想法？我也有兴趣/愿意修改该功能，因为它似乎可能更有效。

Answer 1

我检查了链接的帖子并稍微修改了这个想法。我认为使用apply() 对于大型数据集可能不是一个好主意。所以我宁愿使用与 data.table 相关的方法。首先，我将示例数据转换为 SpatialPointsDataFrame。然后，我按组变量（即组）拆分数据。正如 Eddie 所建议的，我将lapply() 与 data.table 函数结合使用。当您使用gDistance() 时，您有一个二维向量作为输出。我将其转换为 data.table 对象，以便后续数据处理可能更快。我用melt() 重新塑造了dt 对象，并删除了距离= 0 的所有数据点。最后，我为每个Var1 取了第一行。请注意这里的Var1代表样本数据的每一行，mydf。最后一项工作是将新的距离向量添加到原始数据帧中。我希望这会对你有所帮助。

数据

   group user_id  latitude longitude
1    B23   85553 -34.44011  172.6954
2    B23   85553 -34.43929  172.6939
3    B23   85553 -34.43929  172.6939
4    B23   85553 -34.43851  172.6924
5    B23   57357 -34.42747  172.6778
6    B23   57357 -34.42747  172.6778
7    B23   57357 -34.42747  172.6778
8    B23   98418 -34.43119  172.7014
9    B23   98418 -34.43225  172.7023
10   B23   98418 -34.43224  172.7023
11   B23   98418 -34.43224  172.7023
12   B24   57357 -34.43647  172.7141
13   B24   57357 -34.43647  172.7141
14   B24   57357 -34.43647  172.7141
15   B24   98418 -34.43904  172.7172
16   B24   98418 -34.43904  172.7172
17   B24   98418 -34.43904  172.7172
18   B24   98418 -34.43925  172.7168
19   B24   98418 -34.43915  172.7169
20   B24   98418 -34.43915  172.7169
21   B24   98418 -34.43915  172.7169
22   B24   98418 -34.43915  172.7169

代码

library(sp)
library(rgeos)
library(data.table)

# Copy the original
temp <- mydf

#DF to SPDF
coordinates(temp) <- ~longitude+latitude

# Split the data by a group variable
mylist <- split(temp, f = temp$group)

#For each element in mylist, apply gDistance, reshape the output of
# gDistance and create a data.table. Then, reshape the data, remove
# rows with distance = 0. Finally, choose the first row for each 
# variable. levels in variable represents rows in mydf.

out <- rbindlist(
        lapply(mylist, function(x){

           d <- setDT(melt(gDistance(x, byid = TRUE)))
           setorder(d, Var1, value)
           d <- d[value > 0]
           d <- d[, .SD[1], by = Var1]
           d 

        })
    )

out <- cbind(mydf, distance = out$value)

#   group user_id  latitude longitude     distance
#1    B23   85553 -34.44011  172.6954 1.743945e-03
#2    B23   85553 -34.43929  172.6939 1.661118e-03
#3    B23   85553 -34.43929  172.6939 1.661118e-03
#4    B23   85553 -34.43851  172.6924 1.661118e-03
#5    B23   57357 -34.42747  172.6778 1.836642e-02
#6    B23   57357 -34.42747  172.6778 1.836642e-02
#7    B23   57357 -34.42747  172.6778 1.836642e-02
#8    B23   98418 -34.43119  172.7014 1.369100e-03
#9    B23   98418 -34.43225  172.7023 1.456022e-05
#10   B23   98418 -34.43224  172.7023 1.456022e-05
#11   B23   98418 -34.43224  172.7023 1.456022e-05
#12   B24   57357 -34.43647  172.7141 3.862696e-03
#13   B24   57357 -34.43647  172.7141 3.862696e-03
#14   B24   57357 -34.43647  172.7141 3.862696e-03
#15   B24   98418 -34.43904  172.7172 3.245705e-04
#16   B24   98418 -34.43904  172.7172 3.245705e-04
#17   B24   98418 -34.43904  172.7172 3.245705e-04
#18   B24   98418 -34.43925  172.7168 1.393162e-04
#19   B24   98418 -34.43915  172.7169 1.393162e-04
#20   B24   98418 -34.43915  172.7169 1.393162e-04
#21   B24   98418 -34.43915  172.7169 1.393162e-04
#22   B24   98418 -34.43915  172.7169 1.393162e-04

dput() 中的数据

mydf <- structure(list(group = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("B23", 
"B24"), class = "factor"), user_id = c(85553L, 85553L, 85553L, 
85553L, 57357L, 57357L, 57357L, 98418L, 98418L, 98418L, 98418L, 
57357L, 57357L, 57357L, 98418L, 98418L, 98418L, 98418L, 98418L, 
98418L, 98418L, 98418L), latitude = c(-34.440114, -34.43929, 
-34.43929, -34.438507, -34.427467, -34.427467, -34.427467, -34.431187, 
-34.432254, -34.43224, -34.43224, -34.436472, -34.436472, -34.436472, 
-34.439038, -34.439038, -34.439038, -34.439246, -34.439149, -34.439149, 
-34.439149, -34.439149), longitude = c(172.695443, 172.693906, 
172.693906, 172.692441, 172.677763, 172.677763, 172.677763, 172.701413, 
172.702284, 172.702288, 172.702288, 172.71411, 172.71411, 172.71411, 
172.717203, 172.717203, 172.717203, 172.716798, 172.716898, 172.716898, 
172.716898, 172.716898)), .Names = c("group", "user_id", "latitude", 
"longitude"), row.names = c(NA, -22L), class = "data.frame")

Answer 2

这是一个解决方案，可将链接帖子中的解决方案调整为按区域分隔组。首先，我们定义两个函数：

library(sp)
library(rgeos)

nearest.neighbor <- function(lon,lat) {
  df <- data.frame(lon,lat)
  coordinates(df) <- ~lon+lat
  d <- gDistance(df, byid=TRUE)
  # remove the self distance from being considered and use which.min to find the nearest neighbor
  d[cbind(1:nrow(d),1:nrow(d))] <- NA
  min.d <- rbind(apply(d,1,function(x) {ind <- which.min(x); list(ind=ind,distance=x[ind])}))
}

order.by.ind <- function (x,ind) x[ind]

nearest.neighbor 函数与链接帖子中的代码非常相似，只是它返回一个列表向量。每个列表都包含到最近邻居的索引和到该邻居的距离。这里的关键是我们只想计算一次距离以返回最小距离和相应的索引。请注意，我们通过将d 的对角线设置为NA，然后使用which.min 来定位最近的邻居，从而避免了必须进行完整排序，从而将自距离排除在外。

order.by.ind 函数只是根据索引ind 对输入列x 重新排序。

通过这两个函数，我们可以使用dplyr 包中的mutate 来计算按areaname 分组的所需列：

library(dplyr)

result <- ncbaby %>% group_by(areaname) %>%
                     mutate(min.d=nearest.neighbor(longitude, latitude)) %>%
                     mutate_each(vars=c(id, printid, latitude, longitude),
                                 funs(order.by.ind, "order.by.ind", order.by.ind(.,ind=unlist(min.d)[c(TRUE,FALSE)]))) %>%
                     mutate(distance=unlist(min.d)[c(FALSE,TRUE)]) %>%
                     mutate(.Areaname=areaname) %>%
                     select(-min.d)

newvars <- c('n.ID', 'n.printid', 'n.latitude', 'n.longitude', 'distance', '.Areaname')
colnames(result) <- c(colnames(ncbaby), newvars)

注意事项：

1. 第一个mutate 创建一个临时列min.d，其中包含最近邻的ind 和distance 列表。这是该地区最近的邻居，因为我们是group_by areaname。
2. 第二个mutate_each 为vars 中的每个变量创建一个新列，方法是根据ind 对该列重新排序。请注意，我们通过取消列出然后使用[c(TRUE,FALSE)] 提取奇数元素，从min.d 中提取ind。
3. 第三个mutate 通过从min.d 中提取distance 创建distance 列。同样，这是通过取消列出，然后使用 [c(FALSE,TRUE)] 提取偶数元素。
4. 实际上不需要第四个mutate，因为.Areaname 列对于areaname 在结果中是多余的。
5. 最后，我们从结果中删除min.d 列，并根据需要设置结果数据框的列名。

使用您的数据的结果是：

print(result)
##Source: local data frame [7 x 11]
##Groups: areaname [3]
##
##        id   printid   areaname latitude longitude     n.ID n.printid n.latitude n.longitude    distance  .Areaname
##     <int>     <int>     <fctr>    <dbl>     <dbl>    <int>     <int>      <dbl>       <dbl>       <dbl>     <fctr>
##1  7912048 233502729        073 36.06241 -80.44229  7912049 233502730   36.06251   -80.44329 0.001004988        073
##2   735253 171241999 Area 12-06 35.54452 -78.75388 13444458 536073775   35.53987   -78.74922 0.006583168 Area 12-06
##3  4325564  85564887 Area 12-04 35.49328 -78.73756 11997754 356053648   35.49328   -78.73756 0.000000000 Area 12-04
##4  4222241  85461255 Area 12-06 35.53621 -78.75553 13444458 536073775   35.53987   -78.74922 0.007294635 Area 12-06
##5 11997754 356053648 Area 12-04 35.49328 -78.73756  4325564  85564887   35.49328   -78.73756 0.000000000 Area 12-04
##6 13444458 536073775 Area 12-06 35.53987 -78.74922   735253 171241999   35.54452   -78.75388 0.006583168 Area 12-06
##7  7912049 233502730        073 36.06251 -80.44329  7912048 233502729   36.06241   -80.44229 0.001004988        073

我为areaname="073" 添加了一个新行，这样每个区域至少有两行。