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

17

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
Roman Luštrik
źródło
1
Dzięki za dostarczenie odtwarzalnego przykładu. Zawsze pomaga, gdy próbujesz zrozumieć problem!
fdetsch

Odpowiedzi:

21

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

Spacedman
źródło
Świetnie, działa to również w przypadku wielu wielokątów (dodałem przykład do zabawy z odpowiedzią Jozuego).
Roman Luštrik,
W przypadku dużych zestawów danych wieloboków przymus w obiekcie SpatialPolygons wiąże się z dużym nakładem pracy i nie jest konieczny. Zastosowanie „over” do SpatialPolygonsDataFrame zwraca indeks wiersza, którego można użyć do podzestawu punktów. Zobacz mój przykład poniżej.
Jeffrey Evans,
Dużo napowietrznych? Zasadniczo pobiera szczelinę @polygons z obiektu SpatialPolygonsDataFrame. Możesz go nawet „sfałszować”, przypisując klasę SpatialPolygonsDataFrame do „SpatialPolygons” (chociaż jest to hacking i nie jest zalecane). Wszystko, co będzie wykorzystywało geometrię, będzie musiało uzyskać ten slot na pewnym etapie, więc mówiąc stosunkowo nie ma w ogóle narzutów. To i tak nie ma znaczenia w żadnej aplikacji w świecie rzeczywistym, w której będziesz przeprowadzał mnóstwo testów wielokątów punktowych.
Spacedman
W rozliczaniu kosztów ogólnych istnieje więcej niż tylko kwestia prędkości. Tworząc nowy obiekt w przestrzeni nazw R, używasz niezbędnej pamięci RAM. Jeśli nie jest to problemem w małych zestawach danych, wpłynie to na wydajność dużych danych. R nie wykazuje liniowej wydajności. Gdy dane stają się coraz bardziej wydajne, potrzeba ding. Jeśli nie musisz tworzyć dodatkowego obiektu, dlaczego miałbyś?
Jeffrey Evans,
1
Nie wiedzieliśmy tego, dopóki tego nie przetestowałem.
Spacedman,
13

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 winietę, aby uzyskać szczegółowe informacje, np. Przypadek wielokrotnego dopasowania nakładających się wielokątów.

Edzer Pebesma
źródło
Dobrze wiedzieć! Nie byłem świadomy metody geometrii.
Jeffrey Evans
2
Witamy na naszej stronie, Edzer - miło cię tu widzieć!
whuber
1
Dzięki Bill - robi się coraz ciszej na stat.ethz.ch/pipermail/r-sig-geo , a może powinniśmy opracować oprogramowanie, które powoduje więcej problemów! ;-)
Edzer Pebesma
6

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) 
Jeffrey Evans
źródło
Źle - nazwy nazw zwracanego obiektu odpowiadają indeksowi wierszy w tym przypadku - generalnie nazwy wierszy wydają się być nazwami wierszy punktów - które mogą nawet nie być liczbowe. Możesz zmodyfikować swoje rozwiązanie, aby dopasować znak, co może uczynić go nieco bardziej niezawodnym.
Spacedman,
@Sapcedman, Nie bądź taki dogmatyczny. Rozwiązanie nie jest nieprawidłowe. Jeśli chcesz podgrupować punkty do zestawu wielokątów lub przypisać wartości wielokątów do punktów, funkcja over działa bez przymusu. Po uzyskaniu wynikowego obiektu można przejść wiele kroków. Rozwiązanie wymuszania na obiekcie SpatialPolygon powoduje znaczny konieczny narzut, ponieważ operację tę można wykonać bezpośrednio na obiekcie SpatialPolygonDataFrame. Przy okazji, zanim edytujesz post, upewnij się, że masz rację. Termin biblioteka i pakiet są używane zamiennie w R.
Jeffrey Evans,
Dodałem kilka testów porównawczych do mojego postu i zauważyłem inny problem z twoją funkcją. Również „Pakiety to zbiory funkcji R, danych i skompilowanego kodu w dobrze zdefiniowanym formacie. Katalog, w którym przechowywane są pakiety, nazywa się biblioteką”
Spacedman,
Chociaż technicznie masz rację, jeśli chodzi o „pakiet” kontra „biblioteka”, argumentujesz semantyką. Właśnie otrzymałem prośbę redaktora ds. Modelowania ekologicznego, abyśmy zmienili użycie „pakietu” (co jest właściwie moją preferencją) na „bibliotekę”. Chodzi mi o to, że stają się one wymiennymi terminami i kwestią preferencji.
Jeffrey Evans,
1
„technicznie poprawne”, jak zauważył kiedyś dr Sheldon Cooper, „jest najlepszym rodzajem poprawności”. Ten edytor jest technicznie zły, co jest najgorszym rodzajem błędu.
Spacedman,
4

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)
Josh O'Brien
źródło
Nie plydziała okropnie, jeśli ma więcej niż jedną funkcję, ponieważ gIntersectszwraca macierz z jednym wierszem dla każdej cechy. Prawdopodobnie możesz zamienić wiersze na wartość PRAWDA.
Spacedman,
@Spacedman - Bingo. Musisz to zrobić apply(gIntersects(pts, ply, byid=TRUE), 2, any). Właściwie przejdę do odpowiedzi i przestawię odpowiedź na to, ponieważ obejmuje to również przypadek pojedynczego wielokąta.
Josh O'Brien,
Ach any. To może być nieznacznie szybsze niż wersja, którą właśnie przetestowałem.
Spacedman,
@Spacedman - Z moich szybkich testów wygląda obrieni rowlings2biegnie szyja z szyją, obrien może z 2% szybszym.
Josh O'Brien,
@ JoshO'Brien, jak można użyć tej odpowiedzi na wielu wielokątach? To pppowinno mieć znak IDwskazujący, w którym wielokącie znajdują się punkty.
code123
4

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
fdetsch
źródło
1
Powinno tmpbyć pts.intersect? Także parsowanie zwróconych nazw jest zależne od nieudokumentowanego zachowania.
Spacedman,
@Spacedman, masz rację tmp, zapomniałeś go usunąć po zakończeniu kodu. Masz rację także podczas analizowania dimnames. To było szybkie rozwiązanie, które zapewniło pytającemu szybką odpowiedź, a na pewno są lepsze (i bardziej uniwersalne) podejścia, na przykład twoje :-)
fdetsch
1

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
rafa.pereira
źródło