bagaimana cara overlay poligon di atas SpatialPointsDataFrame dan melestarikan data SPDF?

Saya punya SpatialPointsDataFramedengan beberapa data tambahan. Saya ingin mengekstrak titik-titik itu di dalam poligon dan pada saat yang sama, menyimpan SPDFobjek dan data terkait.

Sejauh ini saya kurang beruntung dan terpaksa melakukan pencocokan dan penggabungan melalui ID umum, tetapi ini hanya berfungsi karena saya memiliki data grid dengan masing-masing IDS.

Ini contoh cepat, saya mencari poin di dalam kotak merah.

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

Pendekatan yang paling jelas adalah menggunakan over, tetapi ini mengembalikan data dari poligon.

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

Dari sp::overbantuan:

 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.

Jadi, jika Anda mengonversi SpatialPolygonsDataFrameke SpatialPolygonsAnda mendapatkan kembali vektor indeks dan Anda dapat mengelompokkan poin Anda pada 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
> 

Bagi yang ragu, inilah bukti bahwa overhead konversi bukan masalah:

Dua fungsi - metode pertama Jeffrey Evans, lalu yang asli, lalu konversi yang diretas, kemudian versi berdasarkan gIntersectsjawaban 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,]
}

Sekarang untuk contoh dunia nyata, saya telah menyebarkan beberapa poin acak ke columbusset data:

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

Kelihatan bagus

plot(columbus)
points(pts)

Periksa fungsi melakukan hal yang sama:

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

Dan jalankan 500 kali untuk pembandingan:

> 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 

Sesuai intuisi saya, ini bukan overhead yang bagus, sebenarnya mungkin lebih sedikit overhead daripada mengubah semua indeks baris menjadi karakter dan kembali, atau menjalankan na.omit untuk mendapatkan nilai yang hilang. Yang kebetulan mengarah ke mode kegagalan evansfungsi lainnya ...

Jika satu baris bingkai data poligon adalah semua NA(yang benar-benar valid), maka hamparan dengan SpatialPolygonsDataFrameuntuk titik-titik dalam poligon tersebut akan menghasilkan bingkai data keluaran dengan semuaNA s, yang evans()kemudian akan turun:

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

TAPI gIntersects lebih cepat, bahkan dengan harus menyapu matriks untuk memeriksa persimpangan di R daripada dalam kode C. Saya menduga ini adalah prepared geometryketerampilan GEOS, membuat indeks spasial - ya, dengan prepared=FALSEitu membutuhkan waktu sedikit lebih lama, sekitar 5,5 detik.

Saya terkejut tidak ada fungsi untuk langsung mengembalikan indeks atau poin. Ketika saya menulis splancs20 tahun yang lalu fungsi point-in-polygon memiliki keduanya ...

sp menyediakan formulir yang lebih pendek untuk memilih fitur berdasarkan persimpangan spasial, mengikuti contoh OP:

pts[ply,]

pada:

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

Di balik layar ini adalah kependekan dari

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

Yang perlu diperhatikan adalah bahwa ada geometrymetode yang menjatuhkan atribut: overmengubah perilaku jika argumen keduanya memiliki atribut atau tidak (ini adalah kebingungan OP). Ini berfungsi di semua kelas Spasial * sp, meskipun beberapa overmetode memerlukan rgeos, lihat sketsa ini untuk detail, misalnya kasus beberapa kecocokan untuk poligon yang tumpang tindih.

Anda berada di jalur yang benar dengan berakhir. Rownames dari objek yang dikembalikan sesuai dengan indeks baris poin. Anda dapat menerapkan pendekatan tepat Anda hanya dengan beberapa baris kode tambahan.

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) 

Apakah ini yang Anda cari?

Satu catatan, saat diedit: Panggilan pembungkus ke apply()diperlukan untuk membuat pekerjaan ini dengan SpatialPolygonsobjek sewenang-wenang , mungkin berisi lebih dari satu fitur poligon. Terima kasih kepada @Spacedman karena mendorong saya untuk mendemonstrasikan bagaimana menerapkan ini pada kasus yang lebih umum.

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)

Berikut ini pendekatan yang mungkin menggunakan rgeospaket. Pada dasarnya, ini menggunakan gIntersectionfungsi yang memungkinkan Anda untuk memotong dua spobjek. Dengan mengekstraksi ID dari titik-titik yang terletak di dalam poligon, Anda kemudian dapat mengelompokkan bagian asli Anda SpatialPointsDataFrame, menyimpan semua data yang sesuai. Kode ini hampir menjelaskan sendiri, tetapi jika ada pertanyaan, jangan ragu untuk bertanya!

# 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

Ada solusi yang sangat sederhana menggunakan spatialEcoperpustakaan.

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)

Periksa hasilnya:

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