是否有更好的方法来处理穿过反子午线（日界线）的空间多边形？

Question

长话短说

\n\n

R中处理在纬度+/-180\xc2\xb0处与反子午线相交/重叠的SpatialPolygons并将它们沿该子午线切割成两部分的最佳方法是什么？

\n\n

前言

\n\n

这将是一篇很长的文章，但只是因为我将包含大量代码和图形来进行说明。我将向您展示我的目标是什么以及我通常如何实现该目标，然后演示这一切如何在字面上的边缘情况下分解在一起。正如标题所示，我已经找到了解决问题的一种可能的解决方案，因此我也将其包括在内。但它并不是 100% 干净，我想看看是否有人能想出更优雅的东西。无论如何，我认为这是一个有趣的问题，因为就在几天前，我在我最疯狂的梦想中都不会怀疑这甚至可能成为 2019 年的一个问题。

\n\n

R 中的常规工作流程

\n\n

首先，创建一个有效的示例数据集

\n\n

library(sp)\nlibrary(rgdal)\nlibrary(rgeos)\nlibrary(dismo)\nlibrary(maptools) # this is just for plotting a simple world map in the background\ndata("wrld_simpl")\n\n\n# create a set of locations\nlocations <- 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"))\nplot(wrld_simpl, border="grey50")\npoints(locations, pch=19, col="blue")\n

\n\n

看起来像这样：\n\n然后，我使用 dismo 包中的 Circles() 在这些位置周围创建循环缓冲区。我使用这个函数，因为它考虑到地球不是平的：

\n\n

buffr <- circles(p = locations, d = 1500000, lonlat=TRUE, dissolve=FALSE)\nplot(wrld_simpl, border="grey50")\nplot(buffr, add=TRUE, border="red", lwd=2)\npoints(locations, pch=19, col="blue")\n

\n\n

看起来像这样：\n

\n\n

然后，将单个缓冲区合并为一个大（多）多边形：

\n\n

buffr <- buffr@polygons # extract the SpatialPolygons object from the "CirclesRange" object\nbuffr <- gUnaryUnion(buffr) # merge\n\nplot(wrld_simpl, border="grey50")\nplot(buffr, add=TRUE, border="red", lwd=2)\npoints(locations, pch=19, col="blue")\n

\n\n

这正是我所需要的：\n

\n\n

问题

\n\n

现在观察当我们引入非常接近反子午线（经度的+/-180\xc2\xb0）的位置时会发生什么，缓冲区必须穿过该线：

\n\n

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"))\nbuffr <- circles(p = locations, d = 1500000, lonlat=TRUE, dissolve=FALSE)\n\nplot(wrld_simpl, border="grey50")\nplot(buffr, add=TRUE, border="red", lwd=2)\npoints(locations, pch=19, col="blue")\n

\n\n

Circles() 命令确实设法在反子午线的另一侧创建多边形线段（如果溶解 = FALSE）：

\n\n

\n\n

但多边形穿过整个地球，而不是正确环绕（与 0\xc2\xb0 相交，而不是 180\xc2\xb0）。这会导致自相交和

\n\n

buffr <- gUnaryUnion(buffr@polygons)\n

\n\n

将会失败

\n\n

\n

gUnaryUnion(buffr@polygons)中的错误：TopologyException：输入\n几何0无效：在点或附近的自相交\n 170.08604674698876 12.562175561621103在170.08604674698876 12.562175561621103

\n

\n\n

快速但有点脏的解决方案

\n\n

首先，我们需要检测多边形是否穿过反子午线。然而，它们实际上都不与 +/-180\xc2\xb0 相交。相反，我使用两条伪反子午线，它们靠近真实子午线，但距离东部和西部足够远，可能与所讨论的多边形相交。如果一个多边形与它们两者相交，它也必须穿过反子午线。

\n\n

