나는이 SpatialPointsDataFrame몇 가지 추가 데이터. 다각형 내에서 해당 점을 추출하고 동시에 SPDF객체와 해당 데이터를 유지하고 싶습니다 .
지금까지 운이 거의 없었고 공통 ID를 통해 일치하고 병합하는 데 의존했지만 개별 ID로 데이터를 그리드 화했기 때문에 작동합니다.
다음은 간단한 예입니다. 빨간색 사각형 안에있는 점을 찾고 있습니다.
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답변
로부터 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
> 의심스러운 사람들에게는 변환 오버 헤드가 문제가되지 않는다는 증거가 있습니다.
두 가지 기능-먼저 Jeffrey Evans의 방법, 내 원본, 해킹 된 변환, gIntersectsJosh 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,]
}실제 예제에서는 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
> 그러나 gIntersectsC 코드가 아닌 R의 교차점을 검사하기 위해 매트릭스를 스위프해야하더라도 더 빠릅니다. 나는 그것의 의심prepared geometry공간 인덱스를 생성하는 GEOS 기술을 예 prepared=FALSE, 약 5.5 초가 조금 더 걸립니다.
지수 나 포인트를 직선으로 돌려주는 함수가 없다는 것에 놀랐습니다. 내가 splancs20 년 전에 쓰면서 다각형 함수는 둘 다 가지고있었습니다.
답변
sp OP 예에 따라 공간 교차를 기반으로 피쳐를 선택하는 짧은 형식을 제공합니다.
pts[ply,]현재:
points(pts[ply,], col = 'red')무대 뒤에서 이것은 짧은
pts[!is.na(over(pts, geometry(ply))),]주목해야 할 것은 geometry속성을 삭제 하는 메소드 가 있다는 것입니다. over두 번째 인수에 속성이 있는지 여부에 따라 동작이 변경됩니다 (이것은 OP의 혼란이었습니다). 이 방법은의 모든 Spatial * 클래스에서 작동 sp하지만, 일부 over메소드에는 필요 하지만 , 겹치는 다각형에 대해 여러 개의 일치하는 경우와 같은 자세한 내용 rgeos은 이 비네팅을 참조하십시오 .
답변
당신은 올바른 길을 가고있었습니다. 반환 된 객체의 행 이름은 포인트의 행 인덱스에 해당합니다. 몇 줄의 코드 만 추가하면 정확한 접근 방식을 구현할 수 있습니다.
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) 답변
이것이 당신이 무엇을하고 있습니까?
편집시 참고 사항 : 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)답변
rgeos패키지를 사용하는 가능한 방법은 다음과 같습니다 . 기본적으로 gIntersection두 sp객체 를 교차시킬 수 있는 기능을 사용 합니다. 다각형 내에있는 점의 ID를 추출하면 이후 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답변
라이브러리를 사용 하는 매우 간단한 솔루션 이 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