Как наложить многоугольник на SpatialPointsDataFrame и сохранить данные SPDF?

17

у меня есть 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
Роман Луштрик
источник
1
Спасибо за предоставленный воспроизводимый пример. Всегда помогает при попытке понять проблему!
fdetsch

Ответы:

21

Из sp::overсправки:

 x = "SpatialPoints", y = "SpatialPolygons" returns a numeric
      vector of length equal to the number of points; the number is
      the index (number) of the polygon of ‘y’ in which a point
      falls; NA denotes the point does not fall in a polygon; if a
      point falls in multiple polygons, the last polygon is
      recorded.

Так что если вы конвертируете SpatialPolygonsDataFrameв SpatialPolygonsвас, вы получите вектор индексов, и вы можете установить свои баллы на NA:

> over(pts,as(ply,"SpatialPolygons"))
  [1] NA  1  1 NA  1  1 NA NA  1  1  1 NA NA  1  1  1  1  1 NA NA NA  1 NA  1 NA
 [26]  1  1  1 NA NA NA NA NA  1  1 NA NA NA  1  1  1 NA  1  1  1 NA NA NA  1  1
 [51]  1 NA NA NA  1 NA  1 NA  1 NA NA  1 NA  1  1 NA  1  1 NA  1 NA  1  1  1  1
 [76]  1  1  1  1  1 NA NA NA  1 NA  1 NA NA NA NA  1  1 NA  1 NA NA  1  1  1 NA

> nrow(pts)
[1] 100
> pts = pts[!is.na(over(pts,as(ply,"SpatialPolygons"))),]
> nrow(pts)
[1] 54
> head(pts@data)
         var1 var2
2  0.04001092    v
3  0.58108350    v
5  0.85682609    q
6  0.13683264    y
9  0.13968804    m
10 0.97144627    o
> 

Для сомневающихся, вот доказательство того, что затраты на преобразование не являются проблемой:

Две функции - сначала метод Джеффри Эванса, затем мой оригинал, затем мое взломанное преобразование, а затем версия, основанная на gIntersectsответе Джоша О'Брайена:

evans <- function(pts,ply){
  prid <- over(pts,ply)
  ptid <- na.omit(prid) 
  pt.poly <- pts[as.numeric(as.character(row.names(ptid))),]
  return(pt.poly)
}

rowlings <- function(pts,ply){
  return(pts[!is.na(over(pts,as(ply,"SpatialPolygons"))),])
}

rowlings2 <- function(pts,ply){
  class(ply) <- "SpatialPolygons"
  return(pts[!is.na(over(pts,ply)),])
}

obrien <- function(pts,ply){
pts[apply(gIntersects(columbus,pts,byid=TRUE),1,sum)==1,]
}

Теперь для реального примера я разбросал несколько случайных точек по columbusнабору данных:

require(spdep)
example(columbus)
pts=data.frame(
    x=runif(100,5,12),
    y=runif(100,10,15),
    z=sample(letters,100,TRUE))
coordinates(pts)=~x+y

Выглядит хорошо

plot(columbus)
points(pts)

Проверьте, что функции делают то же самое:

> identical(evans(pts,columbus),rowlings(pts,columbus))
[1] TRUE

И запустить 500 раз для бенчмаркинга:

> system.time({for(i in 1:500){evans(pts,columbus)}})
   user  system elapsed 
  7.661   0.600   8.474 
> system.time({for(i in 1:500){rowlings(pts,columbus)}})
   user  system elapsed 
  6.528   0.284   6.933 
> system.time({for(i in 1:500){rowlings2(pts,columbus)}})
   user  system elapsed 
  5.952   0.600   7.222 
> system.time({for(i in 1:500){obrien(pts,columbus)}})
  user  system elapsed 
  4.752   0.004   4.781 

Согласно моей интуиции, это не большие издержки, на самом деле это может быть меньше, чем преобразование всех индексов строк в символьные и обратно, или запуск na.omit для получения пропущенных значений. Что случайно приводит к другому отказу режима работы evans...

Если строка фрейма данных полигонов - это все NA(что совершенно правильно), то наложение с SpatialPolygonsDataFrameточками для этого полигона создаст фрейм выходных данных со всеми NAs, который evans()затем будет отброшен:

> columbus@data[1,]=rep(NA,20)
> columbus@data[5,]=rep(NA,20)
> columbus@data[17,]=rep(NA,20)
> columbus@data[15,]=rep(NA,20)
> set.seed(123)
> pts=data.frame(x=runif(100,5,12),y=runif(100,10,15),z=sample(letters,100,TRUE))
> coordinates(pts)=~x+y
> identical(evans(pts,columbus),rowlings(pts,columbus))
[1] FALSE
> dim(evans(pts,columbus))
[1] 27  1
> dim(rowlings(pts,columbus))
[1] 28  1
> 

