#' ---
#' title: "Preparation of rasters with sdmpredictors"
#' author: "Miguel Gandra"
#' date: "Abril 2015"
#' update: "February 2025"
#' script: runs on Rstudio and R-terminal
#' ---


# Script to plot a distribution map of a chosen species from GBIF.
# http://www.gbif.org

#   ---> no prior file manipulation required
#   ---> automatically distinguishes points on land and sea, plotting them with different colors
#   ---> automatically saves a pdf file with the map on the chosen directory
#   ---> exports a csv file with the fetched coordinates
#   ---> prints a summary of the records


######################################################################################
# Automatically install required libraries  ##########################################
# (check http://tlocoh.r-forge.r-project.org/mac_rgeos_rgdal.html
# if rgeos and rgdal installation fails on a Mac)

#----------------------------------------------------------------------
#---------------------------------------------------------------------------#
# install and load packages
#---------------------------------------------------------------------------
check.packages <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

packages<-c("dismo", "XML","jsonlite", "graphics",
            "maps", "maptools", "tibble",
            "mime", "robis", "rfishbase", "dplyr")

check.packages(packages)

root <- setwd("~/Desktop/MBE")


######################################################################################
# Variables ##########################################################################

genus <- "XX" # write genus name here
species <- "XX" # write species name here
genusspecies <- paste(genus,species)

# Create results directory
results.directory <- "R.results"
dir.create(paste(root,results.directory,sep="/"))
results.directory <- paste(root,results.directory,sep="/")
getwd()


# "marginal.tolerance" 
# this is important to understand
# variable sets the margin at which records are considered to be in land.
# a positive value will consider records to be on land, even if they are up to x degrees outside the nearest coast.
# a negative value will constraint land records to the ones that are at least x degrees within the coastline.
marginal.tolerance <- 0

# user can automatically delete points in land by setting the "remove.land.points" variable to TRUE,
# or delete points in water by setting the "remove.ocean.points" variable to TRUE.
remove.land.points <- TRUE
remove.ocean.points <- FALSE
export.csv <- TRUE

#######################################################################################
# download GBIF them directly from the portal ############

gbif.data <- try(gbif(genus,species,geo=TRUE,removeZeros=TRUE))
coordinates<-data.frame(matrix(ncol=2, nrow=0))
coordinates <- data.frame(gbif.data$lon,gbif.data$lat)
coordinates <- na.omit(coordinates)

total <- nrow(coordinates)

# remove duplicate coordinates
dups <- duplicated(coordinates[,1:2])
dups <- dups[dups==TRUE]
coordinates <- unique(coordinates)
coordinates
records <- coordinates


#######################################################################################
# Set geographical area ###############################################################

x <- coordinates[,1]
y <- coordinates[,2]

xmin=round(min(x)-5,0)
xmax=round(max(x)+5,0)
ymin=round(min(y)-5,0)
ymax=round(max(y)+5,0)

xlim <- c(xmin, xmax)
ylim <- c(ymin, ymax)

########################################################################################
# Plot Map  ############################################################################
x_all<-coordinates$gbif.data.lon
y_all<-coordinates$gbif.data.lat

#plot an empty map
plot(x_all,y_all,xlab="Longitude",ylab="Latitude",asp=1,type="n")
#add continents mass
map(resol=0,fill=TRUE,col="gray95",border="gray95",add=TRUE)
#add coastline
map(resol=0,fill=FALSE,interior=FALSE,col="gray40",add=TRUE)
#add sampling points
points(x_all,y_all,col="blue",pch=16,cex=0.8)

########################################################################################
# Compute land and ocean points  #######################################################
data(wrld_simpl)

specie.pts <- SpatialPoints(coordinates, proj4string=CRS("+proj=longlat +datum=WGS84"))
projection(wrld_simpl) <- CRS("+proj=longlat +datum=WGS84")

min.distances <- c()

pts.on.land <- !is.na(over(specie.pts, wrld_simpl)$FIPS)

if (marginal.tolerance==0){
  x.ocean <- x[pts.on.land==FALSE]
  y.ocean <- y[pts.on.land==FALSE]
  x.land <- x[pts.on.land==TRUE]
  y.land <- y[pts.on.land==TRUE]
} else if (marginal.tolerance>0){
  distances <- gDistance(specie.pts,wrld_simpl,byid=TRUE)
  for (i in 1:length(specie.pts)) {min.distances[i] <- min(distances[,i])}
  x.ocean <- x[min.distances>marginal.tolerance]
  x.land <- x[min.distances<=marginal.tolerance]
  y.ocean <- y[min.distances>marginal.tolerance]
  y.land <- y[min.distances<=marginal.tolerance]
} else if (marginal.tolerance<0) {
  specie.pts <- specie.pts[pts.on.land==TRUE]
  mp <- map("world", xlim=c(xmin,xmax), ylim=c(ymin,ymax), col="gray60",  border="gray60", fill=TRUE, resolution=0)
  coastline <- cbind(mp$x, mp$y)[!is.na(mp$x),]
  coast.pts <- SpatialPoints(coastline, proj4string=CRS(proj4string(wrld_simpl)))
  distances <- gDistance(specie.pts,coast.pts,byid=TRUE)
  for (i in 1:length(specie.pts)) {min.distances[i] <- min(distances[,i])}
  x.ocean <- c(x[pts.on.land==FALSE],specie.pts@coords[min.distances<=abs(marginal.tolerance),1])
  y.ocean <- c(y[pts.on.land==FALSE],specie.pts@coords[min.distances<=abs(marginal.tolerance),2])
  x.land <- specie.pts@coords[min.distances>abs(marginal.tolerance),1]
  y.land <- specie.pts@coords[min.distances>abs(marginal.tolerance),2]
}


########################################################################################
# Plot ocean points ?  #################################################################

if (remove.ocean.points==FALSE){
  points(x.ocean, y.ocean, pch=21, col='green', bg='green', cex=0.5, lwd=1) 
}

########################################################################################
# Plot land points ?  ##################################################################

if (remove.land.points==TRUE){
  points(x.land, y.land, pch=21, col='red', bg='red', cex=0.5, lwd=1) 
}

legend("bottomleft", 
       legend = c("ocean", "land"), 
       col = c("green", "red"), 
       pch = 19, 
       bty = "n", 
       pt.cex = 2, 
       cex = 1.2, 
       text.col = "black", 
       horiz = F , 
       inset = c(0.1, 0.1))

########################################################################################
# Save pdf #############################################################################

#pdf.name <- paste(genus,'_',species,".pdf",sep="")
#pdf.file <- file.path(results.directory,pdf.name)
#dev.copy2pdf(file=pdf.file)


########################################################################################
# Export data as csv ###################################################################

if(export.csv==TRUE){
  occurrences <- data.frame(y.ocean,x.ocean)
  colnames(occurrences) <- c("lat","lon")
  csv.name <- paste(genus,'_',species,".csv",sep="")
  csv.file <- file.path(results.directory,csv.name)
  write.csv(occurrences, file=csv.file, row.names=FALSE)
}

########################################################################################
# Print summary table ##################################################################

r.summary <- data.frame(c(total,
                          length(dups),sum(length(x.ocean),length(x.land)), length(x.ocean),length(x.land)))
rownames(r.summary) <- c("Total","Duplicated","Plotted","Ocean","Land")
colnames(r.summary) <- c("RECORDS")
r.summary


########################################################################################
# END
########################################################################################