Векторизация для l oop - самодельное геокодирование - PullRequest
Новая
       6

Векторизация для l oop - самодельное геокодирование

0 голосов
/ 15 января 2020

Я анализирую свою историю местоположений в Google (сбрасывается с здесь , если кому-то интересно). Теперь в набор данных не входит ни одного поля, содержащего название города, но, учитывая, что для каждой строки есть комбинация широта / долгота, мы можем вычислить ее самостоятельно.
Учитывая, что мой набор данных имеет длину 1,2 млн строк, используя бесплатные API геокодирования вне таблицы (трафик c явно задушен).

Данные

Пара аэропортов 

airport_coords <-
  structure(
    list(
      V1 = c("LIMC", "LIRF"),
      V2 = c("MXP", "FCO"),
      V3 = c("MALPENSA", "FIUMICINO"),
      V4 = c("MILANO", "ROME"),
      V5 = c("ITALY", "ITALY"),
      V6 = c(45L, 41L),
      V7 = c(37L, 48L),
      V8 = c(53L, 46L),
      V9 = c("N", "N"),
      V10 = c(8L, 12L),
      V11 = c(43L, 15L),
      V12 = c(40L, 11L),
      V13 = c("E", "E"),
      V14 = c(234L, 4L),
      V15 = c(45.631, 41.813),
      V16 = c(8.728,
              12.253)
    ),
    row.names = c(NA,-2L),
    class = "data.frame"
  )

А вот несколько строк упрощенной версии история местоположений от Google 

loc_history <- 
  structure(list(latitudeGPS = c(41.8713521, 41.8713478, 41.8714064, 
41.8714201, 41.8713419, 41.8713981, 41.8713237, 41.8714538, 41.8713845, 
41.8714139, 41.8714417, 41.8714538, 41.8714417, 41.8714538, 41.8714538, 
41.8714538, 41.8714538, 41.8714538, 41.8714594, 41.8714594), 
    longitudeGPS = c(12.4414861, 12.441478, 12.4415342, 12.4415539, 
    12.4414757, 12.4415345, 12.4414538, 12.4415871, 12.441514, 
    12.4415466, 12.4415735, 12.4415871, 12.4415735, 12.4415871, 
    12.4415871, 12.4415871, 12.4415871, 12.4415871, 12.4415954, 
    12.4415954)), row.names = c(NA, 20L), class = "data.frame")

Лоскутное решение

Мой подход заключается в том, чтобы вычислить расстояние между координатами широты и долготы и аэропортом городов, которые меня интересуют (от это набор данных), предполагая, что если расстояние <50 км, я смотрю на город, где расположен аэропорт (который должен быть достаточно точным для моих нужд). Я написал следующее <code>for l oop (я знаю ...), которое работает, хотя и очень медленно. Я ищу способы превратить этот мусорный бак в нечто более быстрое, используя векторизованные функции, такие как семейство apply. 

library(raster)  # for pointDistance
library(dplyr)

# Init empty df to store results
dist <- data.frame(
  dist_mt = NA,
  city = NA
)
for (i in 1:nrow(loc_history)) {

  # Tmp df to store computed distances
  tmp <- data.frame(
    dist_mt = NA,
    city = NA
    )

  for (x in 1:nrow(airport_coords)) {
    # Coompute point - airport distance
    v <- pointDistance(c(data[i,]$latitudeGPS,
                         data[i,]$longitudeGPS),
                       c(airport_coords[x,]$V15,
                         airport_coords[x,]$V16),
                       lonlat = TRUE)

    # Append to tmp dataframe
    tmp[x,]$dist_mt <- v
    tmp[x,]$city <- airport_coords[x,]$V4  # Keep city label
  }

  # Append city if distance < 50km
  if (min(tmp$dist_mt) <= 50000) {
    dist[i,] <- filter(tmp, dist_mt == min(dist_mt))
  } else {
    dist[i,]$city <- "other"
  }

}

Performance

Для обработки l oop требуется около 4 секунд ~ 1.0000 строк. При наличии 1,2 млн. Строк для его запуска потребуется ~ 80 минут.

Ответы [ 2 ]

2 голосов
/ 15 января 2020

Попробуйте использовать пакеты sf и lwgeom: 

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
airport_coords = st_as_sf(airport_coords, coords=c('V16', 'V15'), crs=4326)

loc_history = st_as_sf(loc_history, coords=c('longitudeGPS', 'latitudeGPS'), crs=4326)

dist = st_distance(loc_history, airport_coords)
dist
#> Units: [m]
#>           [,1]     [,2]
#>  [1,] 513625.5 16943.33
#>  [2,] 513625.5 16942.53
#>  [3,] 513622.8 16949.33
#>  [4,] 513622.4 16951.42
#>  [5,] 513625.9 16942.10
#>  [6,] 513623.5 16949.00
#>  [7,] 513626.6 16939.65
#>  [8,] 513620.9 16955.40
#>  [9,] 513623.8 16946.85
#> [10,] 513622.6 16950.60
#> [11,] 513621.4 16953.84
#> [12,] 513620.9 16955.40
#> [13,] 513621.4 16953.84
#> [14,] 513620.9 16955.40
#> [15,] 513620.9 16955.40
#> [16,] 513620.9 16955.40
#> [17,] 513620.9 16955.40
#> [18,] 513620.9 16955.40
#> [19,] 513620.8 16956.27
#> [20,] 513620.8 16956.27

