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