TL; DR

¿Cuál es la mejor manera en R para manejar los Polígonos Espaciales que se cruzan / superponen al anti meridiano a +/- 180 ° de latitud y los corta en dos secciones a lo largo de ese meridiano?

Prefaci

Esto será largo, pero solo porque voy a incluir una gran cantidad de código y figuras para ilustrar. Te mostraré cuál es mi objetivo y cómo lo logro normalmente y luego demostraré cómo se rompe todo en un caso literal. Como sugiere el título, ya encontré una posible solución a mi problema, así que también la incluiré. Pero no está 100% limpio y me gustaría ver si alguien puede llegar a algo más elegante. En cualquier caso, creo que este es un problema interesante, ya que hace solo un par de días no habría sospechado en mis sueños más locos que esto podría ser un problema en 2019.

Flujo de trabajo regular en R

Primero, cree un conjunto de datos de ejemplo que funcione

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")

Se ve como esto:Luego, uso los círculos () del paquete de dismo para crear búferes circulares alrededor de esas ubicaciones. Utilizo esta función, porque tiene en cuenta que la Tierra no es plana:

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")

Eso se ve así:

uego, combine las memorias intermedias individuales en un gran polígono (multi):

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")

Esto es exactamente lo que necesito:

El problem

Ahora observe lo que sucede cuando introducimos ubicaciones que están tan cerca del anti-meridiano (+/- 180 ° de longitud) que el búfer tiene que cruzar esa línea:

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")

El comando circles () logra crear segmentos de polígono del otro lado del antimeridiano (si se disuelve = FALSO):

pero el polígono cruza todo el globo en lugar de envolverse adecuadamente (intersectando con 0 ° en lugar de 180 °). Eso lleva a auto-intersecciones y

buffr <- gUnaryUnion(buffr@polygons)

fallará con

Error en gUnaryUnion (buffr @ polygons): TopologyException: la entrada geom 0 no es válida: Auto-intersección en o cerca del punto 170.08604674698876 12.562175561621103 en 170.08604674698876 12.562175561621103

La solución rápida y ligeramente sucia

Primero, necesitamos detectar si un polígono cruza el anti meridiano. Sin embargo, ninguno de ellos se cruza realmente +/- 180 °. En cambio, estoy usando dos pseudo anti meridianos que se encuentran cerca del real, pero lo suficientemente al este y al oeste como para probablemente intersecar los polígonos en cuestión. Si un polígono se cruza con ambos, también debe cruzar el anti meridiano.

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)

Después de detectar y separar los polígonos "malos", simplemente los dividí en dos secciones separadas mirando las coordenadas longitudinales. Cada par de coordenadas que tiene un valor negativo allí va al nuevo polígono occidental, los positivos al oriental. Luego, combino todo de nuevo, hago mi gunaryUnion y tengo casi lo que necesito:

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)

El resultado final:

La pregunta real

Entonces, esta solución funciona para mí para mi caso de uso actual, pero tiene algunos problemas:

s probable que se rompa completamente tan pronto como uno de los amortiguadores cruce tanto el meridiano anti como el meridiano primario (lo cual no es tan improbable si las ubicaciones de los puntos originales se encuentran cerca de los polos).it no es del todo exacto, ya que las dos secciones de polígono no se cortan a +/- 180 ° sino a los valores negativos / positivos más altos de latitud que estaban presentes en el polígono original.e resulta difícil creer que no haya una forma "adecuada" de hacerlo.

Así que la pregunta a la que todo se reduce es: Hay una mejor manera de hacer esto?

Mientras intentaba resolver esto, me encontré con lanowrapRecenter() ynowrapSpatialPolygons() funciones demaptools package, que a primera vista parecía que hacían exactamente lo que yo quería. Tras una inspección más cercana, apuntan al caso de uso opuesto (centrando un mapa en el anti meridiano y, por lo tanto, cortando polígonos a lo largo del meridiano principal). Jugué con ellos, pero no pude hacer que funcionaran para mí; de hecho, solo lograron empeorar las cosas.

¡Gracias por su atencion