closest = apply(dist, 1, 
            function(r) ifelse(min(r)<=50000, airport_coords$V4[which.min(r)], NA))

Создано в 2020-01-15 пакетом prex (v0.3.0)

1 голос
/ 15 января 2020

Вам нужно создать матрицу из ваших данных, а не передавать одно значение за раз в соответствии с pointDistance справочным документом:

Аргументы

p1 координаты x и y первой (набора) точек (точек), либо c (x, y), матрицы (ncol = 2) или SpatialPoints *.

p2 x и координата y второй (набора) второй (ых) точки (точек) (как для p1). Если этот аргумент отсутствует, матрица расстояний вычисляется для p1

. Таким образом, чтобы получить все это за один проглот, вы должны сделать 

pointDistance(  
  matrix(c(loc_history$longitudeGPS, loc_history$latitudeGPS), ncol=2),
  matrix(c(airport_coords$V16, airport_coords$V15), ncol =2), 
  lonlat = TRUE) -> distmat

distmat
#>           [,1]     [,2]
#>  [1,] 513625.5 16943.33
#>  [2,] 513625.5 16942.53
#>  [3,] 513622.8 16949.33
#>  [4,] 513622.4 16951.42
#>  [5,] 513625.9 16942.10
#>  [6,] 513623.5 16949.00
#>  [7,] 513626.6 16939.65
#>  [8,] 513620.9 16955.40
#>  [9,] 513623.8 16946.85
#> [10,] 513622.6 16950.60
#> [11,] 513621.4 16953.84
#> [12,] 513620.9 16955.40
#> [13,] 513621.4 16953.84
#> [14,] 513620.9 16955.40
#> [15,] 513620.9 16955.40
#> [16,] 513620.9 16955.40
#> [17,] 513620.9 16955.40
#> [18,] 513620.9 16955.40
#> [19,] 513620.8 16956.27
#> [20,] 513620.8 16956.27

, так как каждый столбец представляет расстояние для каждого аэропорта (в порядке их появления в кадре данных вашего аэропорта), если вы найдете минимальное значение в каждой строке, вы найдете индекс аэропорта. Вы можете сделать это с помощью apply

loc_history$nearest_airport <- apply(distmat, 1, function(x) 
         { if(x[which.min(x)] < 50000) airport_coords$V4[which.min(x)] else NA })
loc_history$distance_to_nearest_airport <- apply(distmat, 1, min)

, и это должно быть результатом, который вы искали: 

loc_history
#>    latitudeGPS longitudeGPS nearest_airport distance_to_nearest_airport
#> 1     41.87135     12.44149            ROME                    16943.33
#> 2     41.87135     12.44148            ROME                    16942.53
#> 3     41.87141     12.44153            ROME                    16949.33
#> 4     41.87142     12.44155            ROME                    16951.42
#> 5     41.87134     12.44148            ROME                    16942.10
#> 6     41.87140     12.44153            ROME                    16949.00
#> 7     41.87132     12.44145            ROME                    16939.65
#> 8     41.87145     12.44159            ROME                    16955.40
#> 9     41.87138     12.44151            ROME                    16946.85
#> 10    41.87141     12.44155            ROME                    16950.60
#> 11    41.87144     12.44157            ROME                    16953.84
#> 12    41.87145     12.44159            ROME                    16955.40
#> 13    41.87144     12.44157            ROME                    16953.84
#> 14    41.87145     12.44159            ROME                    16955.40
#> 15    41.87145     12.44159            ROME                    16955.40
#> 16    41.87145     12.44159            ROME                    16955.40
#> 17    41.87145     12.44159            ROME                    16955.40
#> 18    41.87145     12.44159            ROME                    16955.40
#> 19    41.87146     12.44160            ROME                    16956.27
#> 20    41.87146     12.44160            ROME                    16956.27

Вы должны получить NA в столбце near_airport, если есть нет аэропорта в пределах 50 км.

Другими словами, вы можете заменить весь свой «пожарный контейнер» на: 

distmat <- pointDistance(  
  matrix(c(loc_history$longitudeGPS, loc_history$latitudeGPS), ncol=2),
  matrix(c(airport_coords$V16, airport_coords$V15), ncol =2), 
  lonlat = TRUE)

loc_history$nearest_airport <- apply(distmat, 1, function(x) 
         { if(x[which.min(x)] < 50000) airport_coords$V4[which.min(x)] else NA })

loc_history$distance_to_nearest_airport <- apply(distmat, 1, min)
...