НО gIntersectsбыстрее, даже при необходимости развернуть матрицу для проверки пересечений в R, а не в C-коде. Я подозреваю, что это prepared geometryнавыки GEOS, создание пространственных индексов - да, сprepared=FALSE это занимает немного больше времени, около 5,5 секунд.

Я удивлен, что нет функции, которая бы напрямую возвращала индексы или очки. Когда я писал splancs20 лет назад, функции точка-многоугольник имели оба ...

Spacedman
источник
Отлично, это также работает для нескольких полигонов (в ответ Джошуа я добавил пример для игры).
Роман Луштрик
При использовании больших наборов данных полигонов приведение к объекту SpatialPolygons приводит к большим накладным расходам и не является необходимым. Применение «over» к SpatialPolygonsDataFrame возвращает индекс строки, который можно использовать для поднабора точек. Смотрите мой пример ниже.
Джеффри Эванс
Много накладных расходов? По сути, он просто берет слот @polygons из объекта SpatialPolygonsDataFrame. Вы даже можете «подделать» его, переназначив класс SpatialPolygonsDataFrame как «SpatialPolygons» (хотя это и не рекомендуется). Все, что будет использовать геометрию, в любом случае должно будет получить этот слот на некоторой стадии, так что, условно говоря, это вообще не накладные расходы. В любом случае, это незначительно в любом реальном приложении, где вы собираетесь выполнять множество тестов с точечными полигонами.
Spacedman
При учете накладных расходов важны не только соображения скорости. При создании нового объекта в пространстве имен R вы используете необходимую оперативную память. Если это не проблема для небольших наборов данных, это повлияет на производительность при работе с большими данными. R действительно демонстрирует линейное снижение производительности. По мере того как данные становятся больше, производительность возрастает. Если вам не нужно создавать дополнительный объект, зачем вам?
Джеффри Эванс
1
Мы не знали этого, пока я не проверил это прямо сейчас.
Spacedman
13

sp обеспечивает более короткую форму для выбора объектов на основе пространственного пересечения, следуя примеру OP:

pts[ply,]

по состоянию на:

points(pts[ply,], col = 'red')

За кулисами это сокращение от

pts[!is.na(over(pts, geometry(ply))),]

Следует отметить, что существует geometryметод, который отбрасывает атрибуты: overизменяет поведение, если у его второго аргумента есть атрибуты или нет (это была путаница OP). Это работает для всех классов Spatial * в sp, хотя некоторые overметоды требуют rgeos, см. Эту виньетку для деталей, например, в случае множественных совпадений для перекрывающихся полигонов.

Эдзер Пебесма
источник
Хорошо знать! Я не знал о методе геометрии.
Джеффри Эванс
2
Добро пожаловать на наш сайт, Эдзер - приятно видеть вас здесь!
whuber
1
Спасибо, Билл, на stat.ethz.ch/pipermail/r-sig-geo становится все тише , или, может быть, мы должны разработать программное обеспечение, которое вызывает больше проблем! ;-)
Эдзер Пебесма
6

Вы были на правильном пути с более чем. Имена строк возвращаемого объекта соответствуют индексу строк точек. Вы можете реализовать свой точный подход с помощью всего лишь нескольких строк кода.

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

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))

# Subset points intersecting polygon
prid <- over(pts,ply)
  ptid <- na.omit(prid) 
    pt.poly <- pts[as.numeric(as.character(row.names(ptid))),]  

plot(pts)
  axis(1); axis(2)
    plot(ply, add=TRUE, border="red")
      plot(pt.poly,pch=19,add=TRUE) 
Джеффри Эванс
источник
Неправильно - имена строк возвращаемого объекта соответствуют индексу строк in_this_case - в общем, имена строк кажутся именами строк точек - которые могут даже не быть числовыми. Вы можете изменить свое решение, чтобы оно соответствовало характеру, что могло бы сделать его немного более устойчивым.
Spacedman
@ Sapcedman, не будь таким догматичным. Решение не является неправильным. Если вы хотите поместить точки в набор полигонов или назначить значения полигонов для точек, то функция over работает без принуждения. Есть несколько, чтобы выполнить пешеходный переход, как только вы получите результирующий объект. Решение по принуждению к объекту SpatialPolygon создает значительные необходимые накладные расходы, поскольку эту операцию можно выполнить непосредственно над объектом SpatialPolygonDataFrame. Кстати, прежде чем редактировать пост, убедитесь, что вы правы. Термин библиотека и пакет используются взаимозаменяемо в R.
Джеффри Эванс
Я добавил в свой пост несколько тестов и обнаружил еще одну проблему с вашей функцией. Также «Пакеты - это коллекции функций R, данных и скомпилированного кода в четко определенном формате. Каталог, в котором хранятся пакеты, называется библиотекой»
Spacedman
Хотя вы технически верны в отношении «пакета» и «библиотеки», вы спорите о семантике. У меня только что был запрос редактора Ecological Modeling, чтобы мы изменили наше использование «package» (что на самом деле является моим предпочтением) на «library». Я хочу сказать, что они становятся взаимозаменяемыми терминами и предметом предпочтения.
Джеффри Эванс
1
«технически правильно», как однажды заметил д-р Шелдон Купер, «лучший вариант». Этот редактор технически не прав, что является худшим видом ошибок.
Spacedman
4

