Цикл while внутри цикла for для вычисления геопространственного расстояния между 2 наборами данных в R - PullRequest
Новая
       60

Цикл while внутри цикла for для вычисления геопространственного расстояния между 2 наборами данных в R

0 голосов
/ 11 февраля 2019

У меня есть data.table с 957 геокодами.Я хочу сопоставить его с другим набором данных с 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

Я использую функцию distm из library(geosphere) для вычисления расстояния.

Я решил, что для решения этой проблемы используется цикл while.Идея состоит в том, чтобы взять первое наблюдение из muni и измерить расстояние до первого наблюдения в stations.Если расстояние составляет 5000 метров или меньше, присвойте station_number первого наблюдения в station первому наблюдению в 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 ', кажется, небыли перезаписаны после выполнения этого цикла for (i in 1:3).Я предполагаю, что в моем цикле есть ошибка, которая игнорирует части station_number := и distance :=.

Я бы ожидал, что этот цикл перезапишет muni, так что весь столбец будет иметь station_number.

1 Ответ

0 голосов
/ 12 февраля 2019

Я прочитал ваши несколько выборочных точек как data.frames и преобразовал их в sf ниже для ответа.Если вы привязаны к geosphere, простите за каламбур, все должно применяться одинаково, учитывая, что geosphere::distm также возвращает матрицу расстояний.

Сначала мы получаем ваши данные в формате 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

...