TL; DR
Как лучше всего использовать R для обработки пространственных полигонов, пересекающих / перекрывающих анти-меридиан на широте +/- 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")
Похоже на это: 
Затем я использую circle () из пакета dismo для созданиякольцевые буферы вокруг этих мест.Я использую эту функцию, потому что она учитывает, что Земля не плоская:
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")
Команда circle () делаетудается создать сегменты многоугольника на другой стороне антимеридиана (если распустить = ЛОЖЬ):
, но многоугольник пересекает весь земной шар вместо переносавокруг (пересекается с 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() и nowrapSpatialPolygons() из пакета maptools, которые на первый взгляд выглядели так, как будто они делают именно то, что я хочу. При ближайшем рассмотрении они нацелены на почти противоположный вариант использования (центрирование карты на анти-меридиане и, таким образом, разрезание полигонов вдоль основного меридиана). Я играл с ними, но не смог заставить их работать на меня - на самом деле, они только усугубили ситуацию.
Спасибо за внимание!