Рассчитать все расстояния между двумя наборами точек, используя st_distance

Question

У меня есть два набора точек, хранящихся в R как объекты sf.Точечный объект x содержит 204 467, а точка y содержит 5 297 точек.

Теоретически я хотел бы рассчитать расстояние от всех точек в x до всех точек в y.Я понимаю, что это создаст зверя из матрицы, но это возможно при использовании st_distance (x, y, by_element = FALSE) в пакете sf примерно за 40 минут на моем рабочем столе i7.

Что я хочу сделатьdo - рассчитать расстояние от всех точек в x до всех точек в y, затем я хочу преобразовать это в data.frame, который содержит все переменные для соответствующей пары x и y точек.Это потому, что я хочу гибкости с точки зрения агрегирования с использованием dplyr, например, я хочу найти количество точек в y, которое находится в пределах 10, 50, 100 км от x и где x $ year

Я успешно создал матрицу расстояний, которая содержит около 1 083 061 699 ячеек.Я знаю, что это очень неэффективный способ сделать это, но он дает гибкость с точки зрения агрегирования.Другие предложения приветствуются.

Ниже приведен код для создания двух точечных объектов и измерения расстояния между ними.Затем я хотел бы преобразовать это в data.frame со всеми переменными из x и y, но здесь я не могу продолжить.

Если предложенный мной рабочий процесс неосуществим, может ли кто-нибудь предложить альтернативное решение для измерения расстояния до всех точек в пределах заранее определенного радиуса и создать data.frame результата со всеми переменными из x и y?

# Create two sf point objects 
set.seed(123)
library(sf)


pts1 <- st_as_sf(x = data.frame(id=seq(1,204467,1),
                                year=sample(seq(from = 1990, to = 2018, by = 1), size = 204467, replace = TRUE),
                                xcoord=sample(seq(from = -180, to = 180, by = 1), size = 204467, replace = TRUE),
                                ycoord=sample(seq(from = -90, to = 90, by = 1), size = 204467, replace = TRUE)),
                 coords=c("xcoord","ycoord"),crs=4326)

pts2 <- st_as_sf(x = data.frame(id=seq(1,5297,1),
                                year=sample(seq(from = 1990, to = 2018, by = 1), size = 5297, replace = TRUE),
                                xcoord=sample(seq(from = -180, to = 180, by = 1), size = 5297, replace = TRUE),
                                ycoord=sample(seq(from = -90, to = 90, by = 1), size = 5297, replace = TRUE)),
                 coords=c("xcoord","ycoord"),crs=4326)

distmat <- st_distance(pts1,pts2,by_element = FALSE)

Answer 1

Я бы рассмотрел подход к этому по-другому.Если у вас есть матрица distmat, вы можете выполнять описанные вами типы расчетов без использования data.frame.Вы можете использовать стандартное подмножество, чтобы найти, какие точки соответствуют указанным критериям.

Например, чтобы найти комбинации точек, в которых pts1$year больше pts2$year, мы можем сделать:

subset_points = outer(pts1$year, pts2$year, `>`)

Затем, чтобы узнать, сколько из них разделено более чем на 100 км, мы можем сделать

library(units)
sum(distmat[subset_points] > (100 * as_units('km', 1)))

Замечание об использовании памяти

Однако вы подходите к этомус объектами sf или data.frame есть вероятность, что вы начнете увеличивать ограничения ОЗУ с 1e9 числами с плавающей запятой в каждой матрице или столбце таблицы data.table.Вместо этого вы можете подумать о преобразовании вашей матрицы расстояний в raster.Затем растр можно сохранить на диске, а не в памяти, и вы можете использовать функции безопасности памяти в пакете raster, чтобы ускорить ваш путь.

Как мы можем использовать растры для работыс диска и сохраните ОЗУ

Мы можем использовать безопасные для памяти растровые операции для таких очень больших матриц, например:

library(raster)

# convert our matrices to rasters, so we can work on them from disk
r = raster(matrix(as.numeric(distmat), length(pts1$id), length(pts2$id)))
s = raster(subset_points)
remove('distmat', 'subset_points')

# now create a raster equal to r, but with zeroes in the cells we wish to exclude from calculation
rs = overlay(r,s,fun=function(x,y){x*y}, filename='out1.tif')     

# find which cells have value greater than x (1e6 in the example)
Big_cells = reclassify(rs, matrix(c(-Inf, 1e6, 0, 1e6, Inf, 1), ncol=3, byrow=TRUE), 'out.tiff', overwrite=T)

# and finally count the cells
N = cellStats(Big_cells, sum)