# This script builds distribution maps
# 
library(plyr)
library(rgeos)
library(rgdal)
library(ggplot2)
library(gridExtra)
library(maps)

# theme -------------------------------------------------------------------

theme_nothing2 <- function (base_size = 12, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace% 
    theme(rect = element_rect(fill = "transparent", 
                              colour = NA,
                              color = NA, 
                              size = 0, 
                              linetype = 0), 
          line = element_blank(), 
          plot.title = element_text(size=14, face='italic',hjust = 0.5),
          axis.text = element_blank(),  
          axis.title = element_blank(),
          panel.background = element_blank(), axis.ticks.length = grid::unit(0,"lines"), legend.position = "none", 
          panel.spacing = grid::unit(c(0,0, 0, 0), "lines"), 
          plot.margin = grid::unit(c(0.1,0.1, 0.1,0.1), "lines"))
}




# load master table -------------------------------------------------------

setwd("~/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/Final files/to submit/DataS1") 
setwd("C:/Users/Fernando/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/NeoBiota/2nd submission/DataS1")
dir()

master.yes <- read.csv("Citizen science.csv", head=T)
master.no <- read.csv("All records.csv", head=T)

head(master.no)


# table with species names
setwd("~/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/Final files/to submit/DataS1")
sp.list <- read.csv("Species List.csv", head=T)

# load GIS ---------------------------------------------------------------

setwd("~/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/Final files/to submit/DataS1/GIS")
setwd("C:/Users/Fernando/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/NeoBiota/2nd submission/DataS1/GIS")

# Iberia
iberia = readOGR(dsn=".", layer="iberia")
iberia <- spTransform(iberia, CRS("+init=epsg:4326"))
iberia <- iberia[,"CNTRY_NAME"]

# grid
grid = readOGR(dsn=".", layer="UTM_iberia_May2016")
grid <- spTransform(grid, CRS("+init=epsg:4326"))


# Urban areas
setwd("~/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/Final files/to submit/GIS")
iberian.urban <- readOGR(dsn=".", layer="Urban areas")

plot(iberia)
plot(iberian.urban, add=TRUE, col=1)


# prepare table
coord <- data.frame(coordinates(grid))
names(coord) <- c("x","y")
grid@data <- cbind(grid@data,coord)

head(grid@data)

master.yes1 <- merge(master.yes, grid@data, by.x="UTM", by.y="CODIGO")
master.no1 <- merge(master.no, grid@data, by.x="UTM", by.y="CODIGO")

head(master.yes1)
head(master.no1)


grid2 <- fortify(grid)


# Plots PDF -------------------------------------------------------------

# template

# Iberia
PI <- spTransform(iberia, CRS("+init=epsg:23029"))
PI <- gBuffer(PI, width=500) # if some points are outside PI
PI <- spTransform(PI, CRS("+init=epsg:4326"))
PI$name <- "Iberia"

PI@data$id = rownames(PI@data)
PI.points = fortify(PI, region="id")
PI.df = join(PI.points, PI@data, by="id")

# urban areas
urban <- gBuffer(iberian.urban, width = 0, byid = TRUE)
urban@data$id = rownames(urban@data)
urban.points = fortify(urban, region="id")
urban.df = join(urban.points, urban@data, by="id")


# plot template
(PI.plot <- ggplot(PI.df) + 
    aes(long,lat,group=group) + 
    geom_polygon(fill="white") +
    geom_path(color="black") +
    geom_polygon(data=urban.df, aes(group=group), fill="grey") +
    theme_minimal() +
    coord_map()) 



# loop over species 

setwd("~/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/Final files/to submit/MAPS")

for(myClass in c("Amphibians", "Reptiles","Birds","Mammals")){
  
  pdf(paste(myClass,".pdf", sep=""),  onefile = TRUE)
  
  for(sp in names(master.no1)[grep(myClass, names(master.no1))]){
    
    m.sub <- master.no1[,c(sp,"UTM", "y", "x")]
    m.sub <- m.sub[m.sub[,1]==1,]
    
    class <- strsplit(sp, "_")[[1]][1]
    sp.name <- strsplit(sp, "_")[[1]][2]
    sp.name <- gsub("[.]"," ",sp.name) 
    
    common.name <- unique(sp.list[sp.list$Species_in.use==sp.name, "common_name"])
    subtitle <- as.character(sp.name)
    
    (p <- PI.plot +
      theme_nothing2() +
      # geom_polygon(aes(x = long, y = lat, group = group), data = grid2,
      #              fill = 'white', color="grey", size = 0.1, alpha=.5) +
      geom_point(data=m.sub, aes(x=x, y=y, group=m.sub[,1]), 
                 size =1.6, fill="darkred", shape=22) +
      labs(title=sp.name) + 
      annotate("text", x = 3, y = 36.6, label = common.name, size = 4.5, hjust=1) +
      #annotate("text", x = 3, y = 37, label = class, size = 4.5, hjust=1) +
      annotate("text", x = 3, y = 36.2, label = paste("N = ", nrow(m.sub)), size = 4.5, hjust=1))
    
    
    if(sp %in% names(master.yes1)){
      
      m.sub2 <- master.yes1[,c(sp,"UTM", "x", "y")]
      m.sub2 <- m.sub2[m.sub2[,1]==1,]
      
      (p <- p +
         geom_point(data=m.sub2, aes(x=x, y=y, group=m.sub2[,1]), 
                   size =1.3, fill="white", shape=22))
    }
    
    grid.arrange(p)
    
  }
  
  dev.off()
}


# Only present in Citizen Science platforms

