有没有更好的方法来处理跨越反子午线（日期变更线）的 SpatialPolygons？
Is there a better way for handling SpatialPolygons that cross the antimeridian (date line)?
TL;DR
R 中处理 SpatialPolygons intersecting/overlapping 纬度 +/-180° 处的反子午线并沿该子午线将它们分成两部分的最佳方法是什么？
前言
这将是一篇很长的文章，但这只是因为我将包含大量代码和图表以供说明。我将向您展示我的目标是什么以及我通常如何实现该目标，然后演示如何在字面上的边缘案例中将其分解在一起。正如标题所暗示的那样，我已经找到了一种可能的解决方案来解决我的问题，因此我也将其包括在内。但它不是 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")
这正是我需要的：
问题
现在观察当我们引入靠近缓冲区必须穿过那条线的 anti-meridian（+/-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() 命令确实设法在反子午线的另一侧创建多边形线段（如果 dissolve=FALSE）：
但是多边形穿过整个地球而不是正确环绕（与 0° 相交而不是 180°）。这导致 self-intersections 和
buffr <- gUnaryUnion(buffr@polygons)
将失败
Error in gUnaryUnion(buffr@polygons) : TopologyException: Input
geom 0 is invalid: Self-intersection at or near point
170.08604674698876 12.562175561621103 at 170.08604674698876 12.562175561621103
快速但有点脏的解决方案
首先，我们需要检测多边形是否穿过反子午线。然而，其中 none 实际上相交了 +/-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)
检测并分离 "bad" 多边形后，我通过查看纵坐标将它们分成两个独立的部分。每个具有负值的坐标对进入新的西部多边形，正值进入东部多边形。然后我将它们全部合并回去，做我的 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)
最终结果：
实际问题
所以，这个解决方案适用于我当前的用例，但它有一些问题：
- 一旦其中一个缓冲区同时越过反面和prime-meridian，它可能会完全破裂（如果原始点位置靠近两极，则不太可能）。
- 这不是很准确，因为两个多边形部分不是在 +/-180° 处切割，而是在原始多边形中存在的最高 negative/positive 纬度值处切割。
- 我很难相信没有 "proper" 方法可以做到这一点。
所以所有这些问题归结为：有没有更好的方法？
当我试图解决这个问题时，我发现了 maptools 包中的 nowrapRecenter() and nowrapSpatialPolygons() 函数，乍一看它们完全符合我的要求。仔细检查后，它们的目标几乎是相反的用例（将地图以反子午线为中心，从而沿本初子午线切割多边形）。我和他们一起玩，但没能让他们为我工作——事实上，他们只会让事情变得更糟。
感谢关注！
你是对的，现在是今年，你的问题有解决方案。 sf-package 具有 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 将跨越国际日期变更线或 "antimeridian" 的多边形转换为 MULTIPOLYGON。仅此而已。
这能解决您的问题吗？至少它大大缩短了到达你现在所在位置的方式。 ^^
TL;DR
R 中处理 SpatialPolygons intersecting/overlapping 纬度 +/-180° 处的反子午线并沿该子午线将它们分成两部分的最佳方法是什么？
前言
这将是一篇很长的文章，但这只是因为我将包含大量代码和图表以供说明。我将向您展示我的目标是什么以及我通常如何实现该目标，然后演示如何在字面上的边缘案例中将其分解在一起。正如标题所暗示的那样，我已经找到了一种可能的解决方案来解决我的问题，因此我也将其包括在内。但它不是 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")
看起来像这样：
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")
这正是我需要的：
问题
现在观察当我们引入靠近缓冲区必须穿过那条线的 anti-meridian（+/-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() 命令确实设法在反子午线的另一侧创建多边形线段（如果 dissolve=FALSE）：
但是多边形穿过整个地球而不是正确环绕（与 0° 相交而不是 180°）。这导致 self-intersections 和
buffr <- gUnaryUnion(buffr@polygons)
将失败
Error in gUnaryUnion(buffr@polygons) : TopologyException: Input geom 0 is invalid: Self-intersection at or near point 170.08604674698876 12.562175561621103 at 170.08604674698876 12.562175561621103
快速但有点脏的解决方案
首先，我们需要检测多边形是否穿过反子午线。然而，其中 none 实际上相交了 +/-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)
检测并分离 "bad" 多边形后，我通过查看纵坐标将它们分成两个独立的部分。每个具有负值的坐标对进入新的西部多边形，正值进入东部多边形。然后我将它们全部合并回去，做我的 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)
最终结果：
实际问题
所以，这个解决方案适用于我当前的用例，但它有一些问题：
- 一旦其中一个缓冲区同时越过反面和prime-meridian，它可能会完全破裂（如果原始点位置靠近两极，则不太可能）。
- 这不是很准确，因为两个多边形部分不是在 +/-180° 处切割，而是在原始多边形中存在的最高 negative/positive 纬度值处切割。
- 我很难相信没有 "proper" 方法可以做到这一点。
所以所有这些问题归结为：有没有更好的方法？
当我试图解决这个问题时，我发现了 maptools 包中的 nowrapRecenter() and nowrapSpatialPolygons() 函数，乍一看它们完全符合我的要求。仔细检查后，它们的目标几乎是相反的用例（将地图以反子午线为中心，从而沿本初子午线切割多边形）。我和他们一起玩，但没能让他们为我工作——事实上，他们只会让事情变得更糟。
感谢关注！
你是对的，现在是今年，你的问题有解决方案。 sf-package 具有 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 将跨越国际日期变更线或 "antimeridian" 的多边形转换为 MULTIPOLYGON。仅此而已。
这能解决您的问题吗？至少它大大缩短了到达你现在所在位置的方式。 ^^