Mengekstraksi titik data dalam kdes yang tumpang tindih menggunakan paket R 'ks'

Saya menggunakan paket ks di R, dan ingin menentukan titik data lokasi mana yang termasuk dalam area kontur kernel 2d yang tumpang tindih (saya membandingkan UD dari dua wilayah jelajah spesies yang berbeda). Ada contoh di bawah ini (dimodifikasi dari: http://www.rdocumentation.org/packages/ks/versions/1.5.3/topics/plot.kde?).

Apa yang saya coba buat adalah daftar semua titik y yang termasuk dalam kontur fhatx (misalnya titik kuning di dalam garis hitam). Dan sebaliknya, saya ingin daftar koordinat x yang berada dalam garis kontur yang tepat.

library(ks)
x <- rmvnorm.mixt(n=100, mus=c(0,0), Sigmas=diag(2), props=1)
Hx <- Hpi(x)
fhatx <- kde(x=x, H=Hx) 
y <- rmvnorm.mixt(n=100, mus=c(0.5,0.5), Sigmas=0.5*diag(2), props=1)
Hy <- Hpi(y)
fhaty <- kde(x=y, H=Hy)
contourLevels(fhatx, cont=c(75, 50, 25))
contourSizes(fhatx, cont=25, approx=TRUE)
plot(fhatx, cont=c(50,95), drawpoints=TRUE)
plot(fhaty, cont=c(50,95), col=3, drawpoints=TRUE,col.pt="yellow", add=TRUE)

r kde polygon
person Erin    schedule 02.09.2016    source sumber
comment
Anda harus memberikan semacam contoh yang dapat direproduksi dengan data sampel yang dapat digunakan untuk menguji solusi yang mungkin. Akan lebih mudah untuk membantu Anda dan memperjelas apa hasil yang diinginkan.   -  person MrFlick    schedule 02.09.2016
comment
Anda benar, terima kasih. Saya telah memperbaruinya dengan contoh yang dapat direproduksi.   -  person Erin    schedule 03.09.2016

Jawaban (1)


arrow_upward
0
arrow_downward

Output dari kde dapat diubah menjadi raster dan kemudian dari sana Anda dapat mengekstrak kontur apa pun menggunakan fungsi rasterToPolygons. Setelah titik Anda dikonversi ke format yang dikenali oleh paket sp, Anda dapat melihat perpotongan apa pun antara objek spasial menggunakan fungsi gIntersection.

Anda mendapatkan dua objek SpatialPoints x.inY dan y.inX yang berisi titik x yang termasuk dalam kontur 50% dari fhaty dan sebaliknya. Koordinat titik-titik ini dapat diekstraksi dalam array menggunakan coordinates(...).

Ini mungkin bukan solusi yang paling elegan tetapi akan berfungsi dengan baik jika array yang dikeluarkan oleh fungsi kde tidak terlalu besar.

Semoga membantu.

LANGKAH 1: ubah output dari kde menjadi objek raster

# for the x kde
arrayX <- expand.grid(list(fhatx$eval.points[[1]],fhatx$eval.points[[2]]))
arrayX$z <- as.vector(fhatx$estimate)
rasterX <- rasterFromXYZ(arrayX)
# for the y kde
arrayY <- expand.grid(list(fhaty$eval.points[[1]],fhaty$eval.points[[2]]))
arrayY$z <- as.vector(fhaty$estimate)
rasterY <- rasterFromXYZ(arrayY)

LANGKAH 2: ubah skala raster antara 0 dan 100, lalu ubah semua sel dalam kontur 50 menjadi 1. Tentu saja kontur dapat diubah menjadi 95 atau nilai lainnya

#for raster x
rasterX <- rasterX*100/rasterX@data@max
rasterX[rasterX[]<=50,] <- NA
rasterX[rasterX[]>50,] <- 1
#[enter image description here][1]for raster y
rasterY <- rasterY*100/rasterY@data@max
rasterY[rasterY[]<=50,] <- NA
rasterY[rasterY[]>50,] <- 1

LANGKAH 3: ekstrak poligon yang sesuai dengan kontur 50%.

polyX50<-rasterToPolygons(rasterX, n=16, na.rm=T, digits=4, dissolve=T)
polyY50<-rasterToPolygons(rasterY, n=16, na.rm=T, digits=4, dissolve=T)

LANGKAH 4: ubah poin Anda menjadi objek spasial untuk menggunakan perpustakaan sp

x.points <- SpatialPoints(x)
y.points <- SpatialPoints(y)

LANGKAH 5: temukan titik-titik yang berpotongan dengan satu poligon atau poligon lainnya

#x points falling in fhatx 50 contour
x.inY <- gIntersection(x.points, polyY50)
#y points falling in fhatx 50 contour
y.inX <- gIntersection(y.points, polyX50)

Merencanakan

par(mfrow=c(1,2))
plot(fhatx, cont=c(50,95), col="red")
plot(fhaty, cont=c(50,95), col="green",add=TRUE)
plot(x.points, col="red", add=T)
plot(y.points, col="green", add=T)

plot(fhatx, cont=c(50,95), col="red")
plot(fhaty, cont=c(50,95), col="green",add=TRUE)
plot(x.inY, col="red", add=T)
plot(y.inX, col="green", add=T)
person S.Derville    schedule 05.09.2016
comment
sempurna. Terima kasih. - person Erin; 07.09.2016