Modifying maps from EEMS

[jessicarick]

Is it possible to modify the maps in any way-- and how I would go about doing that? For example, I would like to produce a map with the EEMS surface and my sample points on it. I know how to produce such a map in R with a plain map, but the EEMS program just outputs a PDF of the maps, and therefore they aren't modifiable as they are. Would I need to go in and edit the R code for the program to be able to do this? Or is there an easier way?

Thanks!

[dipetkov]

How to add points to the estimated effective migration surface (the contour plot produced by eems.plots)

library(dplyr,warn.conflicts=FALSE) ## Useful package for manipulating data frames
library(rEEMSplots) ## Defines eems.plots, which visualizes EEMS results

I will use results from running EEMS once with nDemes = 540. [But eems.plots can also combine the results from multiple EEMS runs.]

mcmcpath = "./elephants/V3b-nohyb-nIndiv1124-nDemes540-chain1"

There is no way to get the original sampling coordinates from the EEMS output. The genotypes are not copied either. [I didn't want runeems to create copies of the valuable real data.] Instead the coordinates must be read in from the original input file.

datapath = "V3b-nohyb-nIndiv1124"
coord = read.table(paste0(datapath,".coord"))
plot(coord,xlab="longitude",ylab="latitude")

The scatter of points has the shape of Sub-Saharan Africa, approximately. The goal is to add them the contour plot of the estimated migration rates.

eems.plots(mcmcpath = mcmcpath,
           plotpath = mcmcpath,
           longlat = TRUE,
           m.plot.xy = { points(coord,col="purple") },
           q.plot.xy = { points(coord,col="purple") })

Above I have specified two optional arguments which are NULL by default: m.plot.xy adds graphical elements -- points, lines -- to the effective migration sufrace; q.plot.xy adds to the effective diversity surfaces. The two arguments are not required to have the same value (but have to be specified separately). Also note that I have used the points function, not the plot function, because plot creates a new graphical frame. We want to add to an existing one instead.

You can add more graphical elements to the plot, inside the curly brackets:

eems.plots(mcmcpath = mcmcpath,
           plotpath = mcmcpath,
           longlat = TRUE,
           m.plot.xy = { points(coord,col="purple"); 
                         lines( c(-12.1,38.7), c(9.7, -6.0) ); },
           q.plot.xy = { points(coord,col="purple"); })

This adds a single line. The functionality to add a line is not working properly: now the line is straight but it would be better if it curves with the surface of the earth instead.

It is possible to add the geographic map (say, in the Mercator projection) as well as the sampling locations.

projection.in = "+proj=longlat +datum=WGS84"
projection.out = "+proj=merc +datum=WGS84"
eems.plots(mcmcpath = mcmcpath,
           plotpath = mcmcpath,
           longlat = TRUE,
           projection.in = projection.in,
           projection.out = projection.out,
           add.map = TRUE,
           m.plot.xy = { points(coord,col="purple") },
           q.plot.xy = { points(coord,col="purple") })

In this first attempt the sampling locations are obviously not added correctly - note the single purple dot in the Atlantic ocean. The reason is that the colored surface and the map have been transformed to Mercator but the sampling locations have not.

So let's transform the sampling coordinates using the spTransform function in the sp package.

library(sp)
## Coord is a matrix, with two coordinates (longitude and latitude) per row
## Declare each row to be a "Spatial Point"
coord = SpatialPoints(coord,proj4string=CRS(projection.in))
## Transform the spatial points into the Mercator projection
transformed.coord = spTransform(coord,CRSobj=CRS(projection.out))
eems.plots(mcmcpath = mcmcpath,
           plotpath = mcmcpath,
           longlat = TRUE,
           projection.in = projection.in,
           projection.out = projection.out,
           add.map = TRUE,
           m.plot.xy = { points(transformed.coord,col="purple") },
           q.plot.xy = { points(transformed.coord,col="purple") })