R

Как лучше всего использовать 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 () удается создать сегменты многоугольника на другой стороне антимеридиана (если disolve = 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() а такжеnowrapSpatialPolygons() функции отmaptools пакет, который с первого взгляда выглядел так, как будто они делают именно то, что я хочу. При ближайшем рассмотрении они нацелены на почти противоположный вариант использования (центрирование карты на анти-меридиане и, таким образом, разрезание полигонов вдоль основного меридиана). Я играл с ними, но не смог заставить их работать на меня - на самом деле, они только усугубили ситуацию.

Спасибо за внимание!