antimeridian <- SpatialLines(list(Lines(slinelist=list(Line(coords=cbind(c(179,179), c(90,-90)))), ID="1"),\n              Lines(slinelist=list(Line(coords=cbind(c(-179,-179), c(90,-90)))), ID="2")),\n                                   proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))\nintrscts <- gIntersects(antimeridian, buffr, byid = TRUE)\nany(intrscts[,1] & intrscts[,2])\nintrscts <- which(intrscts[,1] & intrscts[,2])\nbuffr.bad <- buffr[intrscts,]\nbuffr.good <- buffr[-intrscts,]\n\nplot(wrld_simpl)\nplot(buffr.good, border="blue", add=TRUE)\nplot(buffr.bad, border="red", add=TRUE)\n

\n\n

在检测并分离“坏”多边形后，我只需通过查看纵向坐标将它们分成两个单独的部分。每个具有负值的坐标对进入新的西部多边形，正值进入东部多边形。然后我将它们全部合并在一起，执行我的 gUnaryUnion 并得到几乎我需要的东西：

\n\n

buffr.fixed <- buffr.good\nfor(i in 1:length(buffr.bad)){\n  thispoly <- buffr.bad[i,]                              # select first problematic polygon\n  crds <- thispoly@polygons[[1]]@Polygons[[1]]@coords   # extract coordinates\n  crds.west <- subset(crds, crds[,1] < 0)               # western half of the polygon\n  crds.east<- subset(crds, crds[,1] > 0)\n  # turn into Spatial*, merge back together, re-add original crs\n  sppol.east <- SpatialPolygons(list(Polygons(list(Polygon(crds.east)), paste0("east_", i))))\n  sppol.west <- SpatialPolygons(list(Polygons(list(Polygon(crds.west)), paste0("west_", i))))\n  sppol <- spRbind(sppol.east, sppol.west)\n  proj4string(sppol) <- proj4string(thispoly)\n\n  buffr.fixed <- spRbind(buffr.fixed, sppol)\n}\nbuffr.final <- gUnaryUnion(buffr.fixed)\nplot(wrld_simpl, border="grey50")\npoints(locations, pch=19, col="blue")\nplot(buffr.final, add=TRUE, border="red", lwd=2)\n

\n\n

最终结果：\n

\n\n

实际问题

\n\n

因此，这个解决方案适用于我当前的用例，但它有一些问题：

\n\n

\n
• 一旦其中一个缓冲区同时穿过反子午线和本初子午线，它可能就会完全破坏（如果原始点位置靠近两极，则这种情况不太可能发生）。
\n
• 它不太精确，因为两个多边形部分不是在 +/-180\xc2\xb0 处切割的，而是在原始多边形中存在的最高负/正纬度值处切割的。
\n
• 我很难相信没有“正确”的方法可以做到这一点。
\n

\n\n

因此，所有这一切都归结为一个问题：是否有更好的方法来做到这一点？

\n\n

当我试图弄清楚这一点时，我发现了包中的nowrapRecenter()和nowrapSpatialPolygons()函数maptools，乍一看它们似乎完全符合我的要求。经过仔细检查，它们的目标几乎是相反的用例（将地图以反子午线为中心，从而沿着本初子午线切割多边形）。我和他们一起玩过，但没能让他们为我工作 \xe2\x80\x93 事实上，他们只会让事情变得更糟。

\n\n

感谢您的关注！

\n

Answer 1

你是对的，今年是你的问题的解决方案。-packagesf具有 function st_wrap_dateline()，这正是您所需要的。

library(dismo)
library(sf)


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)

buffr2 <- as(buffr@polygons, Class = "sf") %>%
            st_wrap_dateline(options = c("WRAPDATELINE=YES")) %>%
             st_union()


plot(wrld_simpl, border="grey50")
plot(buffr2, add=TRUE, border="red", lwd=2)
points(locations, pch=19, col="blue")

st_wrap_dateline将跨越国际日期变更线或“反子午线”的多边形转换为MULTIPOLYGON. 就是这样。

这能解决你的问题吗？至少它缩短了你到达现在所在位置的路程。^^