Более быстрый способ оценки точек внутри списка полигонов в R

У меня есть два набора данных, один из которых содержит более 13 миллионов прямоугольных полигонов (набор из 4 точек широты и долготы), а другой — 10 тысяч точек, относящихся к ценам в этом месте.

> polygons
     id                                 pol_lat                                 pol_lng
 1: 148 -4.250236,-4.250236,-4.254640,-4.254640 -49.94628,-49.94494,-49.94494,-49.94628
 2: 149 -4.254640,-4.254640,-5.361601,-5.361601 -49.94494,-49.07906,-49.07906,-49.94494
 3: 150 -5.361601,-5.361601,-5.212208,-5.212208 -49.07906,-49.04469,-49.04469,-49.07906
 4: 151 -5.212208,-5.212208,-5.002878,-5.002878 -49.04469,-48.48664,-48.48664,-49.04469
 5: 152 -5.002878,-5.002878,-5.080018,-5.080018 -48.48664,-48.43699,-48.43699,-48.48664
 6: 153 -5.080018,-5.080018,-5.079819,-5.079819 -48.43699,-48.42480,-48.42480,-48.43699
 7: 154 -5.079819,-5.079819,-5.155606,-5.155606 -48.42480,-47.53891,-47.53891,-48.42480
 8: 155 -5.155606,-5.155606,-4.954156,-4.954156 -47.53891,-47.50354,-47.50354,-47.53891
 9: 156 -4.954156,-4.954156,-3.675864,-3.675864 -47.50354,-45.39022,-45.39022,-47.50354
10: 157 -3.675864,-3.675864,-3.706356,-3.706356 -45.39022,-45.30724,-45.30724,-45.39022
11: 158 -3.706356,-3.706356,-3.705801,-3.705801 -45.30724,-45.30722,-45.30722,-45.30724
> points
    longitude  latitude  price
 1: -47.50308 -4.953936 3.0616
 2: -47.50308 -4.953936 3.2070
 3: -47.50308 -4.953936 3.0630
 4: -47.50308 -4.953936 3.0603
 5: -47.50308 -4.953936 3.0460
 6: -47.50308 -4.953936 2.9900
 7: -49.07035 -5.283658 3.3130
 8: -49.08054 -5.347284 3.3900
 9: -49.08054 -5.347284 3.3620
10: -49.21726 -5.338270 3.3900
11: -49.08050 -5.347255 3.4000
12: -49.08042 -5.347248 3.3220
13: -49.08190 -5.359508 3.3130
14: -49.08046 -5.347277 3.3560

Я хочу сгенерировать среднюю цену среди всех точек, которые помещаются внутри каждого многоугольника, для каждого многоугольника.

Прямо сейчас я использую sp::point.in.polygon, чтобы получить индексы для всех точек, которые вписываются в данный полигон, а затем получить его среднюю цену.

w <- lapply(1:nrow(polygons),
            function(tt) {
              ind <- point.in.polygon(points$latitude, points$longitude,
                                      polygons$pol_lat[[tt]], polygons$pol_lng[[tt]]) > 0
              med <- mean(points$price[ind])
              return(med)
            }
)
> unlist(w)
 [1]      NaN 3.361857 3.313000      NaN      NaN      NaN      NaN      NaN 3.071317      NaN      NaN

Однако это, очевидно, медленно. Любые идеи о том, как сделать это быстрее, возможно, используя data.table или dplyr (или любые другие средства)?

Данные следуют

