Reference and Disclaimer: For this exercise, I have used explanation and codes from the book Spatial Ecology and Conservation Modeling by Robert Fletcher & Marie-Josée Fortin. For analysis, I have used codes from (Landscape Genetic Data Analysis with R)[https://bookdown.org/hhwagner1/LandGenCourse_book/] and (How To Do Archaeological Science Using R)[https://benmarwick.github.io/How-To-Do-Archaeological-Science-Using-R/index.html]

I have got help from Carleisha Hanns and Tyler Putt in completing this assignment.

Load the necessary packages

rm(list = ls())

## First specify the packages of interest
packages = c("raster", "rgeos","rgdal","gdistance","igraph","fitdistrplus","fdrtool","ggplot2")

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)
par(mfrow = c(1, 1))

setwd("D:/UGA/Academic Course Work/Fall-2022/FANR-8400-AdvGIS/Assignments/Ex5")
fgdb <- "D:/UGA/Academic Course Work/Fall-2022/FANR-8400-AdvGIS/Assignments/Ex5/BoulderElk.gdb"
ogrListLayers(fgdb)

First - Creating cost layers using R

Load Vector and Raster Data

Load the vector data and evaluate their centroid

area_loc <- readOGR(fgdb,layer="elk_areas")
projection(area_loc)

#get the centroids of plots
area_loc_centroids <- gCentroid(area_loc, byid=T)
#area_loc_centroids@coords

