Извлечение точек данных в перекрывающихся kdes с использованием пакета R 'ks'

Я использую пакет ks в R и хочу определить, какие точки данных о местоположении попадают в области перекрывающихся двухмерных контуров ядра (я сравниваю UD двух разных ареалов обитания видов). Ниже приведен пример (измененный с: http://www.rdocumentation.org/packages/ks/versions/1.5.3/topics/plot.kde?).

Я пытаюсь создать список всех точек y, которые попадают в контуры fhatx (например, желтые точки внутри черных линий). И наоборот, мне нужен список координат x, попадающих в контурные линии fhaty.

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 источник
comment
Вы должны предоставить своего рода воспроизводимый пример с примерами данных, которые можно использовать для проверки возможных решений. Было бы легче помочь вам и сделать более ясным, каков желаемый результат.   -  person MrFlick    schedule 02.09.2016
comment
Вы правы, спасибо. Я обновил его воспроизводимым примером.   -  person Erin    schedule 03.09.2016

Ответы (1)


arrow_upward
0
arrow_downward

Выходные данные 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)
person S.Derville    schedule 05.09.2016
comment
идеально. Спасибо. - person Erin; 07.09.2016