Это то, что вы после?

Одно примечание по редактированию: вызов обтекания apply()необходим для работы с произвольными SpatialPolygonsобъектами, возможно, содержащими более одного полигонального объекта. Спасибо @Spacedman за то, что он подтолкнул меня продемонстрировать, как применить это в более общем случае.

library(rgeos)
pp <- pts[apply(gIntersects(pts, ply, byid=TRUE), 2, any),]


## Confirm that it works
pp[1:5,]
#              coordinates       var1 var2
# 2 (-0.583205, -0.877737) 0.04001092    v
# 3   (0.394747, 0.702048) 0.58108350    v
# 5    (0.7668, -0.946504) 0.85682609    q
# 6    (0.31746, 0.641628) 0.13683264    y
# 9   (-0.469015, 0.44135) 0.13968804    m

plot(pts)
plot(ply, border="red", add=TRUE)
plot(pp, col="red", add=TRUE)
Джош О'Брайен
источник
Сбоит ужасно, если plyимеет более одной функции, потому что gIntersectsвозвращает матрицу с одной строкой для каждой функции. Вы, вероятно, можете подмести строки для значения ИСТИНА.
Космонавт
@ Spacepman - Бинго. Нужно сделать apply(gIntersects(pts, ply, byid=TRUE), 2, any). На самом деле, я пойду дальше и переключу ответ на этот вопрос, поскольку он охватывает также случай одного многоугольника.
Джош О'Брайен
Ах, any. Это может быть немного быстрее, чем версия, которую я только что тестировал.
космонавт
@Spacedman - Из моих быстрых тестов, это выглядит , как obrienи rowlings2запустить шеи и декольте, с , obrien возможно , 2% быстрее.
Джош О'Брайен
@ JoshO'Brien, как можно использовать этот ответ на многих полигонах? То есть ppдолжен иметь, IDкоторый указывает, в каком полигоне расположены точки.
code123
4

Вот возможный подход с использованием rgeosпакета. По сути, он использует gIntersectionфункцию, которая позволяет пересекать два spобъекта. Извлекая идентификаторы тех точек, которые лежат в пределах многоугольника, вы впоследствии сможете поднастроить свой исходныйSpatialPointsDataFrame , сохранив все соответствующие данные. Код почти самоочевиден, но если есть какие-либо вопросы, пожалуйста, не стесняйтесь спрашивать!

# Required package
library(rgeos)

# Intersect polygons and points, keeping point IDs
pts.intersect <- gIntersection(ply, pts, byid = TRUE)

# Extract point IDs from intersected data
pts.intersect.strsplit <- strsplit(dimnames(pts.intersect@coords)[[1]], " ")
pts.intersect.id <- as.numeric(sapply(pts.intersect.strsplit, "[[", 2))

# Subset original SpatialPointsDataFrame by extracted point IDs
pts.extract <- pts[pts.intersect.id, ]

head(coordinates(pts.extract))
              x          y
[1,] -0.5832050 -0.8777367
[2,]  0.3947471  0.7020481
[3,]  0.7667997 -0.9465043
[4,]  0.3174604  0.6416281
[5,] -0.4690151  0.4413502
[6,]  0.4765213  0.6068021

head(pts.extract)
         var1 var2
2  0.04001092    v
3  0.58108350    v
5  0.85682609    q
6  0.13683264    y
9  0.13968804    m
10 0.97144627    o
fdetsch
источник
1
Должно tmpбыть pts.intersect? Кроме того, подобный синтаксический анализ возвращаемых dimnames зависит от недокументированного поведения.
Космонавт
@Spacedman, ты прав tmp, забыл удалить его при завершении кода. Кроме того, вы правы в разборе dimnames. Это было своего рода быстрое решение, позволяющее быстро ответить на вопрос, и, конечно, есть лучшие (и более универсальные) подходы, например, ваш :-)
fdetsch
1

Существует чрезвычайно простое решение с использованием spatialEcoбиблиотеки.

library(spatialEco)

# intersect points in polygon
  pts <- point.in.poly(pts, ply)

# check plot
  plot(ply)
  plot(a, add=T)

# convert to data frame, keeping your data
  pts<- as.data.frame(pts)

Проверьте результат:

pts

>             x          y       var1 var2 polyvar
> 2  -0.5832050 -0.8777367 0.04001092    v     357
> 3   0.3947471  0.7020481 0.58108350    v     357
> 5   0.7667997 -0.9465043 0.85682609    q     357
> 6   0.3174604  0.6416281 0.13683264    y     357
> 9  -0.4690151  0.4413502 0.13968804    m     357
> 10  0.4765213  0.6068021 0.97144627    o     357
rafa.pereira
источник