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)
polygon
或 polypath
完成,因为这些函数只能填充多边形的内部。然而,您可以使用 polypath
填充两个多边形之间的区域。因此,您可以添加第二个多边形,该多边形包含(或超过)绘图边界。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")
我的上一个回答只展示了如何不要在多边形外部绘制点。为了实际上确定多边形外部的点,你可以使用包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::pip2d
和DescTools::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
DescTools::PtInPoly
和ptinpoly::pip2d
两种方法,但速度都太慢了吗?作为预处理步骤,你可以移除所有明显在多边形外部的点(即在多边形外部的最大边界框之外的所有点)和明显在多边形内部的点(即在多边形内部的最大矩形之内的所有点),但这可能仍然会留下太多的点。 - cdalitz如@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))
ptinpoly::pip2d
对您的需求来说太慢了,我想知道为什么 DescTools::PtInPoly
更快。在我的测试中,情况恰恰相反。 - cdalitz看起来你可以通过将每个点连接到其最近的邻居和相反方向上的最近邻居来重构黑色多边形的边缘。然后进行点在多边形内测试。