Note: In order to calculate costs, all of the rasters need to have exactly the same resolution, dimensions, and co-ordinate locations (i.e., the cells need to match up perfectly. The method of resampling (bilinear vs nearest-neighbor) depends on the type of data. Nearest-neighbor is for categorical data whereas bilinear interpolation is for continuous data.

dem <- raster("dem_30")
road <- raster("./Cost_rasters/RoadCost.tif")
stream <- raster("./Cost_rasters/streamcost.tif")
trail <- raster("./Cost_rasters/trailcost.tif")

#create slope from elevation
slope <- terrain(dem, opt=c("slope"))
#make a multilayered file for extraction
layers <- stack(dem, slope, road, stream, trail) #merge the slope and aspect layers into a single raster object that holds all raster layers
names(layers) <- c("elev", "slope","road","stream","trail")
projection(layers)
rm(list=c("dem","road","stream","trail","slope"))

Plot the layers

plot(layers)

Plot the summer and winter sites - Base raster I took slope here

#plot
plot(layers$slope)
plot(area_loc,add=T)
points(area_loc_centroids, col="grey30", pch=20)

Creating Cost and transition layers

layers$elev <- (1/layers$elev)

tot_cost <- sum(layers[[1:5]])

tr.cost <- transition(tot_cost, 
                        transitionFunction = function(x) 1/mean(x), directions=8)
#Correct for geometric distortion
tr.cost1 <- gdistance::geoCorrection(tr.cost,type = "c",multpl=FALSE)

#get lcp
sumr_wntr_lcp <- shortestPath(tr.cost1, area_loc_centroids@coords[1,], area_loc_centroids@coords[2,], output="SpatialLines")
plot(raster::raster(tr.cost1))
lines(sumr_wntr_lcp, col="red", lwd=2)
points(area_loc_centroids, col="black", pch=20)

Next - Cost calculated in ArcGISpro

#Cost calculated in Arcpro and the final layer uploaded for further analysis in R
rm(list=c("tot_cost","tr.cost","tr.cost1","sumr_wntr_lcp"))
surface <- raster("./Cost_rasters/Surfacecost.tif")

plot(surface)
plot(area_loc,add=T)
points(area_loc_centroids, col="grey30", pch=20)

Effective least cost path

#create a conductance transition layer: inverse of resistance data
land_cond <- transition(surface, 
                        transitionFunction = function(x) 1/mean(x), 8)

#make correction; type=c for lcps; type = r for circuit (identical results for this example, so just use c)
land_cond <- geoCorrection(land_cond, type="c", multpl=F)

#get lcp
sumr_wntr_lcp <- shortestPath(land_cond, area_loc_centroids@coords[1,], area_loc_centroids@coords[2,], output="SpatialLines")

Plotting just the least cost path

#plotting this to see on the overall land data (for clarity purpose)
raster::plot(raster::raster(land_cond),main="Conductance")
plot(area_loc,add=T)
points(area_loc_centroids, col="grey30", pch=20)
lines(sumr_wntr_lcp)

Least-cost corridor

#get cumulative costs from each PA
sumr.cost <- accCost(land_cond, area_loc_centroids@coords[1,])
wntr.cost <- accCost(land_cond, area_loc_centroids@coords[2,])

#get least-cost corridor
leastcost_corridor <- overlay(sumr.cost, wntr.cost, fun=function(x, y){return(x + y)})

LCC

#get lower quantile
quantile10 <- quantile(leastcost_corridor, probs=0.10, na.rm=TRUE)

#make new truncated layer
leastcost_corridor10 <- leastcost_corridor
values(leastcost_corridor10) <- NA
leastcost_corridor10[leastcost_corridor < quantile10] <- 1 #truncate to identify corridor

#plot
par(mfrow=c(1,1))
plot(leastcost_corridor)
plot(area_loc, add=T)
plot(leastcost_corridor10, legend=F,axes=F, add=T, col="purple")
points(area_loc_centroids, col="grey30")
lines(sumr_wntr_lcp, col="red", lw=3)

Step-wise plot

par(mfrow=c(2,2))
plot(sumr.cost, legend=F, main="Summer Pref loc")
plot(wntr.cost, axes=F, main="Winter Pref loc")

plot(leastcost_corridor, legend=F, main="LC path")
plot(area_loc, add=T)
points(area_loc_centroids, col="grey30", pch=20)
lines(sumr_wntr_lcp, col="red", lw=3)

plot(leastcost_corridor, axes=F, main = "LC Corridor")
plot(area_loc, add=T)
plot(leastcost_corridor10, legend=F,  add=T, col="purple")
points(area_loc_centroids, col="grey30")
lines(sumr_wntr_lcp, col="red", lw=3)

par(mfrow=c(1,1))

Simple Arbitrary Additive Weight Sensitivity Model

Construct Weights - I am assigning arbitrary weights to classes of each raster variable. These weights should be based on educated guesses, empirical evidence, or to test a theory. The weights and splits should be based on regional literature. Known site locations obviously influence this understanding.

### Slope weighting Models###
slp_from <- c(0, 0.2, 0.4, 0.6, 0.8)
slp_to   <- c(0.2, 0.4, 0.6, 0.8, 10)
slp_wght <- c(50, 30, 15, 5, 0)
slp_rcls<- cbind(slp_from, slp_to, slp_wght) 
print(slp_rcls)
### Elevation weighting Models###
elev_from <- c(1500, 2000, 2500, 3000, 3500)
elev_to   <- c(2000, 2500, 3000, 3500,4100)
elev_wght <- c(50, 30, 15, 5, 0)
elev_rcls<- cbind(elev_from, elev_to, elev_wght) 
print(elev_rcls)

Note: In order to calculate costs, all of the rasters need to have exactly the same resolution, dimensions, and co-ordinate locations (i.e., the cells need to match up perfectly. The method of resampling (bilinear vs nearest-neighbor) depends on the type of data. Nearest-neighbor is for categorical data whereas bilinear interpolation is for continuous data.

layers
layers[["slope"]] <- reclassify(layers[["slope"]], slp_rcls)
layers[["elev"]] <- reclassify(layers[["elev"]], elev_rcls)
# plot(layers2)

model_sum <- sum(layers[[1:5]])
# plot(model_sum)
coords <- coordinates(model_sum)
x <- data.frame(x = coords[,"x"], y = coords[,"y"], value = as.vector(model_sum))
ggplot(x, aes(x = x, y = y, fill = value)) +
  geom_raster(hjust = 0, vjust = 0) +
  theme_minimal() +
  viridis::scale_fill_viridis() +
  labs(title = "Summed Sensitivity Weights")

Note: We can redo the entire analysis with this new layer

Circuitscape

To be able to use the data in circuitscape, the polygon data is converted to rasters. Then surface cost raster and elk location rasters were both converted into .asc format. the data is then put in circuitscape.

From Circuitscape