Edge Density for Selected Sites
(250 - 1500 meter buffers)
Set Directory
This is the folder with all the data.
setwd("~/Michigan_State/FragStats")Load the Raster
Naming the raster “Wooded.cvr” for this analysis.
#install.packages("sp")
library(sp)
#install.packages("raster")
library(raster)
#install.packages("lattice")
library(lattice)
#install.packages("latticeExtra")
library(latticeExtra)## Loading required package: RColorBrewer#install.packages("RColorBrewer")
library(RColorBrewer)
#install.packages("RasterVis")
library(rasterVis)
#install.packages("rgeos")
library(rgeos)## rgeos version: 0.3-28, (SVN revision 572)
## GEOS runtime version: 3.6.1-CAPI-1.10.1 r0
## Linking to sp version: 1.3-1
## Polygon checking: TRUE#install.packages("dplyr")
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Loaction of the raster on computer (Don't forget the ".tif" at the end)
Wooded.cvr = raster("./MiCCAP_Reclass_Wooded_Cover1/MiCCAP_Reclass_Wooded_Cover1.tif")
#Raster Info
Wooded.cvr## class : RasterLayer
## dimensions : 24826, 23142, 574523292 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 424785, 1119045, 2092365, 2837145 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## data source : C:\Users\humph173\Documents\Michigan_State\FragStats\MiCCAP_Reclass_Wooded_Cover1\MiCCAP_Reclass_Wooded_Cover1.tif
## names : MiCCAP_Reclass_Wooded_Cover1
## values : 0, 1 (min, max)Getting state boundaries (for later plotting).
#install.packages("maps")
library(maps)
#install.packages("maptools")
library(maptools)## Checking rgeos availability: TRUEStates = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"),
function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
States = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(LL84))
#Add a dataframe
pid = sapply(slot(States, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(States),
row.names = pid)
States = SpatialPolygonsDataFrame(States, p.df)
# Change to projection of raster
Michigan = spTransform(subset(States, ID == 21), CRS(proj4string(Wooded.cvr)))
#Also, getting Counties
Counties = map("county",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(Counties$names, ":"),
function(x) x[1])
Counties = map2SpatialPolygons(Counties, IDs = IDs,
proj4string = CRS(LL84))
#Add a dataframe
pid = sapply(slot(Counties, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(Counties),
row.names = pid)
Counties = SpatialPolygonsDataFrame(Counties, p.df)
Counties = spTransform(Counties, CRS(proj4string(Wooded.cvr)))
Counties = gBuffer(Counties, byid = TRUE, width = 0)
Counties = crop(Counties, Michigan)Load Site Locations
Turing coordinates to points and then re-projecting coordinates to match the raster.
Sites = read.csv("./camera_locations.csv")
#take a peek
head(Sites)## Site.ID Latitude Longitude
## 1 GS02 43.28716 -84.27902
## 2 GS03 43.27228 -84.30752
## 3 GS04 43.27293 -84.35764
## 4 GS05 43.28710 -84.35831
## 5 GS06 43.25835 -84.37042
## 6 GS07 43.27619 -84.41386# Convert to points
Site.pnts = SpatialPointsDataFrame(Sites[,c("Longitude", "Latitude")], Sites)
proj4string(Site.pnts) = proj4string(States)
#Change projection
Site.pnts = spTransform(Site.pnts, CRS(proj4string(Wooded.cvr)))
#Adding new coordinates to dataframe (attribute table)
Site.pnts@data = Site.pnts@data %>%
mutate(pLong = Site.pnts@coords[,1],
pLat = Site.pnts@coords[2])
#Check data
head(Site.pnts)## Site.ID Latitude Longitude pLong pLat
## 1 GS02 43.28716 -84.27902 944100.0 942019.9
## 2 GS03 43.27228 -84.30752 942019.9 942019.9
## 3 GS04 43.27293 -84.35764 937993.3 942019.9
## 4 GS05 43.28710 -84.35831 937746.2 942019.9
## 5 GS06 43.25835 -84.37042 937167.5 942019.9
## 6 GS07 43.27619 -84.41386 933442.0 942019.9Quick Map
rng = seq(0, 1, 0.5)
mCols = brewer.pal(11, "RdYlBu")[-6]
cr0 = (colorRampPalette((mCols))(n = 256))
cr = colorRampPalette(c("tan", "green"),
bias = 1, space = "rgb")
M.plt = levelplot(Wooded.cvr,
margin = FALSE,
xlab = "Site Locations",
ylab = NULL,
col.regions = cr, at = rng,
colorkey = FALSE,
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Counties, col = "black", lwd = 1)) +
latticeExtra::layer(sp.polygons(Michigan, col = "black", lwd = 2)) +
latticeExtra::layer(sp.polygons(Site.pnts, col = "red", pch=19, cex = 1))
M.plt
As above, but Zoom in some
Buff.pnt = gBuffer(Site.pnts, byid = TRUE, width = 10000)
Focal.cnty = crop(Counties, extent(Buff.pnt))
Focal.r = crop(Wooded.cvr, extent(Buff.pnt))
M.plt2 = levelplot(Focal.r,
margin = FALSE,
xlab = "Site Locations",
ylab = NULL,
col.regions = cr, at = rng,
colorkey = FALSE,
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1.25)) +
latticeExtra::layer(sp.polygons(Focal.cnty, col = "black", lwd = 1)) +
latticeExtra::layer(sp.polygons(Site.pnts, col = "red", pch=19, cex = 1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(),
offset = c(930000, 2257000),
scale = 25000)})
M.plt2
Calculate Edge Density in Buffers
Here’s the guide: https://cran.r-project.org/web/packages/spatialEco/spatialEco.pdf
#install.packages("spatialEco")
library(spatialEco)
Buffer_250m = land.metrics(
Site.pnts, #Points
Wooded.cvr, #Raster
bkgd = 0, #Background value
bw = c(250), #Buffer widths
metrics = c("landscape.shape.index"), #Stat for Edge Density
echo = FALSE) #Don't need to automatically print results
Buffer_500m = land.metrics(
Site.pnts,
Wooded.cvr,
bkgd = 0,
bw = c(500),
metrics = c("landscape.shape.index"),
echo = FALSE)
Buffer_750m = land.metrics(
Site.pnts,
Wooded.cvr,
bkgd = 0,
bw = c(750),
metrics = c("landscape.shape.index"),
echo = FALSE)
Buffer_1000m = land.metrics(
Site.pnts,
Wooded.cvr,
bkgd = 0,
bw = c(1000),
metrics = c("landscape.shape.index"),
echo = FALSE)
Buffer_1500m = land.metrics(
Site.pnts,
Wooded.cvr,
bkgd = 0,
bw = c(1500),
metrics = c("landscape.shape.index"),
echo = FALSE)
#Add outputs to points
Site.pnts@data = Site.pnts@data %>%
mutate(m250 = Buffer_250m[,2],
m500 = Buffer_500m[,2],
m750 = Buffer_750m[,2],
m1000 = Buffer_1000m[,2],
m1500 = Buffer_1500m[,2])
#Quick check
head(Site.pnts)Quick Plot
library(ggplot2)##
## Attaching package: 'ggplot2'## The following object is masked from 'package:latticeExtra':
##
## layerlibrary(reshape2)
Pltm = melt(Site.pnts@data[, c(1, 6:10)], "Site.ID")
Pltm$value = as.numeric(Pltm$value)
ggplot(Pltm, aes(factor(Site.ID), value)) +
geom_bar(stat="identity") +
ylab(expression("Edge Density "~m/m^2)) +
xlab("Site") +
facet_wrap(~variable, ncol = 5) +
theme_classic() +
theme(axis.title.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size=20),
axis.text.y = element_text(face="bold", size=12),
strip.text = element_text(face="bold", size = 18),
axis.text.x = element_blank(),
plot.title = element_text(face="bold", size=12, hjust=0.5))## Warning: Removed 1 rows containing missing values (position_stack).
Write Data to CSV
#write.csv(Site.pnts@data, "./Edge_density_Jul07_2018.csv")