###############################
##Code modified from Fletcher and Fortin 2019, Chapter 9
##Textbook - Spatial Ecology and Conservation Modeling
###############################
#gdistance pdf: https://mran.microsoft.com/snapshot/2017-02-04/web/packages/gdistance/vignettes/gdistance1.pdf
#load packages
library(raster)
## Loading required package: sp
library(gdistance) #for calculating connectivity metrics
## Loading required package: igraph
##
## Attaching package: 'igraph'
## The following object is masked from 'package:raster':
##
## union
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
## Loading required package: Matrix
##
## Attaching package: 'gdistance'
## The following object is masked from 'package:igraph':
##
## normalize
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.4 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::as_data_frame() masks tibble::as_data_frame(), igraph::as_data_frame()
## x purrr::compose() masks igraph::compose()
## x tidyr::crossing() masks igraph::crossing()
## x tidyr::expand() masks Matrix::expand()
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::groups() masks igraph::groups()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x dplyr::select() masks raster::select()
## x purrr::simplify() masks igraph::simplify()
## x tidyr::unpack() masks Matrix::unpack()
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(mapview)
#pull in RVA NLCD
lc<-raster("GrRVA_NLCD_2019.tif")
#label projection for later use
crs.land <- projection(lc)
#Richmond City Parks in need of connections
parks <- st_read("Parks.shp")
## Reading layer `Parks' from data source
## `/Users/joshuasementelli/Downloads/Connectivity/Parks.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 169 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 11748520 ymin: 3694416 xmax: 11805200 ymax: 3744511
## Projected CRS: NAD83 / Virginia South (ftUS)
parks<-st_transform(parks, crs.land)
#look at them
mapview(parks)
#identify a few parks to use for assessing landscape connectivity
#I selected Hugeunot FW, Larus park, Powhite park, and PP wetlands,
#b/c I am most familiar with these, but you can choose
#different ones if you like.
#I found OBJECTID by clicking on each park in mapview
parks<-filter(parks, OBJECTID == 186 | OBJECTID == 245 |
OBJECTID == 169 | OBJECTID == 241)
#crop NLCD to match this smaller extent where parks are located
lc_crop<-crop(lc, parks)
#identify the centroid of each park
park_cent <- st_centroid(parks, byid=T)
## Warning in st_centroid.sf(parks, byid = T): st_centroid assumes attributes are
## constant over geometries of x
#plot parks and centroids
plot(lc_crop)
plot(park_cent, add=T, color = "black", pch=20)
## Warning in plot.sf(park_cent, add = T, color = "black", pch = 20): ignoring all
## but the first attribute
#reclassify NLCD as resistance values:
lc_crop<-as.factor(lc_crop)
#Consider the resistance values (1-4, 1=low resistance and 4=high resistance)
#for each of the following cover types
#I was considering white tailed deer and/or red fox which I see in my
#neighborhood.
#you can consider other species of interest and change these values.
#typically this step would require a literature review so that it can be
#based on empirical data and expert opinion
#water = 4
#urban_low = 2
#urban_medLow = 2
#urban_Med = 3
#urban_high = 4
#barren = 2
#forest_dec = 1
#forest_eg = 1
#forest_mix = 1
#shrub = 1
#grass = 2
#pasture = 2
#crops = 2
#woody wetland = 1
#herbaceous wetland = 2
#create reclassification matrix based on above
reclass <- c(rep(4,1), rep(2,2), rep(3,1), rep(4,1),
rep(2,1) ,rep(1,4), rep(2,3), rep(1,1), rep(2,1))
reclass.mat <- cbind(levels(lc_crop)[[1]], reclass)
#we can use levels with raster package instead of unique for the terra package
#convert to factor and reclassify (these are different for raster than terra)
lc_crop<-as.factor(lc_crop)
resist <- reclassify(lc_crop, reclass.mat)
#plot resistance raster
plot(resist)
plot(park_cent, add=T, pch=20, color="black")
## Warning in plot.sf(park_cent, add = T, pch = 20, color = "black"): ignoring all
## but the first attribute
#create a conductance transition layer (inverse of resistance data)
#as mean of conductance
land_cond <- transition(1/resist, transitionFunction=mean, 8)
#also make a correction for the transition layer to account for diagonal
#connections from 8 neghborhood rule
#type=c for lcps; type = r for circuit - results same for this example
land_cond <- geoCorrection(land_cond, type="c", multpl=F)
?geoCorrection
#look at it by transforming the transition matrix to a raster
plot(raster(land_cond))#represents conductance
#geographic (Euclidean) distance matrix - as the crow flies
geo.dist <- pointDistance(park_cent, lonlat=FALSE)
?pointDistance
geo.dist <- as.dist(geo.dist)
?as.dist
## Help on topic 'as.dist' was found in the following packages:
##
## Package Library
## proxy /Library/Frameworks/R.framework/Versions/4.0/Resources/library
## stats /Library/Frameworks/R.framework/Versions/4.0/Resources/library
##
##
## Using the first match ...
?pointDistance
#convert park centroids to sp
park_cent_sp<-as_Spatial(park_cent)
#least-cost distance matrix
lc.dist <- costDistance(land_cond, park_cent_sp)
?costDistance
#commute distance matrix calculates the resistance distance (circuit theory)
#can take a long time with larger rasters
circuit.dist <- commuteDistance(land_cond, park_cent_sp)
#randomized shortest-paths distance matrix
# theta is the degree from which the path randomly deviates from the shortest path
##0 < theta < 20
#can take a long time on larger rasters
rSP.dist_t0001 <- rSPDistance(land_cond,
from=park_cent_sp,
to=park_cent_sp,
theta=0.0001)
?rSPDistance
#inspect
geo.dist
## 1 2 3
## 2 4390.036
## 3 2942.696 2854.240
## 4 4164.740 2411.543 4323.109
lc.dist
## 1 2 3
## 2 5844.601
## 3 3892.082 3298.356
## 4 5843.113 3048.600 5824.163
circuit.dist
## 1 2 3
## 2 194237.6
## 3 169481.5 182767.0
## 4 178508.2 154758.5 187499.1
rSP.dist_t0001
## [,1] [,2] [,3] [,4]
## [1,] 0.000 9652.480 6247.312 9273.446
## [2,] 9652.480 0.000 5767.176 4973.674
## [3,] 6247.312 5767.176 0.000 9348.709
## [4,] 9273.446 4973.674 9348.709 0.000
#take lower triangle of rSP.dist_t0001
rSP.dist_t0001.tri <- rSP.dist_t0001
rSP.dist_t0001.tri[upper.tri(rSP.dist_t0001.tri, diag=TRUE)] <- NA
#inspect
rSP.dist_t0001.tri
## [,1] [,2] [,3] [,4]
## [1,] NA NA NA NA
## [2,] 9652.480 NA NA NA
## [3,] 6247.312 5767.176 NA NA
## [4,] 9273.446 4973.674 9348.709 NA
#make data frame of distances
all.dist <- data.frame(Euclidean=as.vector(geo.dist),
lcd=as.vector(lc.dist),
circuit=as.vector(circuit.dist),
rSP=na.omit(as.vector(rSP.dist_t0001.tri)))
#correlation
pairs(all.dist[,1:4], pch = 19, lower.panel = NULL)
cor<-round(cor(all.dist),3)
library(corrplot)
## corrplot 0.90 loaded
corrplot(cor, method="number", type = "lower")
#are these measures correlated? Do relationships appear linear?
#any other observations?
##There does seem to be some correlation between the metrics calculated, especially between euclidean and lcd.
#now that we have quantified these metrics, now let's map them:
#First Least Cost Path:
#map shortest path between each park, start with path between park 1 and park 2
parks_lcp12 <- shortestPath(land_cond,
park_cent_sp@coords[1,], park_cent_sp@coords[2,],
output = "SpatialLines")
#if you want to get all paths from one park to the others, input range of values of the goal location
parks_lcp1 <- shortestPath(land_cond,
park_cent_sp@coords[1,], park_cent_sp@coords[1:4,],
output = "SpatialLines")
parks_lcp2 <- shortestPath(land_cond,
park_cent_sp@coords[2,], park_cent_sp@coords[1:4,],
output = "SpatialLines")
parks_lcp3 <- shortestPath(land_cond,
park_cent_sp@coords[3,], park_cent_sp@coords[1:4,],
output = "SpatialLines")
parks_lcp4 <- shortestPath(land_cond,
park_cent_sp@coords[2,], park_cent_sp@coords[4,],
output = "SpatialLines")
#plot these shortest paths
plot(resist, axes = F)
plot(park_cent_sp, pch = 20, add = T)
lines(parks_lcp1, col = "red")
lines(parks_lcp2, col = "red")
lines(parks_lcp3, col = "red")
lines(parks_lcp4, col = "red")
#notice that some of the lines are not lying right on top of each other...
#why do you think that might be?
##Lines are not exactly overlaying each because each pathway from either park 1 to park 2 or vice versa will not always have the some obstacles throughout the pathway. This wull cause variation in the shortest pathway lines.
#No Least Cost Corridors:
#mapping least cost corridors are better than paths (single line) which may not be reliable/used.
#this requires a few steps:
#first is to create the cumulative cost/resistance from each location
acc.cost1 <- accCost(land_cond, park_cent_sp@coords[1,])
acc.cost2 <- accCost(land_cond, park_cent_sp@coords[2,])
acc.cost3 <- accCost(land_cond, park_cent_sp@coords[3,])
acc.cost4 <- accCost(land_cond, park_cent_sp@coords[4,])
plot(acc.cost2)
#overlay the costs between two parks by adding the two layers and
#clip new layer based on specific quantile of cost
#(increasing the quantile will increase the width)
leastcost_corridor <- overlay(acc.cost1, acc.cost2, fun = function
(x, y){return (x + y)})
plot(leastcost_corridor)
#get lower 10% quantile of accumulated cost/resistance
quantile10 <- quantile(leastcost_corridor, probs = 0.10, na.rm = TRUE)
#create a blank raster that is the same size as leastcost_corridor
leastcost_corridor10 <- leastcost_corridor
values(leastcost_corridor10) <- NA
#subset only the area with 10% quatile and give it a value of 1
leastcost_corridor10[leastcost_corridor < quantile10] <- 1
#plot it, add points, and LCP line from above
plot(leastcost_corridor10,legend = F,axes = F)
points(park_cent_sp, col = "grey30")
lines(parks_lcp12, col = "red")
# we can now look at land cover along this predicted path for conservation planning
# using the extract function in raster package
#identify land-cover along the lcp
lcp.land <- raster::extract(lc_crop, parks_lcp12)
table(lcp.land)
## lcp.land
## 21 22 23 24 31 41 42 43 90
## 39 7 7 1 2 29 11 92 18
#identify land-cover along the least-cost corridor
corridor.land <- crop(lc_crop, leastcost_corridor10)
corridor.land <- mask(corridor.land, leastcost_corridor10)
plot(corridor.land, axes = F, legend = F)
## Warning in plot.window(...): "legend" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "legend" is not a graphical parameter
## Warning in title(...): "legend" is not a graphical parameter
## Warning in graphics::rasterImage(z, bb[1], bb[3], bb[2], bb[4], interpolate =
## interpolate, : "legend" is not a graphical parameter
table(as.vector(corridor.land))
##
## 21 22 23 24 31 41 42 43 52 81 90 95
## 941 220 129 17 15 337 90 765 1 4 250 1
#what are the most common land-cover types in the least-cost path and least-cost corridor?
#you can repeat this for the other set of two parks if you like
##Least cost path: mixed forest, developed land, and decidous forest
##Least-cost corridor: Mixed forest, Woody wetland, developed land, deciduous forest
#lastly, we can model flow using randomized shortest-path analysis
#must alter the tuning parameter - theta - to interpret this model
#when theta = 0, we are mapping a biased random walk
#movement is diffuse with no hard boundaries
passage.map_t0 <- passage(land_cond,
origin = park_cent_sp@coords[1,],
goal = park_cent_sp@coords[2,], theta = 0)
#alter theta to converge on a lcp
#as theta increases, the prob of passage become more well-deinfed andincrease around the lcp and lc corridor
passage.map_t0001 <- passage(land_cond,
origin = park_cent_sp@coords[1,],
goal = park_cent_sp@coords[2,], theta = 0.0001)
passage.map_t001 <- passage(land_cond,
origin = park_cent_sp@coords[1,],
goal = park_cent_sp@coords[2,], theta = 0.001)
passage.map_t002 <- passage(land_cond,
origin = park_cent_sp@coords[1,],
goal = park_cent_sp@coords[2,], theta = 0.05)
plot(passage.map_t0, axes=F, legend=F)
plot(passage.map_t0001, axes=F, legend=F)
plot(passage.map_t001, axes=F, legend=F)
plot(passage.map_t002, axes=F, legend=F)
#try different values of theta and see which one you think is best - explain
##The final plot displayed above with a theta parameter value of .05 I think bests fits the least cost path. With this value, pathways align where resistance is at the lowest.
#what additional data is needed to better answer the question about connectivity between these parks for wildlife?
##Predators present, slope of landscape could be factors that help answer questions on connectivity for wildlife.
#Do you have any other observations or comments to make about this exercise