R - Cách tìm các điểm trong đường viền cụ thể

Question

Tôi đang tạo các ô mật độ với kde2d (MASS) trên dữ liệu lat và lon. Tôi muốn biết điểm nào từ dữ liệu gốc nằm trong một đường bao cụ thể.

Tôi tạo đường viền 90% và 50% bằng hai cách tiếp cận. Tôi muốn biết điểm nào nằm trong giới hạn 90% và điểm nào nằm trong đường viền 50%. Các điểm trong đường bao 90% sẽ chứa tất cả các điểm trong đường viền 50%. Bước cuối cùng là tìm các điểm trong đường viền 90% không nằm trong đường bao 50% (tôi không nhất thiết cần trợ giúp với bước này).

# bw = data of 2 cols (lat and lon) and 363 rows 
# two versions to do this: 
# would ideally like to use the second version (with ggplot2) 

# version 1 (without ggplot2) 
library(MASS) 
x <- bw$lon 
y <- bw$lat 
dens <- kde2d(x, y, n=200) 

# the contours to plot 
prob <- c(0.9, 0.5) 
dx <- diff(dens$x[1:2]) 
dy <- diff(dens$y[1:2]) 
sz <- sort(dens$z) 
c1 <- cumsum(sz) * dx * dy 
levels <- sapply(prob, function(x) { 
    approx(c1, sz, xout = 1 - x)$y 
}) 
plot(x,y) 
contour(dens, levels=levels, labels=prob, add=T) 

Và đây là phiên bản 2 - sử dụng ggplot2. Tôi lý tưởng muốn sử dụng phiên bản này để tìm các điểm trong vòng 90% và 50% đường nét.

# version 2 (with ggplot2) 
getLevel <- function(x,y,prob) { 
    kk <- MASS::kde2d(x,y) 
    dx <- diff(kk$x[1:2]) 
    dy <- diff(kk$y[1:2]) 
    sz <- sort(kk$z) 
    c1 <- cumsum(sz) * dx * dy 
    approx(c1, sz, xout = 1 - prob)$y 
} 

# 90 and 50% contours 
L90 <- getLevel(bw$lon, bw$lat, 0.9) 
L50 <- getLevel(bw$lon, bw$lat, 0.5) 

kk <- MASS::kde2d(bw$lon, bw$lat) 
dimnames(kk$z) <- list(kk$x, kk$y) 
dc <- melt(kk$z) 

p <- ggplot(dc, aes(x=Var1, y=Var2)) + geom_tile(aes(fill=value)) 
+ geom_contour(aes(z=value), breaks=L90, colour="red") 
+ geom_contour(aes(z=value), breaks=L50, color="yellow") 
+ ggtitle("90 (red) and 50 (yellow) contours of BW") 

Tôi tạo các ô với tất cả các điểm vĩ độ và điểm được vẽ và đường nét 90% và 50%. Tôi chỉ đơn giản muốn biết cách trích xuất các điểm chính xác nằm trong các đường bao 90% và 50%.

Tôi đã cố gắng tìm các giá trị z (độ cao của các ô mật độ từ kde2d) được liên kết với mỗi hàng giá trị lat và lon nhưng không có may mắn. Tôi cũng nghĩ rằng tôi có thể thêm một cột ID vào dữ liệu để gắn nhãn mỗi hàng và sau đó bằng cách nào đó chuyển nó qua sau khi sử dụng melt(). Sau đó, tôi có thể chỉ đơn giản là tập hợp các dữ liệu có giá trị của z phù hợp với mỗi đường viền tôi muốn và xem lat và lon nào được so sánh với dữ liệu BW ban đầu dựa trên cột ID.

Dưới đây là một bức tranh về những gì tôi đang nói về:

Tôi muốn biết điểm đỏ nằm trong đường viền 50% (màu xanh) và đó là trong đường viền 90% (màu đỏ).

Lưu ý: phần lớn mã này là từ các câu hỏi khác. Big hét lên cho tất cả những người đã đóng góp!

Cảm ơn bạn!

Answer 1

Bạn có thể sử dụng point.in.polygon từ sp

## Interactively check points 
plot(bw) 
identify(bw$lon, bw$lat, labels=paste("(", round(bw$lon,2), ",", round(bw$lat,2), ")")) 

## Points within polygons 
library(sp) 
dens <- kde2d(x, y, n=200, lims=c(c(-73, -70), c(-13, -12))) # don't clip the contour 
ls <- contourLines(dens, level=levels) 
inner <- point.in.polygon(bw$lon, bw$lat, ls[[2]]$x, ls[[2]]$y) 
out <- point.in.polygon(bw$lon, bw$lat, ls[[1]]$x, ls[[1]]$y) 

## Plot 
bw$region <- factor(inner + out) 
plot(lat ~ lon, col=region, data=bw, pch=15) 
contour(dens, levels=levels, labels=prob, add=T) 

Answer 2

Tôi nghĩ đây là cách tốt nhất tôi có thể nghĩ đến. Điều này sử dụng một mẹo để chuyển đổi các đường đồng mức thành SpatialLinesDataFrame các đối tượng sử dụng chức năng ContourLines2SLDF() từ gói maptools. Sau đó, tôi sử dụng một mẹo được phác thảo trong Bivand, et al.'s Phân tích dữ liệu không gian được áp dụng với R để chuyển đổi đối tượng SpatialLinesDataFrame thành SpatialPolygons. Những sau đó có thể được sử dụng với over() chức năng để trích xuất điểm trong mỗi đa giác đường viền:

## Simulate some lat/lon data: 
x <- rnorm(363, 45, 10) 
y <- rnorm(363, 45, 10) 

## Version 1 (without ggplot2): 
library(MASS) 
dens <- kde2d(x, y, n=200) 

## The contours to plot: 
prob <- c(0.9, 0.5) 
dx <- diff(dens$x[1:2]) 
dy <- diff(dens$y[1:2]) 
sz <- sort(dens$z) 
c1 <- cumsum(sz) * dx * dy 
levels <- sapply(prob, function(x) { 
    approx(c1, sz, xout = 1 - x)$y 
}) 
plot(x,y) 
contour(dens, levels=levels, labels=prob, add=T) 

## Create spatial objects: 
library(sp) 
library(maptools) 

pts <- SpatialPoints(cbind(x,y)) 

lines <- ContourLines2SLDF(contourLines(dens, levels=levels)) 

## Convert SpatialLinesDataFrame to SpatialPolygons: 
lns <- slot(lines, "lines") 
polys <- SpatialPolygons(lapply(lns, function(x) { 
    Polygons(list(Polygon(slot(slot(x, "Lines")[[1]], 
    "coords"))), ID=slot(x, "ID")) 
    })) 

## Construct plot from your points, 
plot(pts) 

## Plot points within contours by using the over() function: 
points(pts[!is.na(over(pts, polys[1]))], col="red", pch=20) 
points(pts[!is.na(over(pts, polys[2]))], col="blue", pch=20) 

contour(dens, levels=levels, labels=prob, add=T)