> dput(polygons)
structure(list(id = 148:158, pol_lat = list(c(-4.2502356, -4.2502356, 
-4.2546403, -4.2546403), c(-4.2546403, -4.2546403, -5.3616014, 
-5.3616014), c(-5.3616014, -5.3616014, -5.2122078, -5.2122078
), c(-5.2122078, -5.2122078, -5.0028781, -5.0028781), c(-5.0028781, 
-5.0028781, -5.0800181, -5.0800181), c(-5.0800181, -5.0800181, 
-5.0798186, -5.0798186), c(-5.0798186, -5.0798186, -5.1556063, 
-5.1556063), c(-5.1556063, -5.1556063, -4.9541564, -4.9541564
), c(-4.9541564, -4.9541564, -3.6758637, -3.6758637), c(-3.6758637, 
-3.6758637, -3.706356, -3.706356), c(-3.706356, -3.706356, -3.7058011, 
-3.7058011)), pol_lng = list(c(-49.9462826, -49.9449427, -49.9449427, 
-49.9462826), c(-49.9449427, -49.0790599, -49.0790599, -49.9449427
), c(-49.0790599, -49.0446868, -49.0446868, -49.0790599), c(-49.0446868, 
-48.4866355, -48.4866355, -49.0446868), c(-48.4866355, -48.436988, 
-48.436988, -48.4866355), c(-48.436988, -48.4247989, -48.4247989, 
-48.436988), c(-48.4247989, -47.5389072, -47.5389072, -48.4247989
), c(-47.5389072, -47.5035404, -47.5035404, -47.5389072), c(-47.5035404, 
-45.3902168, -45.3902168, -47.5035404), c(-45.3902168, -45.3072392, 
-45.3072392, -45.3902168), c(-45.3072392, -45.3072216, -45.3072216, 
-45.3072392))), row.names = c(NA, -11L), class = c("data.table", 
"data.frame"), .internal.selfref = <pointer: 0x00000000025e1ef0>)
> dput(points)
structure(list(longitude = c(-47.5030772, -47.5030772, -47.5030772, 
-47.5030772, -47.5030772, -47.5030772, -49.0703469, -49.0805422, 
-49.0805422, -49.217259, -49.0804978, -49.0804181, -49.0818997, 
-49.0804625), latitude = c(-4.9539357, -4.9539357, -4.9539357, 
-4.9539357, -4.9539357, -4.9539357, -5.283658, -5.3472839, -5.3472839, 
-5.3382696, -5.3472551, -5.347248, -5.3595084, -5.3472768), price = c(3.0616, 
3.207, 3.063, 3.0603, 3.046, 2.99, 3.313, 3.39, 3.362, 3.39, 
3.4, 3.322, 3.313, 3.356)), row.names = c(NA, -14L), class = c("data.table", 
"data.frame"), .internal.selfref = <pointer: 0x00000000025e1ef0>)

r sp point-in-polygon
person Felipe Alvarenga    schedule 31.07.2018    source источник
comment
Может помочь: data-set/" title="как проверить, находится ли точка в многоугольнике, эффективно используя r для большого набора данных"> stackoverflow.com/questions/36683825/ и stackoverflow.com/questions/49828692 /   -  person Stéphane Laurent    schedule 31.07.2018
comment
Я прошел по второй ссылке, но ptinpoly особого прироста скорости не дает. Посмотрим, помогут ли ответы в первой ссылке. Спасибо   -  person Felipe Alvarenga    schedule 31.07.2018
comment
не являюсь экспертом по пространственной статистике, но ваш lapply выполняет итерацию по всем полигонам, что может быть неэффективным, если количество точек намного меньше, чем количество полигонов (как, по-видимому, в вашем случае). Вы пытались вместо этого перебирать точки, а затем смотреть, принадлежит ли она многоугольнику? Хорошая вещь в этом подходе заключается в том, что если полигоны представляют собой раздел, то, как только вы нашли многоугольник, которому принадлежит точка, вы можете перестать искать в другом месте и перейти к следующей точке.   -  person mbiron    schedule 10.08.2018
comment
Прямо сейчас чувствую себя глупо, потому что не подумал об этом. Спасибо @mbiron   -  person Felipe Alvarenga    schedule 10.08.2018
comment
Нет проблем, чувак. На самом деле, вы, возможно, могли бы сделать лучше. Для каждой точки, вместо случайной проверки каждого полигона (что уже является преимуществом, поскольку в среднем вы будете проверять только половину полигонов), вы можете сначала вычислить расстояния от точки до центроидов всех полигонов (это может быть сделать быстро), а затем оценить полигоны в порядке возрастания расстояния. Таким образом, вы должны оценить всего несколько полигонов, пока не найдете тот, у которого есть точка внутри.   -  person mbiron    schedule 10.08.2018
comment
Кроме того, вместо сортировки всего вектора расстояний длиной 13 миллионов можно использовать аргумент partial функции sort для получения ближайших 10 или 100 полигонов. Было бы действительно странно, если бы точка не принадлежала ни к одному из ближайших 100 элементов, а это, вероятно, невозможно при обычной настройке пространственной статистики.   -  person mbiron    schedule 10.08.2018
comment
Мои полигоны могут иметь перекрывающиеся области. Я уже проверяю только ближайшие полигоны к каждой точке. Инвертирование цикла, т. е. перебор точек вместо этого, помогло.   -  person Felipe Alvarenga    schedule 10.08.2018
comment
Ммм... не знаю, слежу ли я, здесь. Кажется, что если вы не можете сделать какое-то предположение относительно размера полигонов, проверка только ближайших полигонов может привести к ошибкам. Теоретически ничто не препятствует тому, чтобы точка, находящаяся далеко от центроида определенного многоугольника, все еще находилась внутри него, если многоугольник большой.   -  person lbusett    schedule 10.08.2018
comment
они совсем маленькие. Мои многоугольники представляют собой прямоугольники над направлениями Google Maps, поворачивающие точки над всем континентом Южной Америки, поэтому довольно сложно пропустить точку.   -  person Felipe Alvarenga    schedule 10.08.2018

