-
Notifications
You must be signed in to change notification settings - Fork 28
Description
Hi everyone,
I ran EEMS and generated the migration rate map for my project. However, I noticed that the map shows projections over Indiana and Ohio, even though I don't have any coordinates or sampling points in those areas. I'm a bit surprised by this and was wondering if it's normal for EEMS to project migration rates outside the provided .coord file.
If anyone could clarify this, I'd really appreciate it. Below is the code I used:
###set working directory.
setwd("E:/Armadillo_Analysis/Final_Armadillo_September_3/EEMS_armadillo_alfredo/Testing/")
##install devtools
library("devtools")
install_github("dipetkov/reemsplots2")
library("reemsplots2")
##install the library
library(terra)
library(tidyterra)
library(tidyverse)
library(ggspatial)
states <- vect("cb_2018_us_state_20m.shp")
states
occ <- read.csv("Dasy_original_GBIF+92_samples.csv")
occ <- occ[occ$countryCode == "US", ]
occ <- vect(occ, crs = "EPSG:4326", geom = c("Longitude", "Latitude"))
occ_shp <- terra::convHull(occ)
##occ_shp <- terra::buffer(occ_shp, 2e4)
occ_shp <- terra::buffer(occ_shp, 16093) # 10 miles in meters
occ_shp <- terra::crop(occ_shp, states)
plot(occ_shp)
studyarea_points <- as.points(as.lines(occ_shp))
write.csv(crds(studyarea_points), "Armadillo.outer", row.names = F)
plot(studyarea_points)
#Vew(as.data.frame(states))
#Filter in only the states we want to see
#st_in <- c("AR", "AL", "TX", "OK", "MO", "IL", "TN", "MS", "LA", "GA", "FL")
#states <- states[states$STUSPS %in% st_in]
Convert to a data frame to extract centroids for labeling
state_df <- as.data.frame(crds(centroids(states))) # Get centroid coordinates
state_df$NAME <- states$NAME # Replace with the actual column from your shapefile
mcmcpath <- paste0(getwd(), "/mcmc")
Define the paths to the MCMC data and the desired output plot location
mcmcpath <-"E:/Armadillo_Analysis/Final_Armadillo_September_3/EEMS_armadillo_alfredo/aaaa/EEMS_Armadillo_not_FEEMS-output/barrier-schemeZ-nIndiv300-nSites3000-EEMS-nDemes200-chain1_res1"
plots <- make_eems_plots(mcmcpath, longlat = TRUE)
#########Plot the map properly
Create the plot with adjustments
Create the plot with adjustments
final_plot <- plots$mrates01 +
scale_fill_whitebox_c(palette = "muted", na.value = "gray10", guide = guide_colorbar(title= "log(m)")) +
geom_spatvector(data = states, inherit.aes = F, fill = NA, linewidth = 0.5) +
geom_spatvector(data = occ_shp, inherit.aes = F, fill = NA) +
Decreasing font size of state names
geom_text(data = state_df, aes(x = x, y = y, label = NAME), size = 3) +
annotation_north_arrow(pad_y = unit(17, "cm")) +
annotation_scale() +
theme_bw() +
Only "Longitude" and "Latitude" in labels
xlab("Longitude") +
ylab("Latitude") +
Increase axis text size (tick labels) and axis labels
theme(
axis.text.x = element_text(size = 14), # Increase X-axis tick label size
axis.text.y = element_text(size = 14), # Increase Y-axis tick label size
axis.title.x = element_text(size = 18, face = "bold"), # Increase X-axis label size
axis.title.y = element_text(size = 18, face = "bold") # Increase Y-axis label size
)
Save the plot as a high-resolution JPEG
ggsave("EEMS_final_heterozygosity_map.jpeg", final_plot, width = 12, height = 8, dpi = 600, device = "jpeg")
