rpolygongeosr-maptoolsdismo

Is there a better way for handling SpatialPolygons that cross the antimeridian (date line)?

TL;DR

What is the best way in R to handle SpatialPolygons intersecting/overlapping the anti meridian at +/-180° of latitude and cut them into two sections along that meridian?

Preface

This is going to be a long one, but only because I'm going to include a lot of code and figures for illustration. I'll show you what my goal is and how I normally achieve that and then demonstrate how it all breaks together in a literal edge case. As the title suggests, I already found one possible solution to my problem, so I'll include that too. But it is not 100% clean and I'd like to see if somebody can come up with something more elegant. In any case I think this is an interesting problem, as only a couple of days ago I wouldn't have have suspected in my wildest dreams that this could even be an issue in 2019.

Regular work flow in R

First, create an example data set that works

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: locations on a world map 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: locations with their individual buffers

Then, merge the single buffers into one big (multi-) polygon:

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: points with merged buffers

The problem

Now observe what happens when we introduce locations that so are close to the anti-meridian (+/-180° of longitude) that the buffer has to cross that line:

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

The circles() command does manage to create polygon segments on the other side of the antimeridian (if dissolve=FALSE):

enter image description here

but the polygon crosses the entire globe instead of wrapping around properly (intersecting with 0° instead of 180°). That leads to self-intersections and 

buffr <- gUnaryUnion(buffr@polygons)

will fail with

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

The quick and slightly dirty solution

First, we need to detect whether a polygon crosses the anti meridian. However, none of them actually intersects +/-180°. Instead, I'm using two pseudo anti meridians that lie close to the real one, but far enough to the east and west to probably intersect the polygons in question. If a polygon intersects both of them, it must also cross the anti meridian.

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)

After having detected and separated the "bad" polygons I simply split them into two separate sections by looking at the longitudinal coordinates. Every coordinate pair that has a negative value there goes into the new western polygon, positive ones into the eastern one. Then I just merge it all back together, do my gUnaryUnion and have pretty much what I need:

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: merged buffers, problem solved?

The actual question

So, this solution works for me for my current use case, but it has some issues:

So the question all this boils down to is: Is there a better way of doing this?

While I was trying to figure this out, I came across the nowrapRecenter() and nowrapSpatialPolygons() functions from the maptools package, which at first sight looked like they do exactly what I want. Upon closer inspection, they are aimed at pretty much the opposite use case (centering a map on the anti meridian and thus cutting polygons along the prime meridian). I played around with them, but failed to make them work for me – in fact, they only managed to make things worse.

Thanks for you attention!

Solution

  • You're right, it's the current year and there is a solution for your problem. The sf-package has the function st_wrap_dateline(), which is exactly what you need.

    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 converts the polygons which cross the international date line, or "antimeridian", into MULTIPOLYGON. And that's about it.

    Does that solve your problem? At least it shortens the way quite a bit, to get where you are now. ^^

    enter image description here