Ответы (1)


arrow_upward
1
arrow_downward

Если ваши "многоугольники" всегда являются прямоугольниками, как в примере, можно использовать пространственный индекс QuadTree, реализованный в пакете SearchTrees, чтобы повысить скорость определения того, какие точки попадают в каждый многоугольник.

Это может дать вам довольно большой прирост скорости благодаря меньшему количеству "сравнений", которые должны быть выполнены пространственным индексированием, тем больше, чем больше количество точек в вашем наборе данных.

Например:

library(SearchTrees)
library(magrittr)

# Create a "beefier" test dataset based on your data: 14000 pts 
# over 45000 polygons

for (i in 1:10) points   <- rbind(points, points + runif(length(points)))
for (i in 1:12) polygons <- rbind(polygons, polygons)


# Compute limits of the polygons
min_xs <- lapply(polygons$pol_lng , min) %>% unlist()
max_xs <- lapply(polygons$pol_lng , max) %>% unlist()
min_ys <- lapply(polygons$pol_lat , min) %>% unlist()
max_ys <- lapply(polygons$pol_lat, max) %>% unlist()
xlims <- cbind(min_xs, max_xs)
ylims <- cbind(min_ys, max_ys)

# Create the quadtree
tree = SearchTrees::createTree(cbind(points[1],points[2]))

#☻ extract averages, looping over polygons ----
t1 <- Sys.time()
w <- lapply(1:nrow(polygons), 
            function(tt) {
              ind <- SearchTrees::rectLookup(
                tree, 
                xlims = xlims[tt,],
                ylims = ylims[tt,]))
              mean(points$price[ind])

              })
Sys.time() - t1

Разница во времени 2,945789 сек.

w1 <- unlist(w)

На моем старом ноутбуке эта «наивная» реализация работает более чем в 10 раз быстрее, чем ваш первоначальный подход к тестовым данным:

t1 <- Sys.time()
w <- lapply(1:nrow(polygons),
            function(tt) {
              ind <- sp::point.in.polygon(points$latitude, points$longitude,
                                      polygons$pol_lat[[tt]], polygons$pol_lng[[tt]]) > 0
              med <- mean(points$price[ind])
              return(med)
            }
)
Sys.time() - t1
w2 <- unlist(w)

Разница во времени 40,36493 сек.

, с теми же результатами:

> all.equal(w1, w2)
[1] TRUE

Общее улучшение скорости будет зависеть от того, как ваши точки «сгруппированы» по пространственному экстенту и по отношению к полигонам.

Учтите, что вы можете использовать этот подход также, если ваши многоугольники не прямоугольные, сначала извлекая точки, включенные в bbox каждого многоугольника, а затем находя точки «внутри» многоугольника с помощью более стандартного подхода.

Также примите во внимание, что задача удручающе параллельна, так что вы можете легко повысить производительность, используя подход foreach или parlapply к полигонам.

ХТХ!

person lbusett    schedule 09.08.2018
comment
это очень помогло. Странно то, что это было в 4 раза быстрее, а не в 10 раз. Тем не менее, значительное улучшение - person Felipe Alvarenga; 10.08.2018