使用 R 'ks' 包在重叠的 kdes 中提取数据点

Question

我正在使用 R 中的 ks 包，并且想要确定哪些位置数据点落在重叠的 2d 核轮廓区域内（我正在比较两个不同物种的家庭范围的 UD）。下面有一个示例（修改自：http://www.rdocumentation.org/packages/ks/versions/1.5.3/topics/plot.kde?）。

我要生成的是落在 fhatx 轮廓内的所有 y 点的列表（例如黑线内的黄点）。反之亦然，我想要一个落在等高线内的 x 坐标列表。

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)

Answer 1

kde 的输出可以转换为栅格，然后您可以使用 rasterToPolygons 函数从那里提取任何轮廓。一旦您的点被转换为 sp 包可识别的格式，您就可以使用 gIntersection 函数查看空间对象之间的任何类型的交集。

您最终得到两个 SpatialPoints 对象 x.inY 和 y.inX，它们包含 fhaty 的 50% 轮廓内包含的 x 点，反之亦然。这些点的坐标可以使用 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)

使用 R 'ks' 包在重叠的 kdes 中提取数据点

Extracting data points within overlapping kdes using R 'ks' package

r

polygon

kde-plasma