在 for 循环内循环以计算 R 中 2 个数据集之间的地理空间距离

Question

我有一个带有 957 个地理编码的 data.table。我想将它与另一个具有 317 个地理编码的数据集匹配。匹配条件是地理空间接近度。我想将第一个数据集中的每个观察结果与第二个数据集中的每个观察结果相匹配，使得两个观察结果之间的距离不超过 5000 米。

我的数据是这样的：

> muni[1:3]
         mun Lat_Decimal Lon_Decimal
1:      1001    21.76672   -102.2818
2:      1002    22.16597   -102.0657
3:      1003    21.86138   -102.7248
> stations[1:3]
   station_number station_lat station_long
1:          10003      25.100     -106.567
2:          10018      24.944     -106.259
3:          10031      24.523     -105.952

我正在使用 library(geosphere) 中的 distm 函数来计算距离。

我认为解决这个问题的方法是 while 循环。这个想法是从 muni 中获取第一个观测值并测量到 stations 中第一个观测值的距离。如果距离小于或等于5000米，则将station中第一个观测值的station_number分配给muni中的第一个观测值。如果距离大于5000，则尝试muni中的下一次观察，直到距离小于或等于5000米。

本质上，它是一个循环，它在 stations 中找到第一个观测值，即 5000 米或更接近 muni 中的观测值。

这是初步尝试：

for (i in 1:957) {
  j = 1
  while (distm(muni[i, .(Lon_Decimal, Lat_Decimal)],
               stations[j, .(station_long, station_lat)]) > 5000 & j <= 317) {
    muni[i, station_number := as.integer(stations[j, station_number])]
    muni[i, distance := distm(muni[i, .(Lon_Decimal, Lat_Decimal)],
                                   stations[j, .(station_long, station_lat)])]
    j = j + 1
}
}

我可以看出这是行不通的，因为“muni”中的 none 行似乎在运行这个循环 for (i in 1:3) 之后被覆盖了。我想我的循环中有一个错误忽略了 station_number := 和 distance := 部分。

我希望这个循环覆盖 muni，这样整个列都有一个 station_number。

Answer 1

我已经将您的几个示例点读为 data.frames 并将它们转换为下面的 sf 以获得答案。如果您依附于 geosphere，请原谅双关语，考虑到 geosphere::distm 也是 returns 距离矩阵，一切都应该仍然适用。

首先我们将您的数据转换为 sf 格式：


library(sf)

stations_raw <- "station_number station_lat station_long
1:          10003      25.100     -106.567
2:          10018      24.944     -106.259
3:          10031      24.523     -105.952"


mun_raw <- "mun Lat_Decimal Lon_Decimal
1:      1001    21.76672   -102.2818
2:      1002    22.16597   -102.0657
3:      1003    21.86138   -102.7248"

mun_df <- read.table(text = mun_raw)

stations_df <- read.table(text = stations_raw)

mun_sf <- st_as_sf(mun_df, coords = c("Lon_Decimal", "Lat_Decimal"), crs = 4326)
stations_sf <-  st_as_sf(stations_df, 
                          coords = c("station_long", "station_lat"), crs = 4326)

然后，我们找到点之间每次交互的最小值：

closest <- list()

for(i in seq_len(nrow(mun_sf))){
  closest[[i]] <- stations_sf[which.min(
    st_distance(stations_sf, mun_sf[i,])),]
}

最后，我们提取标识符并将它们附加到原始 df，按照您的要求删除 mun_id：


mun_sf$closest_station <- purrr::map_chr(closest, "station_number")

mun_sf <- mun_sf[, c("closest_station", "geometry")]

mun_sf
#> Simple feature collection with 3 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -102.7248 ymin: 21.76672 xmax: -102.0657 ymax: 22.16597
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#>    closest_station                   geometry
#> 1:           10031 POINT (-102.2818 21.76672)
#> 2:           10031 POINT (-102.0657 22.16597)
#> 3:           10031 POINT (-102.7248 21.86138)

下图有助于直观地检查，在这个玩具示例中，我们得到了正确的答案。

ggplot() +
  geom_sf(data = mun_sf, colour = "red") +
  geom_sf_text(data = mun_sf, aes(label = mun), nudge_x = 0.25) +
  geom_sf(data = stations_sf, colour = "blue") +
  geom_sf_text(data = stations_sf, aes(label = station_number), nudge_x = -0.25)
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
#> not give correct results for longitude/latitude data

#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
#> not give correct results for longitude/latitude data

在 for 循环内循环以计算 R 中 2 个数据集之间的地理空间距离

While loop inside a for loop to calculate geospatial distance between 2 datasets in R

r

gis

sf

geosphere