问题描述
为了重现性，让我们简化我的问题，如下所示：我有100个空间多边形，这些多边形代表从总体中抽取的N个随机样本的凸包（100次），以计算模型对单个值的敏感性。 如何计算这些多边形的重叠百分比？理想的解决方案应该是快速的并引入尽可能小的近似值。
除了我认为这可能是解决问题的最简单方法之外，我没有其他理由使用R的GIS功能。
library(sp)
library(raster)
library(sf)
#> Linking to GEOS 3.8.1,GDAL 3.1.1,PROJ 6.3.1
set.seed(11)
dt <- data.frame(x = rnorm(1e3,10,3) + sample(-5:5,1e3,replace = TRUE))
dt$y <- (rnorm(1e3,3,4) + sample(-10:10,replace = TRUE)) + dt$x
dt <- rbind(dt,data.frame(x = -dt$x,y = dt$y))
plot(dt,asp = 1)
dt.chull <- dt[chull(dt),]
dt.chull <- rbind(dt.chull,dt.chull[1,])
lines(dt.chull,col = "green")
uncert.polys <- lapply(1:100,function(i) {
tmp <- dt[sample(rownames(dt),1e2),]
# points(tmp,col = "red")
tmp <- tmp[chull(tmp),]
tmp <- rbind(tmp,tmp[1,])
tmp <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(tmp)),ID = i)))
sp::SpatialPolygonsDataFrame(tmp,data = data.frame(id = i,row.names = i))
# lines(tmp,col = "red")
})
polys <- do.call(rbind,uncert.polys)
plot(polys,add = TRUE,border = "red")
我最初的尝试是使用sf::st_intersection函数：
sf.polys <- sf::st_make_valid(sf::st_as_sf(polys))
all(sf::st_is_valid(sf.polys))
#> [1] TRUE
sf::st_intersection(sf.polys)
#> Error in CPL_nary_intersection(x): Evaluation error: TopologyException: found non-noded intersection between LINESTRING (-9.80706 -0.619557,-7.66331 -3.55177) and LINESTRING (-9.80706 -0.619557,-9.80706 -0.619557) at -9.8070645468969637 -0.61955676978603658.
该错误很可能与多边形线"that are almost coincident but not identical"有关。建议使用多种解决方案（1，2）来解决与GEOS相关的问题，但我都没有设法处理数据：
sf.polys <- sf::st_set_precision(sf.polys,1e6)
sf.polys <- sf::st_snap(sf.polys,sf.polys,tolerance = 1e-4)
sf::st_intersection(sf.polys)
#> Error in CPL_nary_intersection(x): Evaluation error: TopologyException: found non-noded intersection between LINESTRING (-13.7114 32.7341,3.29417 30.3736) and LINESTRING (3.29417 30.3736,3.29417 30.3736) at 3.2941702528617176 30.373627946201278.
因此，我必须使用栅格化来近似多边形重叠：
GT <- sp::GridTopology(cellcentre.offset = c(round(min(dt$x),1),round(min(dt$y),1)),cellsize = c(diff(round(range(dt$x),1))/100,diff(round(range(dt$y),1))/100),cells.dim = c(100,100)
)
SG <- sp::SpatialGrid(GT)
tmp <- lapply(seq_along(uncert.polys),function(i) {
out <- sp::over(SG,uncert.polys[[i]])
out[!is.na(out)] <- 1
out[is.na(out)] <- 0
out
})
tmp <- data.frame(overlapping.n = Reduce("+",lapply(tmp,"[[",1)))
tmp$overlapping.pr <- 100*tmp$overlapping.n/100
uncert.data <- SpatialGridDataFrame(SG,tmp)
## Plot
plot(x = range(dt$x),y = range(dt$y),type = "n"
)
plot(raster::raster(uncert.data),col = colorRampPalette(c("white","red","blue","white"))(100),add = TRUE)
plot(polys,border = adjustcolor("black",alpha.f = 0.2),cex = 0.5)
points(dt,pch = ".",col = "black",cex = 3)
lines(dt.chull,col = "green")
该方法可以得出结果，但是输出是近似值，需要花费很长时间进行处理。必须有一种更好的方法。
出于性能比较的目的，这是我当前的解决方案：
gridOverlap <- function(dt,uncert.polys) {
GT <- sp::GridTopology(cellcentre.offset = c(round(min(dt$x),100)
)
SG <- sp::SpatialGrid(GT)
tmp <- lapply(seq_along(uncert.polys),function(i) {
out <- sp::over(SG,uncert.polys[[i]])
out[!is.na(out)] <- 1
out[is.na(out)] <- 0
out
})
tmp <- data.frame(overlapping.n = Reduce("+",1)))
tmp$overlapping.pr <- 100*tmp$overlapping.n/100
SpatialGridDataFrame(SG,tmp)
}
system.time(gridOverlap(dt = dt,uncert.polys = uncert.polys))
# user system elapsed
# 3.011 0.083 3.105
对于较大的数据集，性能至关重要（此解决方案在实际应用中需要花费几分钟的时间。）
由reprex package（v0.3.0）于2020-09-01创建
解决方法
这是使用spatstat查找内部而没有任何错误的解决方案
以及基础的polyclip软件包。
library(spatstat)
# Data from OP
set.seed(11)
dt <- data.frame(x = rnorm(1e3,10,3) + sample(-5:5,1e3,replace = TRUE))
dt$y <- (rnorm(1e3,3,4) + sample(-10:10,replace = TRUE)) + dt$x
dt <- rbind(dt,data.frame(x = -dt$x,y = dt$y))
# Converted to spatstat classes (`ppp` not strictly necessary -- just a habit)
X <- as.ppp(dt,W = owin(c(-25,25),c(-15,40)))
p1 <- owin(poly = dt[rev(chull(dt)),])
# Plot of data and convex hull
plot(X,main = "")
plot(p1,add = TRUE,border = "green")
# Convex hulls of sampled points in spatstat format
polys <- lapply(1:100,function(i) {
tmp <- dt[sample(rownames(dt),1e2),]
owin(poly = tmp[rev(chull(tmp)),])
})
# Plot of convex hulls
for(i in seq_along(polys)){
plot(polys[[i]],border = "red")
}
# Intersection of all convex hulls plotted in transparent blue
interior <- do.call(intersect.owin,polys)
plot(interior,col = rgb(0,1,0.1))
我不清楚您要从这里做什么，但至少是这种方法 避免了多边形裁剪的错误。
要在spatstat中执行基于网格的解决方案，我会将Windows转换为
二进制图像蒙版，然后从那里开始工作：
Wmask <- as.im(Window(X),dimyx = c(200,200))
masks <- lapply(polys,as.im.owin,xy = Wmask,na.replace = 0)
maskmean <- Reduce("+",masks)/100
plot(maskmean)
速度取决于您选择的分辨率，但是我想这可能很多
比使用sp / raster的当前建议更快（这可能可以
使用与此处相同的逻辑进行了很多改进，因此这将是另一种
坚持使用raster）。
编辑在下面进一步改进了一个可能更好的解决方案。
对这个问题稍作思考，我的想法是三角剖分和动态编程方法可以很好地工作。
- 考虑每个凸包的点和线。将它们标记为它们所属的船体（可能存储在查找中）
- 从所有直线上取点并进行三角测量，这些三角形将被记为它们内有多少个凸包。
- 在这一点上，有很多方法可以确定三角形中有多少个凸包。您显示的示例倾向于一些可能的优化，但是作为一般解决方案，最佳方法可能是在每个三角形上循环并查看其中的船体
O(T*H)。 - 应该有可能记下点/边/三角形，并找出每个其中的船体（特别是每个边沿的左侧和右边是哪些船体，然后可以用来确定每个三角形中的船体（设置其中船体位于行的内侧的并集），然后从中获得三角形所在的船体数的计数。棘手的一点是如何不带
O(T*H)来级联信息。并在稍后回复。
采用更好的方法进行编辑
是否应将它们的交点添加到要进行三角测量的点列表中？ 减少歧义。该技术是一种线扫描算法，尤其是用于在
O(Nlog(N))时间内检测交点的算法，例如https://en.wikipedia.org/wiki/Bentley%E2%80%93Ottmann_algorithm
所以这是一个更新的方法，它更直接一些包括下面的示例图像（看起来比预期的要小...）
上图显示了3个凸包，并具有从左到右横穿每个点的扫掠线编号。尽管实际上Andrew's Algorithm for convex hulls避免了实际扫描线的需要，因为其中一个是算法的一部分。基本上，您可以使用安德鲁（Andrew）的算法一次性构建所有船体，但是具有重复性。
因此基本过程如下：
- 为每个已知船体（G / R / B：绿色，红色，黑色），上，下船体设置空列表。因此，每个点到它们所在的外壳的映射（初始化为空列表）。
- 使用安德鲁算法的排序顺序对所有点（在凸包内）进行排序。
- 使用与安德鲁算法相同的排序顺序，将每个点添加到每个船体（上下）。
- 然后，我们使用安德鲁的算法来考虑点。但诀窍在于，我们已经知道船体将是什么。考虑红色的船体，点2，7和8。其他点4和5（5实际上是2点，忘记了标签）。将添加4作为船体点，但是由于我们专注于红色船体，因此我们将忽略4（因为它不在灰色船体之内）。如果多个船体使用同一个点，则同样适用，因为从技术上讲该点不在任何一个船体之内（除非您要考虑这一点，在这种情况下，所有船体点均在至少1个船体之内，所以这样做可能会有用为了视觉上的好处，我认为这是使交点着色实用的唯一方法。但是，这两个5点都在灰色船体之内，因此我们注意到它们都在红色船体之内。总体表现为O（N * C），其中
N是点数，C是船体数。我想这可能可以降为O(C log N + N log C)之类的东西，或者付出足够的努力，但可能不值得。
您可以运行设置的相交找到所有相交，然后使用它们来构建多边形以获得更精确的着色。但是，这使事情变得更加混乱，而且我仍在尝试制定一个好的解决方案。但是，我怀疑，将一个点数计算为“在其自身的船体之内”可能会对此有所帮助。在这种情况下，您可能只需取组成多边形的点的最小值。因此，如果点在1/2/2/2船壳内，则该区域在1船壳内。
我将首先在没有多个壳体的情况下对此进行测试。然后调整逻辑以支持多个船体。
为了获得最佳性能，我只会在实际的船体点上运行此算法，然后根据需要将结果（如果使用线段路线，则将颜色编码的多边形）覆盖在实际数据集的顶部。如果您没有采用颜色编码的多边形路线，那么我可能会根据它们所在的平均船体数量为多边形着色，或者使用所有点（不仅是船体点）运行算法，但这将是一个巨大的过程性能受到打击。仅对线段进行处理可能更好。
设置时间 控制面板
错误1：Request method ‘DELETE‘ not supported 错误还原：...