Error when attempting to plot runeems_snps results

[mcrossley3]

Hi Dr. Petkova, I came across an error trying to plot my EEMS results in R that I'm not sure how to resolve. Here is some of my last output: Processing the following EEMS output directories : ./run1 Plotting effective migration surface (posterior mean of m rates) ./run1 Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, : TopologyException: Input geom 1 is invalid: Self-intersection at or near point -90.404795806411528 47.0755877112326 at -90.404795806411528 47.0755877112326

I looked into this a bit online, and some have suggested the issue arises when dealing with unprojected coordinates. But I'm not sure how to specify a projection in the data I'm plotting. I tried your code for specifying a projection (page 15 in the EEMS manual), but got the same error. Do you have any insight on how to resolve this error?

-Mike

[dipetkov]

In this case the complex geometry of the habitat (specified in datapath.outer) causes the self-intersection error. For illustration, let's show the habitat outline and the demes in the population grid:

outer = read.table(paste0(mcmcpath,"/outer.txt"))
demes = read.table(paste0(mcmcpath,"/demes.txt"))
plot(outer,type="l")
points(demes)

The self-intersection (or multiple self-intersections) probably occur in the cluster of little islands in the north or the island in the east, where the habitat vertices are spaced very close to each other.

I don't think the geometry can be simplified automatically. Although it will be a fiddly work to do it by hand, this seems most straightforward solution: just remove the cluster of little islands in the north, and I think that there is an island in the east as well. You can remove a whole bunch of points from datapath.outer to simplify the geometry. (The fiddly part is figuring out which sequence of points correspond to the islands.)

Then runeems will run a little faster as well. The MCMC needs to keep the centers of the Voronoi tiles within the habitat. But if the habitat geometry is very complex, it is expensive to check this condition as new Voronoi tiles are added or the centres of existing tiles are updated.

After running runeems to estimate the model parameters, it is easy to plot the complex geographic boundaries on top of the coloured contour plot which represents the effective migration rates. Here is how to do it:

outer = read.table(paste0(datapath,"/outer.txt"))
eems.plots(mcmcpath = mcmcpath,
           plotpath = plotpath,
           longlat = TRUE,
           add.outline = TRUE, col.outline = "black",
           add.demes = TRUE, col.demes = "black",
           add.grid = TRUE, col.grid = "black",
           m.plot.xy = { lines(outer) },    
           q.plot.xy = { lines(outer) })

[dipetkov]

I have modified rEEMSplots to print the following slightly more informative message in cases there is a self-intersection error:

## Error (rgeos):
##  The habitat geometry is not a valid ring (a ring is both simple and closed)
##  Let the habitat be the rectangle defined by (xmin,ymin) and (xmax,yxmax)
##  where xmin,xmax = range(longitude) and ymin,ymax = range(latitude).

As the message explains, if there is a self-intersection error, the plotting script will discard the original habitat geometry and replace it with a rectangle (the bounding box which contains the original habitat).

[mcrossley3]

I made a simpler .outer file by picking 5 vertices that still fully contain all the sample sites, and things are running smoothly now. Thanks!