---
title: 'Climate Change SDMs: Student Template'
author: "Grant Foster"
date: "2025-10-29"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Installing packages: This should only need to be done once on each computer, but including it here just in case. 
```{r, eval=FALSE}
install.packages("knitr")
install.packages("geodata")
install.packages("dismo")
install.packages("dplyr")
install.packages("terra")
install.packages("rJava")
install.packages("viridisLite")
install.packages("predicts")
install.packages("scales")
```

Loading packages: Loading the packages in that contain functions we'll use for grabbing data & training our models
```{r}
library(knitr)
library(geodata)
library(dismo)
library(dplyr)
library(terra)
library(predicts)
library(rJava)
library(dplyr)
library(viridisLite)
library(scales)
```

Loading in climatic data for us to make predictions with

```{r}
#Grabbing a polygon of SC for mapping
SC_boundary_vect <- vect("SC_Boundary.shp")

#Loading in our climatic data + doing a bit of filtering
raster.comb <- terra::rast("SC_Current_predictors_PCA.tif")
raster.comb <- raster::stack(raster.comb)
terra::crs(raster.comb) <- "EPSG:4326"
names(raster.comb) <- paste0("PC", 1:27)
project.comb126 <- terra::rast("SC_Future126_predictors_PCA.tif")
project.comb126 <- raster::stack(project.comb126)
terra::crs(project.comb126) <- "EPSG:4326"
names(project.comb126) <- paste0("PC", 1:27)
project.comb585 <-  terra::rast("SC_Future585_predictors_PCA.tif")
project.comb585 <- raster::stack(project.comb585)
names(project.comb585) <- paste0("PC", 1:27)
terra::crs(project.comb585) <- "EPSG:4326"
```

