แยกจุดข้อมูลภายใน kdes ที่ทับซ้อนกันโดยใช้แพ็คเกจ R 'ks'

ผลลัพธ์ของ kde สามารถแปลงเป็นแรสเตอร์ได้ จากนั้นคุณสามารถแยกรูปร่างใดก็ได้โดยใช้ฟังก์ชัน rasterToPolygons เมื่อคะแนนของคุณถูกแปลงเป็นรูปแบบที่แพ็คเกจ sp รู้จัก คุณจะสามารถดูจุดตัดระหว่างวัตถุเชิงพื้นที่ได้โดยใช้ฟังก์ชัน gIntersection

คุณจะพบกับวัตถุ SpatialPoints สองอัน x.inY และ y.inX ซึ่งมีจุด x ที่รวมอยู่ในรูปร่าง 50% ของ fhaty และในทางกลับกัน พิกัดของจุดเหล่านี้อาจถูกแยกออกมาในอาร์เรย์โดยใช้ coordinates(...)

นี่อาจไม่ใช่วิธีแก้ปัญหาที่หรูหราที่สุด แต่ควรจะทำงานได้ดีหากอาร์เรย์ที่ออกโดยฟังก์ชัน kde นั้นไม่ใหญ่เกินไป

ฉันหวังว่านี่จะช่วยได้.

ขั้นตอนที่ 1: แปลงเอาต์พุตจาก kde เป็นวัตถุแรสเตอร์

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

ขั้นตอนที่ 2: ปรับสเกลแรสเตอร์ระหว่าง 0 ถึง 100 จากนั้นแปลงเซลล์ทั้งหมดภายในรูปร่าง 50 เป็น 1 แน่นอนว่ารูปร่างอาจเปลี่ยนเป็น 95 หรือค่าอื่น

#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

ขั้นตอนที่ 3: แยกรูปหลายเหลี่ยมที่สอดคล้องกับรูปร่าง 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)

ขั้นตอนที่ 4: แปลงคะแนนของคุณเป็นวัตถุอวกาศเพื่อใช้ไลบรารี sp

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

ขั้นตอนที่ 5: ค้นหาจุดที่ตัดกับรูปหลายเหลี่ยมหนึ่งหรืออีกรูปหนึ่ง

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

โครงเรื่อง

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)