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

У меня 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.


r sf gis geosphere
person Arturo Sbr    schedule 11.02.2019    source источник
comment
Вы были бы удовлетворены, если бы смогли сопоставить каждое muni наблюдение с его ближайшим station, а не первое наблюдение из station набора данных на глубине менее 5000 метров?   -  person Fons MA    schedule 11.02.2019
comment
Можете ли вы предоставить набор данных, с которым мы можем работать?   -  person Werner    schedule 11.02.2019
comment
@FonsMA Я подумал, что будет эффективнее запустить его на ‹= 5000 метров, но да, это было бы даже лучше.   -  person Arturo Sbr    schedule 11.02.2019

Ответы (1)


arrow_upward
1
arrow_downward

Я прочитал несколько ваших точек выборки как 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

person Fons MA    schedule 11.02.2019