dipetkov / eems

Estimating Effective Migration Surfaces
GNU General Public License v2.0
102 stars 28 forks source link

Unable to add points in rEEMSplots #38

Open nehasavant opened 5 years ago

nehasavant commented 5 years ago

Hello,

I'm using rEEMSplots to plot results, and I'd like to add the original points to the migration surface graph. I use the following syntax and do not get any errors, but the points are not appearing on my output graph.

plot_p <- rEEMSplots::eems.plots(mcmcpath = eems_results_P, plotpath = names_figures_P, add.grid = TRUE, projection.in = "+proj=longlat +datum=WGS84", projection.out = "+proj=longlat +datum=WGS84", add.outline = TRUE, longlat = FALSE, m.plot.xy = { points(coord_P, col = "purple", pch=18) }, q.plot.xy = { points(coord_P, col = "purple", pch=18) })

In addition, is it possible to overlay a streams/waterbodies shapefile on top of the migration/diversity plots?

Thank you!

tim-oconnor commented 4 years ago

Hi nehasavant -- late reply, but perhaps it'll help you / others. This is absolutely possible, but I've also had issues before. You can just use the "map" or "plot" commands that you want within the m.plot.xy and q.plot.xy arguments. The only trick is to ensure that you first project the points or shapefiles into the same CRS as the EEMS rasters (use spTransform from the sp / rgdal packages).

-Tim

dipetkov commented 4 years ago

Anything inside the m.plot.xy and m.plot.xy is added on top of the raster plot but the statements are not parsed (in order to apply projections, for example). So you have to do it by hand.

Here is the relevant example:

library("rEEMSplots")
extdata_path <- system.file("extdata", package = "rEEMSplots")
eems_results <- file.path(extdata_path, "EEMS-example")
name_figures <- file.path(path.expand("~"), "EEMS-example")

library("rgdal")
projection_none <- "+proj=longlat +datum=WGS84"
projection_mercator <- "+proj=merc +datum=WGS84"

# Start with points specified in (longitude, latitude)
coords <- matrix(c(
  -10, 10,
  10, 10,
  30, 0,
  40, -10,
  30, -20
), ncol = 2, byrow = TRUE)
colors <- c("red", "green", "blue", "purple", "orange")
labels <- LETTERS[1:5]

# Transform the coordinates into the same projection
# that will be used to transform the migration rates
coords_merc <- sp::spTransform(
  SpatialPoints(coords, CRS(projection_none)),
  CRS(projection_mercator)
)

# `coords_merc` is a SpatialPoints structure
# but we only need the coordinates themselves
coords_merc <- coords_merc@coords

# So to summarize start with points in one projection
coords
# And project them to the target projection
coords_merc

eems.plots(
  mcmcpath = eems_results,
  plotpath = paste0(name_figures, "-labels-projected"),
  longlat = TRUE,
  projection.in = projection_none,
  projection.out = projection_mercator,
  m.plot.xy = {
    text(coords_merc, col = colors, pch = labels, font = 2)
  },
  q.plot.xy = {
    text(coords_merc, col = colors, pch = labels, font = 2)
  }
)