jak nałożyć wielobok na SpatialPointsDataFrame i zachować dane SPDF?

Mam SpatialPointsDataFrametrochę dodatkowych danych. Chciałbym wydobyć te punkty wewnątrz wielokąta i jednocześnie zachowaćSPDF obiekt i odpowiadające mu dane.

Do tej pory nie miałem szczęścia i korzystałem z dopasowywania i łączenia za pomocą wspólnego identyfikatora, ale działa to tylko dlatego, że mam dane w siatce z indywidualnymi identyfikatorami.

Oto szybki przykład: szukam punktów wewnątrz czerwonego kwadratu.

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

Najbardziej oczywistym podejściem byłoby użycie over, ale zwraca dane z wielokąta.

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

Z sp::overpomocy:

 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.

Więc jeśli przekonwertujesz SpatialPolygonsDataFramena SpatialPolygons, otrzymasz wektor indeksów i możesz podzielić swoje punkty na 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
> 

Dla wątpiących, oto dowód, że narzut związany z konwersją nie stanowi problemu:

Dwie funkcje - najpierw metoda Jeffreya Evansa, potem moja oryginalna, potem moja zhakowana konwersja, a następnie wersja oparta na gIntersectsodpowiedzi Josha O'Briena:

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

Teraz na przykład w świecie rzeczywistym rozrzuciłem kilka losowych punktów nad columbuszestawem danych:

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

Wygląda dobrze

plot(columbus)
points(pts)

Sprawdź, czy funkcje robią to samo:

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

I uruchom 500 razy na potrzeby testów porównawczych:

> 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 

Według mojej intuicji, nie jest to wielki narzut, w rzeczywistości może być mniejszy narzut niż konwertowanie wszystkich indeksów wierszy na znak iz powrotem lub uruchamianie na.omit, aby uzyskać brakujące wartości. Co nawiasem mówiąc, prowadzi do innego trybu awarii evansfunkcji ...

Jeśli wiersz ramki danych wielokątów jest wszystkim NA(co jest całkowicie poprawne), wówczas nakładka z SpatialPolygonsDataFramepunktami w tym wielokącie wytworzy wyjściową ramkę danych ze wszystkimi NAs, które evans()następnie spadną:

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

ALE gIntersectsjest szybszy, nawet z koniecznością zamiatania macierzy, aby sprawdzić przecięcia w R zamiast w kodzie C. Podejrzewam, że ma prepared geometryumiejętności GEOS, tworząc indeksy przestrzenne - tak, przy prepared=FALSEczym zajmuje to nieco więcej, około 5,5 sekundy.

Dziwię się, że nie ma funkcji ani prostego zwracania wskaźników, ani punktów. Kiedy pisałem splancs20 lat temu, funkcje punkt-w-wielokącie miały zarówno ...

sp zapewnia krótszy formularz do wybierania obiektów w oparciu o przecięcie przestrzenne, zgodnie z przykładem PO:

pts[ply,]

od:

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

Za kulisami jest to skrót

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

Należy zauważyć, że istnieje geometrymetoda, która upuszcza atrybuty: overzmienia zachowanie, jeśli drugi argument ma atrybuty, czy nie (to było zamieszanie OP). Działa to we wszystkich klasach Spatial * sp, chociaż niektóre overmetody tego wymagają rgeos, zobacz tę winietę, aby uzyskać szczegółowe informacje, np. Przypadek wielokrotnego dopasowania nakładających się wielokątów.

Z tobą byłeś na dobrej drodze. Nazwy zwracanego obiektu odpowiadają indeksowi rzędów punktów. Możesz wdrożyć swoje dokładne podejście za pomocą zaledwie kilku dodatkowych wierszy kodu.

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) 

Czy o to ci chodzi?

Jedna uwaga na temat edycji: apply()konieczne jest wywołanie zawijania do dowolnych SpatialPolygonsobiektów, które mogą zawierać więcej niż jeden element wielokąta. Dzięki @Spacedman za nakłonienie mnie do zademonstrowania, jak zastosować to w bardziej ogólnym przypadku.

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)

Oto możliwe podejście przy użyciu rgeospakietu. Zasadniczo korzysta z gIntersectionfunkcji, która umożliwia przecięcie dwóch spobiektów. Wyodrębniając identyfikatory punktów znajdujących się w obrębie wielokąta, możesz następnie podzielić swój oryginał SpatialPointsDataFrame, zachowując wszystkie odpowiednie dane. Kod jest prawie samoobjaśniający, ale jeśli są jakieś pytania, możesz je zadać!

# 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

Istnieje niezwykle proste rozwiązanie z wykorzystaniem spatialEcobiblioteki.

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)

Sprawdź wynik:

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