空间多边形重叠百分比，用于凸包的灵敏度分析

如何解决空间多边形重叠百分比，用于凸包的灵敏度分析

为了重现性，让我们简化我的问题，如下所示：我有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船壳内。

我将首先在没有多个壳体的情况下对此进行测试。然后调整逻辑以支持多个船体。

为了获得最佳性能，我只会在实际的船体点上运行此算法，然后根据需要将结果（如果使用线段路线，则将颜色编码的多边形）覆盖在实际数据集的顶部。如果您没有采用颜色编码的多边形路线，那么我可能会根据它们所在的平均船体数量为多边形着色，或者使用所有点（不仅是船体点）运行算法，但这将是一个巨大的过程性能受到打击。仅对线段进行处理可能更好。

空间多边形重叠百分比，用于凸包的灵敏度分析

如何解决空间多边形重叠百分比，用于凸包的灵敏度分析

解决方法

采用更好的方法进行编辑

相关推荐