Landscape Metrics
Edge Density
Load the Raster
Naming the raster “Wooded.cvr” for this analysis.
#install.packages("sp")
library(sp)
#install.packages("raster")
library(raster)
# Loaction of the raster on computer (Don't forget the ".tif" at the end)
Wooded.cvr = raster("C:/Users/humph173/Documents/Michigan_State/FragStats/MiCCAP_Reclass_Wooded_Cover1/MiCCAP_Reclass_Wooded_Cover1.tif")
# here's what it looks like
plot(Wooded.cvr)
# Take a look at it's properties:
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)Reducing Raster to Save Time
The raster is really big - it has 574,523,292 cells. To speed things up for this example, I’m going to reduce the raster (making the cells 6 times bigger). You won’t have to do this step, but, you’ll have to wait on the full size rasters to run.
Wooded.cvr = aggregate(Wooded.cvr, fact = 6, method = "ngb") # Take a look at it's properties after aggregation
Wooded.cvr## class : RasterLayer
## dimensions : 4138, 3857, 15960266 (nrow, ncol, ncell)
## resolution : 180, 180 (x, y)
## extent : 424785, 1119045, 2092305, 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 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## data source : C:\Users\humph173\Documents\Michigan_State\FragStats\Wooded_ag.tif
## names : Wooded_ag
## values : 0, 1 (min, max)# To see how different it looks
plot(Wooded.cvr)
Open the spatialEco package
Here’s the guide: https://cran.r-project.org/web/packages/spatialEco/spatialEco.pdf
#install.packages("spatialEco")
library(spatialEco)## spatialEco 1.1-0## Type se.news() to see new features/changes/bug fixes.Calculate Edge Density
This function calculates Edge Density for the whole landscape (entire raster), so is kind of time consumming. In this case, I name the output/results as “MyMetrics.”
If you’re only interested in locations where cats have been identified, there’s a faster way to get the values within a buffer distance of observation points (see below for example).
MyMetrics = focal.lmetrics(Wooded.cvr, #Name of raster to analyze
bkg = 0, #Background value in raster
land.value = 1, #Raster value for wooded areas
metric = "landscape.shape.index", #The landscape metric that you want is called "landscape.shape.index"
latlong = FALSE) #Your raster is projected in Albers equal-area (meters)# Details for MyMetrics
MyMetrics## class : RasterLayer
## dimensions : 4138, 3857, 15960266 (nrow, ncol, ncell)
## resolution : 180, 180 (x, y)
## extent : 424785, 1119045, 2092305, 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 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## data source : C:\Users\humph173\Documents\Michigan_State\FragStats\MyMetrics.tif
## names : MyMetrics
## values : 1, 3 (min, max)Plot and View Results
Two packages are used for plotting the results (levelplot and latticeExtra).
First, getting state boundaries.
#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("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(MyMetrics))) # Raster version
Michigan.r = rasterize(Michigan,
Wooded.cvr,
field = 0,
background = NA)Then, plotting Edge Density. Recall, this raster has been reduced some.
#install.packages("RColorBrewer")
library(RColorBrewer)
#install.packages("lattice")
library(lattice)
#install.packages("latticeExtra")
library(latticeExtra)
#install.packages("RasterVis")
library(rasterVis)
rng = seq(1, 3, 0.01)
cr = colorRampPalette(c("yellow", "red"),
bias = 1, space = "rgb")
M1RE = levelplot(MyMetrics,
margin = FALSE,
sub = "Edge Density \n(meters of edge per square meter)",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(1, 3, 0.01),
labels=list(at=c(1, 2, 3)),
labels=c("0", "1", "2", "3")),
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
M1RE +
latticeExtra::layer(sp.polygons(Michigan, col = "black", lwd = 2)) 
Get Fake Site Points (12 random points)
Fake_sites = sampleRandom(Michigan.r, 12, na.rm = TRUE, sp = TRUE)
#Name Sites
Fake_sites$Site = paste("Site", 1:12, sep=".")
rownames(Fake_sites@data) = Fake_sites$Site
#Details
Fake_sites## class : SpatialPointsDataFrame
## features : 12
## extent : 478515, 1022475, 2162235, 2709435 (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 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 2
## names : MichiganR, Site
## min values : 0, Site.1
## max values : 0, Site.9View points:
M1RE +
latticeExtra::layer(sp.polygons(Michigan, col = "black", lwd = 2)) +
latticeExtra::layer(sp.polygons(Fake_sites, col = "black", pch = 19, cex=2))
Calculate Edge Density in Buffer (5 km buffer)
Buffer_km = land.metrics(
Fake_sites, #Points
Wooded.cvr, #Raster
bkgd = NA, #Value to ignore
bw = 5000, #Buffer width
metrics = c("landscape.shape.index", "total.area", "n.patches"), # Stats for Edge Density, area, wooded patches
echo = FALSE) #Don't need to automatically print results
#Check output
Buffer_km[["0"]]## class n.patches total.area landscape.shape.index
## Site.1 0 21 2170800 5.058824
## Site.2 0 21 7290000 4.400000
## Site.3 0 24 9201600 5.058824
## Site.4 0 34 33598800 10.692308
## Site.5 0 1 63180000 5.426966
## Site.6 0 4 648000 2.333333
## Site.7 0 1 32400 1.000000
## Site.8 0 2 2494800 2.055556
## Site.9 0 25 1846800 5.187500
## Site.10 0 61 13446000 10.512195
## Site.11 0 35 25174800 8.375000
## Site.12 0 31 2980800 6.000000#Add to points
Fake_sites@data = cbind(Fake_sites@data, as.data.frame(Buffer_km[["0"]]))
# Take a look
Fake_sites@data## MichiganR Site class n.patches total.area landscape.shape.index
## Site.1 0 Site.1 0 21 2170800 5.058824
## Site.2 0 Site.2 0 21 7290000 4.400000
## Site.3 0 Site.3 0 24 9201600 5.058824
## Site.4 0 Site.4 0 34 33598800 10.692308
## Site.5 0 Site.5 0 1 63180000 5.426966
## Site.6 0 Site.6 0 4 648000 2.333333
## Site.7 0 Site.7 0 1 32400 1.000000
## Site.8 0 Site.8 0 2 2494800 2.055556
## Site.9 0 Site.9 0 25 1846800 5.187500
## Site.10 0 Site.10 0 61 13446000 10.512195
## Site.11 0 Site.11 0 35 25174800 8.375000
## Site.12 0 Site.12 0 31 2980800 6.000000Compare Sites
library(ggplot2)##
## Attaching package: 'ggplot2'## The following object is masked from 'package:latticeExtra':
##
## layerggplot(Fake_sites@data, aes(factor(Site), landscape.shape.index)) +
geom_bar(stat="identity") +
ylab(expression("Edge Density "~m/m^2)) +
xlab("Site") +
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),
axis.text.x = element_text(face="bold", size=12, vjust=1),
plot.title = element_text(face="bold", size=18, hjust=0.5))