como sobrepor um polígono sobre SpatialPointsDataFrame e preservar os dados SPDF?

Eu tenho um SpatialPointsDataFramecom alguns dados adicionais. Gostaria de extrair esses pontos dentro de um polígono e, ao mesmo tempo, preservar o SPDFobjeto e seus dados correspondentes.

Até agora, tive pouca sorte e recorri à correspondência e à fusão por meio de um ID comum, mas isso só funciona porque eu agrupei os dados com o IDS individual.

Aqui está um exemplo rápido: estou procurando pontos dentro do quadrado vermelho.

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

A abordagem mais óbvia seria usar over, mas isso retorna os dados do polígono.

> over(pts, ply)
    polyvar
1        NA
2       357
3       357
4        NA
5       357
6       357

Da sp::overajuda:

 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.

Portanto, se você converter seu SpatialPolygonsDataFramepara SpatialPolygonsvocê, receberá de volta um vetor de índices e poderá definir seus pontos em 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
> 

Para os que duvidam, aqui estão as evidências de que a sobrecarga da conversão não é um problema:

Duas funções - primeiro o método de Jeffrey Evans, depois meu original, depois minha conversão hackeada, depois uma versão baseada na gIntersectsresposta de Josh O'Brien:

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,]
}

Agora, para um exemplo do mundo real, eu espalhei alguns pontos aleatórios pelo columbusconjunto de dados:

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

Parece bom

plot(columbus)
points(pts)

Verifique se as funções estão fazendo a mesma coisa:

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

E execute 500 vezes para fazer comparações:

> 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 

De acordo com minha intuição, não é uma sobrecarga muito grande; na verdade, pode ser uma sobrecarga menor do que converter todos os índices de linha em caracteres e vice-versa ou executar na.omit para obter valores ausentes. Que, aliás, leva a outro modo de falha da evansfunção ...

Se uma linha do quadro de dados de polígonos for toda NA(o que é perfeitamente válido), a sobreposição SpatialPolygonsDataFramede pontos nesse polígono produzirá um quadro de dados de saída com todos os NAs, que evans()serão eliminados:

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

MAS gIntersectsé mais rápido, mesmo com a varredura da matriz para verificar as interseções no código R e não no código C. Suspeito que sejam as prepared geometryhabilidades do GEOS, criando índices espaciais - sim, com prepared=FALSEisso demora um pouco mais, cerca de 5,5 segundos.

Estou surpreso que não haja uma função para retornar diretamente os índices ou os pontos. Quando escrevi splancshá 20 anos, as funções point-in-polygon tinham ambas ...

sp fornece um formato mais curto para selecionar recursos com base na interseção espacial, seguindo o exemplo do OP:

pts[ply,]

a partir de:

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

Nos bastidores, é a abreviação de

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

O que deve ser observado é que existe um geometrymétodo que descarta atributos: overaltera o comportamento se o segundo argumento tiver atributos ou não (essa foi a confusão do OP). Isso funciona em todas as classes Spatial * sp, embora alguns overmétodos exijam rgeos, consulte esta vinheta para obter detalhes, por exemplo, o caso de várias correspondências para polígonos sobrepostos.

Você estava no caminho certo com mais de. Os nomes de nomes de objeto do objeto retornado correspondem ao índice de linhas dos pontos. Você pode implementar sua abordagem exata com apenas algumas linhas de código adicionais.

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) 

É isso que você procura?

Uma observação, na edição: a chamada de empacotamento para apply()é necessária para fazer isso funcionar com SpatialPolygonsobjetos arbitrários , possivelmente contendo mais de um recurso de polígono. Agradeço ao @Spacedman por me incentivar a demonstrar como aplicar isso ao caso mais geral.

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)

Aqui está uma abordagem possível usando o rgeospacote. Basicamente, utiliza a gIntersectionfunção que permite interceptar dois spobjetos. Ao extrair os IDs desses pontos que estão dentro do polígono, você poderá subconfigurar o original SpatialPointsDataFrame, mantendo todos os dados correspondentes. O código é quase auto-explicativo, mas se houver alguma dúvida, não hesite em perguntar!

# 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

Existe uma solução extremamente simples usando a spatialEcobiblioteca.

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)

Veja o resultado:

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