mcglinnlab / shark-ray-div

Global patterns of shark and ray diversity
0 stars 1 forks source link

distorted indopacific #33

Closed EmmalineSheahan closed 2 years ago

EmmalineSheahan commented 6 years ago

So the indopacific when projected into CEA looks like this:

indopacific.pdf

which makes its mrd when plotted look like this:

indopacific_mrd.pdf

So I'm not certain if we should drop that one portion which is stretching the realm across the map, or if there's a better way to plot it so that doesn't happen. Here's the code I've been using to plot the realms and the mrd:

# creating and transforming realms
regions <- readOGR(dsn = "./data/continent", layer = "meow_ecos")
class(regions)
plot(regions)
regions$REALM
pull1 <- grep("Tropical Atlantic", regions$REALM)
tropical_atlantic <- regions[pull1,]
plot(tropical_atlantic)
pull2 <- grep("Central Indo-Pacific", regions$REALM)
central_indopacific <- regions[pull2,]
plot(central_indopacific)
tropical_atlantic <- spTransform(tropical_atlantic, CRS("+proj=cea +units=km"))
central_indopacific <- spTransform(central_indopacific, CRS("+proj=cea +units=km"))
phylo_bl1 <- compute.brlen(indopacific_tree, 1)
all_dist <- dist.nodes(phylo_bl1)
root_dist <- all_dist[length(indopacific_tree$tip.label) + 1, 
                      1:length(indopacific_tree$tip.label)]
tips_to_root <- data.frame(spp.name=indopacific_tree$tip.label, root_dist)

mrd_test_list <- vector("list", length = length(indopacific_mat)) 
for (i in 1:length(indopacific_mat)) {
  allspecies = colnames(indopacific_mat[[i]])
  mrd_test_list[[i]] = rep(0, nrow(indopacific_mat[[i]]))
  for(j in 1:nrow(indopacific_mat[[i]])) {
    sp_list = data.frame(spp.name = allspecies[indopacific_mat[[i]][j, ] == 1])
    if (nrow(sp_list) > 0) {
      root_dist_tot <- merge(sp_list, tips_to_root, sort = F)
      mrd_test_list[[i]][j] <- mean(root_dist_tot$root_dist)
    }
  }
}

load('./data/raster/indopacific_list.Rdata')
indopacific_mrd <- vector("list", length = 6)
pdf('./figures/indopacific_mrd.pdf')
for (i in 1:6) {
  mrd_raster <- indopacific_list[[i]]
  mrd_raster@data@values <- mrd_test_list[[i]]
  plot(mrd_raster)
  plot(continents, add = T, col = "black")
  indopacific_mrd[[i]] <- mrd_raster
}
dev.off()