R - 如何在特定轮廓内找到点

Posted 2023-03-15

【中文标题】R - 如何在特定轮廓内找到点

【英文标题】：R - How to find points within specific Contour

【发布时间】：2015-08-11 14:05:24

【问题描述】：

我正在用 kde2d (MASS) 在 lat 和 lon 数据上创建密度图。我想知道原始数据中的哪些点在特定轮廓内。

我使用两种方法创建 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)

这是第 2 版 - 使用 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% 轮廓内的精确点。

我试图找到与每行 lat 和 lon 值相关联的 z 值（来自 kde2d 的密度图的高度），但没有运气。我还想我可以在数据中添加一个 ID 列来标记每一行，然后在使用 melt() 之后以某种方式将其转移过来。然后我可以简单地将 z 值与我想要的每个轮廓匹配的数据子集，并根据 ID 列将它们与原始 BW 数据进行比较。

这是我正在谈论的图片：

我想知道哪些红点在 50% 等高线（蓝色）内，哪些在 90% 等高线内（红色）。

注意：此代码大部分来自其他问题。向所有做出贡献的人大声喊叫！

谢谢！

【问题讨论】：

当您说“在 90% 和 50% 等高线内”时，您的意思是您想知道 z 值大于 90% 或 50% 的所有点的纬度/经度所有的 z 值？

已编辑问题 - 我想找到 2 个轮廓“圆圈”内的红点。

【参考方案1】：

你可以从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)

【讨论】：

太棒了！简单明了。现在使用 point.in.polygon 的答案非常明显。信息量超级大。

@jenesaisquoi，如果我想使用代码来查找一对新点是否在 95% 的轮廓内，我需要做什么？

【参考方案2】：

我认为这是我能想到的最好方法。它使用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)

【讨论】：

太棒了！感谢您提供所有其他信息。我将不得不接受 6pool 的回答，因为它更直接一些。但是，您的回答让我对各种新的可能性有了很多了解！ :)

嗨，我正在尝试复制上面的代码。有人可以解释这是在做什么吗？ dx

代码正在提取轮廓网格级别中与prob 向量中提供的值相对应的点。查看kde2d() 函数的文档和dens 的数据结构以了解发生了什么。基本上，您正在查看 X/Y 方向上的差分向量和 Z 值的累积和，以找到对应于适当百分位数的网格值。

那么，如果我想获得在 90% 等高线内但不在 50% 等高线内的点，应该 out-inner 给出结果吗？

我有点困惑...over() 函数为您提供所需的一切？要计算某个范围内的这些点（比如在 0.5 和 0.9 等值线之间），您可以执行以下操作：pts[!is.na( over(pts, polys[1]) ) & is.na( over(pts, polys[2]) )] 希望我能理解您的问题？