у меня есть SpatialPointsDataFrame
некоторые дополнительные данные. Я хотел бы извлечь эти точки внутри многоугольника и в то же время сохранить SPDF
объект и соответствующие ему данные.
До сих пор мне не везло, и я прибегал к сопоставлению и объединению с помощью общего идентификатора, но это работает только потому, что у меня есть привязанные к сетке данные с отдельными идентификаторами.
Вот быстрый пример, я ищу точки внутри красного квадрата.
library(sp)
set.seed(357)
pts <- data.frame(x = rnorm(100), y = rnorm(100), var1 = runif(100), var2 = sample(letters, 100, replace = TRUE))
coordinates(pts) <- ~ x + y
class(pts)
plot(pts)
axis(1); axis(2)
ply <- matrix(c(-1,-1, 1,-1, 1,1, -1,1, -1,-1), ncol = 2, byrow = TRUE)
ply <- SpatialPolygons(list(Polygons(list(Polygon(ply)), ID = 1)))
ply <- SpatialPolygonsDataFrame(Sr = ply, data = data.frame(polyvar = 357))
plot(ply, add = TRUE, border = "red")
Наиболее очевидный подход заключается в использовании over
, но это возвращает данные из многоугольника.
> over(pts, ply)
polyvar
1 NA
2 357
3 357
4 NA
5 357
6 357
Ответы:
Из
sp::over
справки:Так что если вы конвертируете
SpatialPolygonsDataFrame
вSpatialPolygons
вас, вы получите вектор индексов, и вы можете установить свои баллы наNA
:Для сомневающихся, вот доказательство того, что затраты на преобразование не являются проблемой:
Две функции - сначала метод Джеффри Эванса, затем мой оригинал, затем мое взломанное преобразование, а затем версия, основанная на
gIntersects
ответе Джоша О'Брайена:Теперь для реального примера я разбросал несколько случайных точек по
columbus
набору данных:Выглядит хорошо
Проверьте, что функции делают то же самое:
И запустить 500 раз для бенчмаркинга:
Согласно моей интуиции, это не большие издержки, на самом деле это может быть меньше, чем преобразование всех индексов строк в символьные и обратно, или запуск na.omit для получения пропущенных значений. Что случайно приводит к другому отказу режима работы
evans
...Если строка фрейма данных полигонов - это все
NA
(что совершенно правильно), то наложение сSpatialPolygonsDataFrame
точками для этого полигона создаст фрейм выходных данных со всемиNA
s, которыйevans()
затем будет отброшен:НО
gIntersects
быстрее, даже при необходимости развернуть матрицу для проверки пересечений в R, а не в C-коде. Я подозреваю, что этоprepared geometry
навыки GEOS, создание пространственных индексов - да, сprepared=FALSE
это занимает немного больше времени, около 5,5 секунд.Я удивлен, что нет функции, которая бы напрямую возвращала индексы или очки. Когда я писал
splancs
20 лет назад, функции точка-многоугольник имели оба ...источник
sp
обеспечивает более короткую форму для выбора объектов на основе пространственного пересечения, следуя примеру OP:по состоянию на:
За кулисами это сокращение от
Следует отметить, что существует
geometry
метод, который отбрасывает атрибуты:over
изменяет поведение, если у его второго аргумента есть атрибуты или нет (это была путаница OP). Это работает для всех классов Spatial * вsp
, хотя некоторыеover
методы требуютrgeos
, см. Эту виньетку для деталей, например, в случае множественных совпадений для перекрывающихся полигонов.источник
Вы были на правильном пути с более чем. Имена строк возвращаемого объекта соответствуют индексу строк точек. Вы можете реализовать свой точный подход с помощью всего лишь нескольких строк кода.
источник
Это то, что вы после?
Одно примечание по редактированию: вызов обтекания
apply()
необходим для работы с произвольнымиSpatialPolygons
объектами, возможно, содержащими более одного полигонального объекта. Спасибо @Spacedman за то, что он подтолкнул меня продемонстрировать, как применить это в более общем случае.источник
ply
имеет более одной функции, потому чтоgIntersects
возвращает матрицу с одной строкой для каждой функции. Вы, вероятно, можете подмести строки для значения ИСТИНА.apply(gIntersects(pts, ply, byid=TRUE), 2, any)
. На самом деле, я пойду дальше и переключу ответ на этот вопрос, поскольку он охватывает также случай одного многоугольника.any
. Это может быть немного быстрее, чем версия, которую я только что тестировал.obrien
иrowlings2
запустить шеи и декольте, с ,obrien
возможно , 2% быстрее.pp
должен иметь,ID
который указывает, в каком полигоне расположены точки.Вот возможный подход с использованием
rgeos
пакета. По сути, он используетgIntersection
функцию, которая позволяет пересекать дваsp
объекта. Извлекая идентификаторы тех точек, которые лежат в пределах многоугольника, вы впоследствии сможете поднастроить свой исходныйSpatialPointsDataFrame
, сохранив все соответствующие данные. Код почти самоочевиден, но если есть какие-либо вопросы, пожалуйста, не стесняйтесь спрашивать!источник
tmp
бытьpts.intersect
? Кроме того, подобный синтаксический анализ возвращаемых dimnames зависит от недокументированного поведения.tmp
, забыл удалить его при завершении кода. Кроме того, вы правы в разбореdimnames
. Это было своего рода быстрое решение, позволяющее быстро ответить на вопрос, и, конечно, есть лучшие (и более универсальные) подходы, например, ваш :-)Существует чрезвычайно простое решение с использованием
spatialEco
библиотеки.Проверьте результат:
источник