TL;博士
R中处理在纬度+/- 180°处与反子午线相交/重叠的SpatialPolygons并将它们沿该子午线切成两部分的最佳方法是什么？
前言
这将是一篇很长的文章，但这只是因为我将包含大量代码和图表以进行说明。我将向您展示我的目标是什么以及我通常如何实现这一目标，然后在一个字面的边缘案例中演示它们是如何组合在一起的。正如标题所示，我已经找到了一种可能的解决方案来解决我的问题，因此我也会将其包括在内。但它不是 100% 干净的，我想看看是否有人能想出更优雅的东西。无论如何，我认为这是一个有趣的问题，因为就在几天前，我做梦也不会怀疑这甚至可能在 2019 年成为一个问题。
R中的常规工作流程
首先，创建一个有效的示例数据集
library(sp)
library(rgdal)
library(rgeos)
library(dismo)
library(maptools) # this is just for plotting a simple world map in the background
data("wrld_simpl")
# create a set of locations
locations <- SpatialPoints(coords=cbind(c(50,0,0,0), c(10, 30, 50, 70)), proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
plot(wrld_simpl, border="grey50")
points(locations, pch=19, col="blue")
看起来像这样：
然后，我使用 dismo 包中的 circles() 在这些位置周围创建循环缓冲区。我使用这个函数，因为它考虑到地球不是平的：
buffr <- circles(p = locations, d = 1500000, lonlat=TRUE, dissolve=FALSE)
plot(wrld_simpl, border="grey50")
plot(buffr, add=TRUE, border="red", lwd=2)
points(locations, pch=19, col="blue")
看起来像这样：
然后，将单个缓冲区合并为一个大（多）多边形：
buffr <- buffr@polygons # extract the SpatialPolygons object from the "CirclesRange" object
buffr <- gUnaryUnion(buffr) # merge
plot(wrld_simpl, border="grey50")
plot(buffr, add=TRUE, border="red", lwd=2)
points(locations, pch=19, col="blue")
这正是我需要的：
问题
现在观察当我们引入非常接近缓冲区必须越过该线的反子午线（经度+/-180°）的位置时会发生什么：
locations <- SpatialPoints(coords=cbind(c(50,0,0,0, 175, -170), c(10, 30, 50, 70,0,-10)), proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
buffr <- circles(p = locations, d = 1500000, lonlat=TRUE, dissolve=FALSE)
plot(wrld_simpl, border="grey50")
plot(buffr, add=TRUE, border="red", lwd=2)
points(locations, pch=19, col="blue")
circles() 命令确实设法在反子午线的另一侧创建多边形段（如果溶解 = FALSE）：
但多边形穿过整个地球而不是正确环绕（与 0° 相交而不是 180°）。这会导致自相交和
buffr <- gUnaryUnion(buffr@polygons)
将失败
gUnaryUnion(buffr@polygons) 中的错误：TopologyException：输入几何 0 无效：在点 170.08604674698876 12.562175561621103 处或附近自相交 170.08604674698876 12.562175561621103
快速且略显肮脏的解决方案
首先，我们需要检测一个多边形是否穿过反子午线。但是，它们实际上都没有相交 +/-180°。相反，我使用了两条伪经线，它们靠近真实经线，但距离东西方足够远，可能与所讨论的多边形相交。如果一个多边形与它们两者相交，它也必须穿过反子午线。
antimeridian <- SpatialLines(list(Lines(slinelist=list(Line(coords=cbind(c(179,179), c(90,-90)))), ID="1"),
Lines(slinelist=list(Line(coords=cbind(c(-179,-179), c(90,-90)))), ID="2")),
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
intrscts <- gIntersects(antimeridian, buffr, byid = TRUE)
any(intrscts[,1] & intrscts[,2])
intrscts <- which(intrscts[,1] & intrscts[,2])
buffr.bad <- buffr[intrscts,]
buffr.good <- buffr[-intrscts,]
plot(wrld_simpl)
plot(buffr.good, border="blue", add=TRUE)
plot(buffr.bad, border="red", add=TRUE)
在检测并分离出“坏”多边形后，我只需通过查看纵向坐标将它们分成两个单独的部分。每个具有负值的坐标对进入新的西部多边形，正值进入东部多边形。然后我只是将它们重新合并在一起，做我的 gUnaryUnion 并拥有我需要的东西：
buffr.fixed <- buffr.good
for(i in 1:length(buffr.bad)){
thispoly <- buffr.bad[i,] # select first problematic polygon
crds <- thispoly@polygons[[1]]@Polygons[[1]]@coords # extract coordinates
crds.west <- subset(crds, crds[,1] < 0) # western half of the polygon
crds.east<- subset(crds, crds[,1] > 0)
# turn into Spatial*, merge back together, re-add original crs
sppol.east <- SpatialPolygons(list(Polygons(list(Polygon(crds.east)), paste0("east_", i))))
sppol.west <- SpatialPolygons(list(Polygons(list(Polygon(crds.west)), paste0("west_", i))))
sppol <- spRbind(sppol.east, sppol.west)
proj4string(sppol) <- proj4string(thispoly)
buffr.fixed <- spRbind(buffr.fixed, sppol)
}
buffr.final <- gUnaryUnion(buffr.fixed)
plot(wrld_simpl, border="grey50")
points(locations, pch=19, col="blue")
plot(buffr.final, add=TRUE, border="red", lwd=2)
最终结果：
实际问题
所以，这个解决方案适用于我当前的用例，但它有一些问题：
- 一旦其中一个缓冲区穿过反子午线和本初子午线，它可能会完全中断（如果原始点位置靠近两极，这不太可能）。
- 这不是很精确，因为两个多边形部分不是在 +/-180° 处切割，而是在原始多边形中存在的最高负/正纬度值处切割。
- 我很难相信没有“正确”的方式来做到这一点。
所以这一切归结为一个问题：有没有更好的方法来做到这一点？
当我试图弄清楚这一点时，我发现了包中的nowrapRecenter()andnowrapSpatialPolygons()函数maptools，乍一看，它们看起来就像我想要的那样。经过仔细检查，它们几乎针对相反的用例（以反子午线为中心，从而沿本初子午线切割多边形）。我和他们一起玩，但没有让他们为我工作——事实上，他们只会让事情变得更糟。
感谢您的关注！
