NOTE: This is a legacy script that will not work when RGDAL is deprecated in fall 2023.

PostProcessingModels.R

#Post-processing complementary script
#========================================================
#author: Hannah Owens
#date: 13 March, 2018

#Set the working directory
#====================================
setwd("~/Dropbox/ENMSeminar/Labs:Homeworks/Lab9/Data/"); #Change to YOUR working directory

library(raster);

#SOME MODEL POST-PROCESSING JOLLITY   
##################################################################

#Get environmental data used to train your model
setwd("Present/");
envtList <- list.files(pattern = ".asc")
envtStack <- stack(envtList);

#Get occurrence points
setwd("~/Dropbox/ENMSeminar/Labs:Homeworks/Lab9/Data/");
codPoints <- read.csv("cleanCodPointsSPOCC.csv");

# Thresholding
#====================================

# Example thresholding to 95% suitabilty score. You can use any number, though.
gadusMorhuaModel <- raster("Gadus_morhuaModel.asc");
suitabilityScores <- extract(gadusMorhuaModel, codPoints[,2:3]); #Extract suitability scores at occurrences
suitabilityScores <- suitabilityScores[complete.cases(suitabilityScores)] #Clean out no data values
threshold <- sort(suitabilityScores, decreasing = T)[round(length(suitabilityScores)*.95,0)] #Takes 5th percentile suitability score
m <- c(0, threshold, 0,  threshold, 1, 1); #Create a list of thresholds in the following format: (Minimum, maximum, valueToChange)
rclmat <- matrix(m, ncol=3, byrow=TRUE); #Changes the threshold list to a matrix
codDist <- reclassify(gadusMorhuaModel, rcl = rclmat); #Reclassify suitability surface to presence/absence
#NOTE: codDist can be saved like any other raster file (refer to Lab 2);

# Calculating the extent of ENM-suitable habitat
#====================================

#Calculating the extent of ENM-suitable habitat
cellAreas <- area(codDist)*codDist #Makes a raster of cell area
plot(cellAreas);
cellAreaMeasures <- rasterToPoints(cellAreas); #Adds up the area of each cell
paste("The area of the species' distribution is ", sum(cellAreaMeasures[,3]), " km2.", sep=""); #prints an answer


#Basic niche characterization
#====================================

# Extracting descriptive statistics from ENM-predicted presences
gadusMorhuaSuit <- envtStack * codDist;
names(gadusMorhuaSuit) <- names(envtStack);
gadusMorhuaPts <- cbind(codPoints, extract(envtStack, codPoints[,2:3]));
gadusMorhuaPts <- gadusMorhuaPts[complete.cases(gadusMorhuaPts),4:7]

#Shows difference in inferred niche characteristics
cellStats(envtStack, stat = 'max'); #Characteristics across training region
cellStats(gadusMorhuaSuit, stat = 'max'); #Characteristics at ENM-predicted presences
noquote(summary(gadusMorhuaPts[,5:9], digits = 2)[6,]); #Characteristics at occurrence points