用于凸包敏感性分析的空间多边形重叠百分比
Percentage overlap of spatial polygons for a sensitivity analysis of convex hull
为了可重复性，让我们将我的问题简化如下：我有 100 个空间多边形，代表从总体（100 次）中抽取的 N 个随机样本的凸包，以计算模型对单个值的敏感性。 如何计算这些多边形的重叠百分比？理想的解决方案应该是快速的并且引入尽可能少的近似值。
我没有特别的理由使用 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, 1e3, 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". Multiple solutions (1, 2) 已被建议解决此 GEOS 相关问题，none 我已设法使用我的数据：
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, add = TRUE, 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),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
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, 1e3, 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]], add = TRUE, border = "red")
}
# Intersection of all convex hulls plotted in transparent blue
interior <- do.call(intersect.owin, polys)
plot(interior, add = TRUE, col = rgb(0,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)。
- 应该可以注意到 points/edges/triangles 并计算出每个包在哪些壳内（特别是每个边缘的左侧和右侧在哪个壳内，然后可以用来确定每个边缘内有哪些壳三角形（设置线内侧的船体的并集），并从中计算三角形所在的船体数量。棘手的一点是如何在不采用
O(T*H) 的情况下级联信息。想多了再回复。
用更好的方法编辑
Should their intersection be added to the list of points to be triangulated?
Reducing the ambiguity. That technique is a linescan algorithm especially that for detecting intersections in O(Nlog(N)) time, such as the https://en.wikipedia.org/wiki/Bentley%E2%80%93Ottmann_algorithm
所以这里有一个更直接的更新方法包括下面的示例图像（看起来比预期的要小...）
上图显示了 3 个凸包，并为从左到右穿过每个点的扫掠线编号。尽管 Andrew's Algorithm for convex hulls 确实避免了对实际扫描线的需要，因为它是算法的一部分。基本上，您使用安德鲁的算法一次性构建所有船体，但有重复。
所以基本过程如下所示：
- 为每个已知船体（G/R/B：绿色、红色、黑色）、上下船体设置空列表。所以每个点到它们所在的船体的映射（初始化为空列表）。
- 使用安德鲁算法的排序顺序对所有点（在凸包内）进行排序。
- 使用与 Andrew 算法相同的排序顺序，将每个点添加到每个船体（上和下）。
- 我们再用安德鲁的算法来考虑点数。不过，诀窍在于我们已经知道船体将是什么。考虑红色船体，点 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 个船体内。
我会首先在多个船体中没有点的情况下进行测试。然后调整逻辑支持多壳。
为了获得最佳性能，我只会 运行 这个算法在实际船体点上，然后将结果（color-coded 多边形，如果你走线段路线）叠加在实际之上数据集，如果你需要的话。如果你没有采用颜色编码的多边形路线，那么我可能会根据它们所在的平均船体数量或者可能 运行 使用所有点（不仅仅是船体点）的算法来为多边形着色，但就是这样成为巨大的性能打击。只为线段做工作可能更好。
为了可重复性，让我们将我的问题简化如下：我有 100 个空间多边形，代表从总体（100 次）中抽取的 N 个随机样本的凸包，以计算模型对单个值的敏感性。 如何计算这些多边形的重叠百分比？理想的解决方案应该是快速的并且引入尽可能少的近似值。
我没有特别的理由使用 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, 1e3, 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". Multiple solutions (1, 2) 已被建议解决此 GEOS 相关问题，none 我已设法使用我的数据：
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, add = TRUE, 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),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
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, 1e3, 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]], add = TRUE, border = "red")
}
# Intersection of all convex hulls plotted in transparent blue
interior <- do.call(intersect.owin, polys)
plot(interior, add = TRUE, col = rgb(0,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)。 - 应该可以注意到 points/edges/triangles 并计算出每个包在哪些壳内（特别是每个边缘的左侧和右侧在哪个壳内，然后可以用来确定每个边缘内有哪些壳三角形（设置线内侧的船体的并集），并从中计算三角形所在的船体数量。棘手的一点是如何在不采用
O(T*H)的情况下级联信息。想多了再回复。
用更好的方法编辑
Should their intersection be added to the list of points to be triangulated? Reducing the ambiguity. That technique is a linescan algorithm especially that for detecting intersections in
O(Nlog(N))time, such as the https://en.wikipedia.org/wiki/Bentley%E2%80%93Ottmann_algorithm
所以这里有一个更直接的更新方法包括下面的示例图像（看起来比预期的要小...）
上图显示了 3 个凸包，并为从左到右穿过每个点的扫掠线编号。尽管 Andrew's Algorithm for convex hulls 确实避免了对实际扫描线的需要，因为它是算法的一部分。基本上，您使用安德鲁的算法一次性构建所有船体，但有重复。
所以基本过程如下所示：
- 为每个已知船体（G/R/B：绿色、红色、黑色）、上下船体设置空列表。所以每个点到它们所在的船体的映射（初始化为空列表）。
- 使用安德鲁算法的排序顺序对所有点（在凸包内）进行排序。
- 使用与 Andrew 算法相同的排序顺序，将每个点添加到每个船体（上和下）。
- 我们再用安德鲁的算法来考虑点数。不过，诀窍在于我们已经知道船体将是什么。考虑红色船体，点 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 个船体内。
我会首先在多个船体中没有点的情况下进行测试。然后调整逻辑支持多壳。
为了获得最佳性能，我只会 运行 这个算法在实际船体点上，然后将结果（color-coded 多边形，如果你走线段路线）叠加在实际之上数据集，如果你需要的话。如果你没有采用颜色编码的多边形路线，那么我可能会根据它们所在的平均船体数量或者可能 运行 使用所有点（不仅仅是船体点）的算法来为多边形着色，但就是这样成为巨大的性能打击。只为线段做工作可能更好。