R中处理在纬度+/-180°处与反子午线相交/重叠的空间多边形并将其沿该子午线切割成两部分的最佳方法是什么？
这将是一篇很长的文章，但只是因为我将包含大量代码和图形来进行说明。我将向您展示我的目标是什么以及我通常如何实现该目标，然后演示这一切如何在字面上的边缘情况下分解在一起。正如标题所示，我已经找到了解决我的问题的一种可能的解决方案，因此我也将其包括在内。但它并不是 100% 干净，我想看看是否有人能想出更优雅的东西。无论如何，我认为这是一个有趣的问题，因为就在几天前，我在最疯狂的梦想中都不会怀疑这甚至可能成为 2019 年的一个问题。
首先，创建一个有效的示例数据集
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")
Looks like this:
Then, I use circles() from the dismo package to create circular buffers around those locations. I use this function, because it takes into account that the Earth is not flat:
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")
That looks like this:
然后，将单个缓冲区合并为一个大（多）多边形：
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")
This is exactly what I need:
现在观察当我们引入非常接近反子午线（经度+/-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：输入 geom 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)
The final outcome:
因此，这个解决方案适用于我当前的用例，但它有一些问题：
所以这一切归结起来的问题是：有更好的方法吗?
当我试图弄清楚这一点时，我遇到了nowrapRecenter() and nowrapSpatialPolygons() https://rdrr.io/cran/maptools/man/nowrapRecenter.html函数从maptools包装，乍一看似乎完全符合我的要求。经过仔细检查，它们的目标几乎是相反的用例（将地图以反子午线为中心，从而沿着本初子午线切割多边形）。我和他们一起玩，但没能让他们为我工作——事实上，他们只会让事情变得更糟。
感谢您的关注！
你是对的，今年是你的问题的解决方案。这sf- 包有功能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。就是这样。
这能解决你的问题吗？至少它缩短了你到达现在所在位置的路程。 ^^