outs <- names(master.yes1)[!(names(master.yes1) %in% names(master.no1))]


pdf("Only citizen science.pdf",  onefile = TRUE)

for(sp in outs){
  
  print(sp)
  
  m.sub <- master.yes1[,c(sp,"UTM", "y", "x")]
  m.sub <- m.sub[m.sub[,1]==1,]
  
  class <- strsplit(sp, "_")[[1]][1]
  sp.name <- strsplit(sp, "_")[[1]][2]
  sp.name <- gsub("[.]"," ",sp.name) 
  
  common.name <- unique(sp.list[sp.list$Species_in.use==sp.name, "common_name"][1])
  
  subtitle <- sp.name
  
  (p <- PI.plot +
      theme_nothing2() +
      # geom_polygon(aes(x = long, y = lat, group = group), data = grid2,
      #              fill = 'white', color="grey", size = 0.1, alpha=.5) +
      geom_point(data=m.sub, aes(x=x, y=y, group=m.sub[,1]), 
                 size =3, fill="red", shape=22) +
      labs(title=sp.name) + 
      annotate("text", x = 3, y = 36.6, label = common.name, size = 4.5, hjust=1) +
      annotate("text", x = 3, y = 37, label = class, size = 4.5, hjust=1) +
      annotate("text", x = 3, y = 36.2, label = paste("N = ", nrow(m.sub)), size = 4.5, hjust=1))
  
  
  
  grid.arrange(p)
  
}

dev.off()



# plot totals -------------------------------------------------------------

head(master.no1)

# plot template
(PI.plot <- ggplot(PI.df) + 
   aes(long,lat,group=group) + 
   geom_polygon(fill="white") +
   geom_path(color="black") +
   geom_polygon(data=urban.df, aes(group=group), fill="grey") +
   theme_minimal() +
   coord_map()) 


# Amphibians
my.sp <- names(master.no1)[grep("Amphibians", names(master.no1))]

myd <- master.no1[,c("UTM", "x", "y", my.sp)]
myd$total <- rowSums(myd[,-c(1:3)])

(p.amphibians <- ggplot(PI.df) +
    aes(long,lat,group=group) + 
    geom_polygon(fill="grey") +
    geom_path(color="black") +
    theme_minimal() +
    coord_map() +
    geom_point(data=myd, aes(x=x, y=y, group=myd[,1], 
                                    colour = factor(total)), 
               shape = 16) +
    scale_colour_viridis_d("Number of\nspecies") +
    geom_path(color="white") +
    labs(title="Amphibians", x="Longitude", y="Latitude") + 
    annotate("text", x = 3, y = 36.2, label = paste("N = ", ncol(myd)-4), size = 4.5, hjust=1))

# Reptiles
my.sp <- names(master.no1)[grep("Reptiles", names(master.no1))]

myd <- master.no1[,c("UTM", "x", "y", my.sp)]
myd$total <- rowSums(myd[,-c(1:3)])

(p.reptiles <- ggplot(PI.df) +
    aes(long,lat,group=group) + 
    geom_polygon(fill="grey") +
    geom_path(color="black") +
    theme_minimal() +
    coord_map() +
    geom_point(data=myd, aes(x=x, y=y, group=myd[,1], 
                             colour = (total)), 
               shape = 16) +
    scale_colour_viridis_c("Number of\nspecies") +
    geom_path(color="white") +
    labs(title="Reptiles", x="Longitude", y="Latitude") + 
    annotate("text", x = 3, y = 36.2, label = paste("N = ", ncol(myd)-4), size = 4.5, hjust=1))

# Birds
my.sp <- names(master.no1)[grep("Birds", names(master.no1))]

myd <- master.no1[,c("UTM", "x", "y", my.sp)]
myd$total <- rowSums(myd[,-c(1:3)])

(p.birds <- ggplot(PI.df) +
    aes(long,lat,group=group) + 
    geom_polygon(fill="grey") +
    geom_path(color="black") +
    theme_minimal() +
    coord_map() +
    geom_point(data=myd, aes(x=x, y=y, group=myd[,1], 
                             colour = log(total+1)), 
               shape = 16) +
    scale_colour_viridis_c("Number of\nspecies\n(log scale)") +
    geom_path(color="white") +
    labs(title="Birds", x="Longitude", y="Latitude") + 
    annotate("text", x = 3, y = 36.2, label = paste("N = ", ncol(myd)-4), size = 4.5, hjust=1))


# Mammals
my.sp <- names(master.no1)[grep("Mammals", names(master.no1))]

myd <- master.no1[,c("UTM", "x", "y", my.sp)]
myd$total <- rowSums(myd[,-c(1:3)])

(p.mammals <- ggplot(PI.df) +
    aes(long,lat,group=group) + 
    geom_polygon(fill="grey") +
    geom_path(color="black") +
    theme_minimal() +
    coord_map() +
    geom_point(data=myd, aes(x=x, y=y, group=myd[,1], 
                             colour = (total)), 
               shape = 16) +
    scale_colour_viridis_c("Number of\nspecies") +
    geom_path(color="white") +
    labs(title="Mammals", x="Longitude", y="Latitude") + 
    annotate("text", x = 3, y = 36.2, label = paste("N = ", ncol(myd)-4), size = 4.5, hjust=1))




setwd("~/OneDrive - Universidade de Lisboa/Documents/Exotic data paper/Final files/to submit/DataS1/MAPS")
ggsave(grid.arrange(p.amphibians,p.reptiles,p.birds,p.mammals),
       filename="Distribution patterns.png", units="cm", width = 20, height = 20)