#' --- #' 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 ########################################################################################