Defining our main analysis function. You'll have to run this once so that "trainMaxent" appears in your environment (to the right), which we can use later. 
```{r}
trainMaxent <- function(Occ, split=0.75, sp="SpeciesName"){
    start_time <- Sys.time()
    Occ <- dplyr::select(Occ, "species", "decimalLongitude", "decimalLatitude") #Select just the Occ columns
    Occ['cells']<-terra::cellFromXY(object= raster.comb[[1]], xy = Occ[,2:3]) #Grab the occurrence cells
    #creating a 200km2 buffer to select the psudo-absence points
    Occ <- unique(Occ) #remove redundant records
    OccSf<-sf::st_as_sf(x=Occ[,2:3], coords = c("decimalLongitude",'decimalLatitude'), 
                        crs=sp::CRS("EPSG:4326")) #crs=sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
    OccBuffer<-sf::st_buffer(OccSf, dist=100000) #Create buffer objects with ~100km buffers
    OccBuffer<-sf::st_union(OccBuffer) #Merge buffer objects into one
    print("Successfully created buffers")

      ##Masking the buffer so that it does not occur in the ocean2deg
      OccBufferVect <- terra::vect(OccBuffer)
      OccBufferVect <- terra::project(OccBufferVect, crs(raster.comb))
      OccBufferRas <- terra::rasterize(x = OccBufferVect, y=terra::rast(raster.comb))
      OccBufferRas<- raster::mask(x = OccBufferRas, mask = terra::rast(raster.comb)[[1]])
    
    ##Selecting the potential pseudo absence sites
    AllPseudo <- terra::as.points(OccBufferRas)
    AllPseudo <- raster::extract(x= terra::rast(raster.comb)[[1]], AllPseudo[,1:2],cell=T)
    AllPseudo2 <- AllPseudo[!(AllPseudo$cell %in% Occ$cells), ] #remove true presence from potential psuedo points
    
    #Actually sample psuedoabsences
    PseudoA <- sample(x = AllPseudo$cell, size = nrow(Occ)) #Sample a number of psuedoA equal to # of occurances
    PseudoA<-as.data.frame(raster::xyFromCell(object = terra::rast(raster.comb)[[1]], cell = PseudoA)) #Make a df
    
    #Now, train a Maxent model
    Occurrences<-cbind(Occ[,2:3], PseudoA[,1:2]) #bind into 4 cols; left is real, right is associated psuedo
    colnames(Occurrences)<-c('Long','Lat','x','y')
    

    id.training <- sample(1:nrow(Occurrences), round(split*nrow(Occurrences),0)) #create test-train split (split=percent train)
    training <- dismo::prepareData(x=raster.comb, p=Occurrences[id.training,c("Long", "Lat")], b=Occurrences[id.training,c("x", "y")])  #Use dismo prepare data to set up df
    testing <-dismo::prepareData(x=raster.comb, p=Occurrences[-id.training,c("Long", "Lat")], b=Occurrences[-id.training,c("x", "y")])    
    training<-na.omit(training)
    testing<-na.omit(testing)
      
    colnames(Occurrences)<-rep(c('Long','Lat'),2) #Change names
    testingPresPA<-rbind.data.frame(Occurrences[-id.training,1:2],Occurrences[-id.training,3:4])
    print("Successfully formatted Occ data")
          ##Maxent
    Sys.setenv(NOAWT=TRUE)
    Maxent.Model <-dismo::maxent(x = as.data.frame(training[,-1]), p = as.data.frame(training[,1])) #Train Maxent Model
    predicts::predict(object=Maxent.Model, x=raster.comb)
    Maxent.0  <-dismo::predict(object = Maxent.Model, x = raster.comb, args='outputformat=cloglog') #Predict suitability according to trained model
    crs(Maxent.0) <- crs(SC_boundary_vect)
    Maxent.eval <-dismo::evaluate(p = testing[testing[,"pb"]==1,-1], a = testing[testing[,"pb"]==0,-1], model = Maxent.Model)

    # Mask the raster to the shape of SC
    Maxent_masked <- mask(terra::rast(Maxent.0), SC_boundary_vect)
    Maxent_masked <- crop(Maxent_masked, c(-84, -78, 32, 35.5))
    cPer <- 100*length(Maxent_masked[Maxent_masked>0.75])/terra::ncell(Maxent_masked)
    cPer <- round(cPer, 3)
     print("Successfully trained model on current climate")
    pal <- viridisLite::viridis(100)
    
    Maxent.future126  <-dismo::predict(object = Maxent.Model, x = project.comb126, args='outputformat=cloglog')
    crs(Maxent.future126) <- crs(SC_boundary_vect)
    Maxent.future126_cut <- mask(terra::rast(Maxent.future126), SC_boundary_vect)
    Maxent.future126_cut <- crop(Maxent.future126_cut, c(-84, -78, 32, 35.5))
    fPer126 <- 100*length(Maxent.future126_cut[Maxent.future126_cut>0.75])/terra::ncell(Maxent.future126_cut)
    fPer126 <- round(fPer126, 3)
    print("Successfully trained model on future climate (SSP 126)")
    
        Maxent.future585  <-dismo::predict(object = Maxent.Model, x = project.comb585, args='outputformat=cloglog')
    crs(Maxent.future585) <- crs(SC_boundary_vect)
    Maxent.future585_cut <- mask(terra::rast(Maxent.future585), SC_boundary_vect)
    Maxent.future585_cut <- crop(Maxent.future585_cut, c(-84, -78, 32, 35.5))
    fPer585 <- 100*length(Maxent.future585_cut[Maxent.future585_cut>0.75])/terra::ncell(Maxent.future585_cut)
    fPer585 <- round(fPer585, 3)
    print("Successfully trained model on future climate (SSP 585)")
    
#Plotting the outputs    
pdf(file= paste(sp, "_modeloutputs.pdf", sep=""), width = 14, height = 6)


# 3 panels side-by-side
op <- par(mfrow = c(1, 3),
          mar = c(4, 4, 3.5, 2),   # bottom, left, top, right
          oma = c(0, 0, 2, 0))     # outer margin (for an overall title if you want)

# Use a shared color scale from 0 to 1
zlim_shared <- c(0, 1)

### Panel 1: Current
terra::plot(Maxent_masked,
            main = paste0(sp, ": Current Conditions \nModel AUC: ", round(Maxent.eval@auc, 3)),
            xlab = "Decimal Longitude",
            ylab = "Decimal Latitude",
            zlim = zlim_shared,
            col = pal)
# Add state boundary (works if SC_boundary_vect is either SpatVector or Spatial)
try(terra::plot(terra::vect(SC_boundary_vect), add = TRUE, lwd = 2, border = "black"), silent = TRUE)
# Add occurrences
points(x = Occ$decimalLongitude, y = Occ$decimalLatitude, pch = 16, cex = 0.5, col=scales::alpha("darkorange", 0.5))
mtext(paste0(cPer, "% of state area has suitability > 0.75"), side = 3, line = -1, cex = 0.85)

### Panel 2: Future SSP126 (2061-2080)
terra::plot(Maxent.future126_cut,
            main = paste0("Future Conditions (2061-2080, SSP1-2.5)\n% area of state > 0.75: ", fPer126, "%"),
            xlab = "Decimal Longitude",
            ylab = "",
            zlim = zlim_shared,
            col = pal)
try(terra::plot(terra::vect(SC_boundary_vect), add = TRUE, lwd = 2, border = "black"), silent = TRUE)
points(x = Occ$decimalLongitude, y = Occ$decimalLatitude, pch=16, cex = 0.5, col=scales::alpha("darkorange", 0.5))
mtext(paste0(fPer126, "% of state area has suitability > 0.75"), side = 3, line = -1, cex = 0.85)

### Panel 3: Future SSP585 (2061-2080)
terra::plot(Maxent.future585_cut,
            main = paste0("Future Conditions (2061-2080, SSP5-8.5)\n% area of state > 0.75: ", fPer585, "%"),
            xlab = "Decimal Longitude",
            ylab = "",
            zlim = zlim_shared,
            col = pal)
try(terra::plot(terra::vect(SC_boundary_vect), add = TRUE, lwd = 2, border = "black"), silent = TRUE)
points(x = Occ$decimalLongitude, y = Occ$decimalLatitude, pch = 16, cex = 0.5, col=scales::alpha("darkorange", 0.5))
mtext(paste0(fPer585, "% of state area has suitability > 0.75"), side = 3, line = -1, cex = 0.85)


# Optional overall title (in outer margin)
mtext(paste0("SDM projections for ", sp, " (Current and 2061–2080 futures)"),
      outer = TRUE, cex = 1.2, font = 2)

# restore graphics and close device
par(op)
dev.off()

  end_time <- Sys.time()
  elapsed <- end_time - start_time
  
  cat("Time elapsed:",
      round(as.numeric(elapsed, units = "secs"), 2), "seconds (",
      round(as.numeric(elapsed, units = "mins"), 2), "minutes )\n")


print(paste("Done, output saved at ", paste(sp, "_modeloutputs.pdf", sep=""), sep=""))
}
```

Alright, now that we have our climate data and our analysis function, we're ready to apply it to some real data. Let's practice using some Cottonmouth viper (Agkistrodon piscivorus) occurrence data I got from GBIF. I've saved its occurrence points in a file called "Cottonmouth.csv".
```{r}
cottonmouth <- read.csv(file="Cottonmouth.csv", header=T, sep="\t") #load data

cottonmouth <- dplyr::filter(cottonmouth, is.na(decimalLatitude)==F) #remove points without latitude (no spatial info)
cottonmouth <- cottonmouth[sample(1:nrow(cottonmouth), 3000, replace=F),] #Optional: subset to a smaller # of points

trainMaxent(Occ=cottonmouth, split=0.75, sp="Cottonmouth2") #Run the model on the "cottonmouth" dataframe
```

If that ran, there should now be a pdf in your "directory" (the folder where this file is) called "Cottonmouth_modeloutputs.pdf." You can open it up to see the actual results of the model.

## Insert Your Species 1 Code Below
```{r}




```

## Insert Your Species 2 Code Below
```{r}



```


