如何在R中排除超出一条直线的点

4
我正在使用R处理两个数据框:一个是"red"数据框,另一个是"black"数据框。它们都有两列表示坐标。
我用一张图来解释我的意思。

我想选择所有在"黑色"线之外的"红色"数据框中的点。例如,所有在由黑色点形成的多边形区域之外的点。

enter image description here

5个回答

1
您可以使用sf包定义凸包,然后将目标点与该多边形相交。
基于black定义一个凸包:
library(sf)

set.seed(99)
red <- data.frame(x = runif(100,-10,10), y = runif(100,-4,4))
black <- data.frame(x = runif(100,-8,8), y = runif(100,-4,3))

# Convert df to point feature
blk <- st_as_sf(black, coords = c("x", "y"))
# Convert to multipoint
blk_mp <- st_combine(blk)
# Define convex hull
blk_poly <- st_convex_hull(blk_mp)

plot(black)
points(red, col = "red")
plot(blk_poly, add = TRUE)

进入图片描述

将凸包与red相交,返回该多边形中的red

rd <- st_as_sf(red, coords = c("x", "y"))
rd_inside <- st_intersection(rd, blk_poly)

plot(black)
points(red)
plot(blk_poly, add = TRUE)
plot(rd_inside, pch = 24, col = "red", bg = "red", add = TRUE)

enter image description here


1
“黑线”对应于一个凹多边形。 - user1196549
我认为使用@Skaqqs建议的凸包可能是正确的方法,问题是我需要尽可能精确地绘制追踪轮廓的线条。实际上,凸包对我不起作用,因为我将在图的左上角形成对角线中包含所有红色点的凸包中,无法将它们排除并检索到它们。 - Giordano
@Giordano,你能分享一下你在问题中使用的数据吗?如果你的点是有序的,我们可以很容易地将它们转换成一个多边形。 - Skaqqs

1
一种可能的解决方案是在点之后绘制多边形,并用白色填充其外部区域。这不能直接使用 polygonpolypath 完成,因为这些函数只能填充多边形的内部。然而,您可以使用 polypath 填充两个多边形之间的区域。因此,您可以添加第二个多边形,该多边形包含(或超过)绘图边界。
以下是在基础 R 中有效的示例:
p.outer <- list(x=c(0,100,100,0), y=c(0,0,100,100))
p.inner <- list(x=c(20,40,80,50,40,30), y=c(30,20,70,80,50,60))
plot(p.outer, type="n")
points(runif(100, min=0, max=100), runif(100, min=0, max=100))
polypath(x=c(p.outer$x, NA, p.inner$x), y = c(p.outer$y, NA, p.inner$y), col ="white", rule="evenodd")

enter image description here


1

我的上一个回答只展示了如何不要在多边形外部绘制点。为了实际上确定多边形外部的点,你可以使用包ptinpoly中的函数pip2d。它返回负值表示点在多边形外部。

例如:

library(ptinpoly)
poly.vertices <- data.frame(x=c(20,40,80,50,40,30), y=c(30,20,70,80,50,60))
p <- data.frame(x=runif(100, min=0, max=100), y=runif(100, min=0, max=100))

outside <- (pip2d(as.matrix(poly.vertices), as.matrix(p)) < 0)

plot(p$x, p$y, col=ifelse(outside, "red", "black"))
polygon(poly.vertices$x, poly.vertices$y, border="blue", col=NA)

使用包DescTools中的函数PtInPoly也可以实现相同的功能,对于多边形外部的点,该函数返回零。然而,ptinpoly的实现具有优势,因为它实现了特别高效的算法,该算法在以下文献中得到描述:

J. Liu, Y.Q. Chen, J.M. Maisog, G. Luta: "A new point containment test algorithm based on preprocessing and determining triangles." Computer-Aided Design, Volume 42, Issue 12, December 2010, Pages 1143-1150

编辑:出于好奇,我使用microbenchmark比较了ptinpoly::pip2dDescTools::PtInPoly的运行时间,并且使用N=50000个点进行测试,结果表明pip2d要快得多:

> microbenchmark(outside.pip2d(), outside.PtInPoly())
Unit: milliseconds
              expr       min        lq      mean   median        uq      max
   outside.pip2d()  3.375084  3.421631  4.459051  3.48939  4.251395 65.97793
outside.PtInPoly() 27.537927 27.666688 28.739288 27.97984 28.514595 90.11313

neval
100
100

在此输入图片描述


非常好的解决方案!我尝试了一下,但似乎计算上有些具有挑战性。我认为我需要一个基于数据框连接或过滤的解决方案。 - Giordano
1
你尝试过使用DescTools::PtInPolyptinpoly::pip2d两种方法,但速度都太慢了吗?作为预处理步骤,你可以移除所有明显在多边形外部的点(即在多边形外部的最大边界框之外的所有点)和明显在多边形内部的点(即在多边形内部的最大矩形之内的所有点),但这可能仍然会留下太多的点。 - cdalitz

1

如@cdalitz所建议的那样,我使用DescTools包中的PtInPoly函数解决了问题。该函数返回点的坐标(在我的情况下是红色坐标)的数据框和第三列“pip”,其中1表示该点在多边形内,0表示在多边形外。

我将使用另一个数据集来展示结果:

 try <- DescTools::PtInPoly(pnts = red[,c("x","y")], poly.pnts = black[,c("x","y")])

 ggplot()+
      geom_point(try, mapping = aes(x = x, y = y, color = as.character(pip))) +
      geom_polygon(data = black, mapping = aes(x,y))

enter image description here


正如您在上面评论的那样,ptinpoly::pip2d 对您的需求来说太慢了,我想知道为什么 DescTools::PtInPoly 更快。在我的测试中,情况恰恰相反。 - cdalitz
我其实不知道。使用ptinpoly时,我的R崩溃了。我以为是计算量的原因。 - Giordano

0

看起来你可以通过将每个点连接到其最近的邻居和相反方向上的最近邻居来重构黑色多边形的边缘。然后进行点在多边形内测试。


听起来是个有趣的想法,你知道如何在R中实现吗? - Giordano
@Giordano: 抱歉,我不会。 :-) - user1196549
@giordano,也许你可以将这个翻译成R语言?我可以提供帮助。https://dev59.com/uGgv5IYBdhLWcg3wHNDG - Skaqqs
@Skaqqs:这个方法不同,不能在这里使用。 - user1196549

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