R - 如何在指定轮廓内查找点

16
我正在使用kde2d (MASS) 对经纬度数据创建密度图。我想知道原始数据中哪些点在特定轮廓线内。
我使用两种方法创建90%和50%的轮廓线。我想知道哪些点在90%轮廓线内,哪些点在50%轮廓线内。90%轮廓线内的点将包含所有50%轮廓线内的点。最后一步是找到在90%轮廓线内但不在50%轮廓线内的点(我不一定需要这一步的帮助)。
# 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)

这里是第二个版本 - 使用 ggplot2。我希望使用这个版本来找到在90%和50%轮廓线内的点。

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

我创建了包含所有纬度和经度点以及90%和50%轮廓的绘图。我只想知道如何提取在90%和50%轮廓内的确切点。
我尝试查找与每行纬度和经度值相关联的z值(从kde2d的密度图中获取的高程),但没有成功。我还考虑过为数据添加一个ID列来标记每一行,然后在使用melt()之后以某种方式将其转移过去。然后,我可以简单地对具有与我想要的每个轮廓匹配的z值的数据进行子集,并根据ID列比较它们与原始BW数据中的哪些纬度和经度。
这是我所说的图片:
我想知道哪些红色点在50%轮廓(蓝色)内,哪些在90%轮廓(红色)内。
注意:很多代码都来自其他问题。向所有做出贡献的人致以最大的赞扬!

当你说“在90%和50%轮廓线内”时,你是指想知道所有z值大于所有z值的90%或50%的点的纬度/经度吗? - eipi10
我想找到在两个轮廓“圆圈”内的红点。 - squishy
2个回答

14
你可以使用 `sp` 中的 `point.in.polygon` 函数。
## 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)

enter image description here


太棒了!简单明了,而且点在多边形内这个答案现在变得非常显而易见了。超级有信息量。 - squishy
@jenesaisquoi,如果我想使用该代码查找新点对是否落在95%轮廓内,我需要做什么? - user1560215

5
我认为这是我能想到的最好方法。它使用一个技巧,使用maptools软件包中的ContourLines2SLDF()函数将等高线转换为SpatialLinesDataFrame对象。然后,我使用Bivand等人在Applied Spatial Data Analysis with R中概述的技巧,将SpatialLinesDataFrame对象转换为SpatialPolygons。 然后可以使用over()函数与每个等高线多边形中的点进行提取。
##  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)

enter image description here


太棒了!感谢您提供的所有额外信息。我将不得不接受6pool的答案,因为它更加直接。但是,您的答案让我对各种新可能性有了很多深入的了解! :) - squishy
嗨,我正在尝试复制上面的代码。有人能解释一下这是在做什么吗?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 }) - user1560215
该代码正在提取与“prob”向量中提供的值相对应的轮廓网格级别中的点。查看“kde2d()”函数的文档和“dens”的数据结构,以了解发生了什么。基本上,您需要查看X / Y方向的差异向量和Z值的累积总和,以找到与适当百分位数相对应的网格值。 - Forrest R. Stevens
那么,如果我想获取位于90%等高线但不在50%等高线上的点,out-inner 应该提供结果吗? - user1560215
我有点困惑... over() 函数是否给你需要的一切?为了计算在某个区间内的点(比如 0.5 和 0.9 等高线之间的点),你可以采取以下方法:pts[!is.na( over(pts, polys[1]) ) & is.na( over(pts, polys[2]) )]。希望我理解了你的问题。 - Forrest R. Stevens

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接