Fragstats was designed to analyse whole landscapes with multiple classes of habitat types. However in reality most of the landscape metrics have rarely been used, due to difficulties in interpretation. In some cases measures of landscape diversity can be useful, but the most popular application of Fragstats has been to extract simpler landscape metrics for a single class of habitat.
Fragstats works reasonably well when the number of patches are small and the pixel resolution is low in comparison to patch size (more than around 100 pixels per patch). However when pixel resolution approaches patch size the analysis of small patches is artefactual as the perimeter measurements are not correct.
If the patches have been digitised as polygons it is much better to analyse them without conversion to a raster format.
Here are some idea about how to start in R.
First load the relevant packages. These are rgdal, rgeos. Spdep is used for some basic analysis of the distance matrix between patches.
rm(list = ls())
library(rgdal)
library(rgeos)
library(spdep)
If the data are held locally as a shapefile in a folder called “shapes” within the working directory this will load them.
setwd("/home/duncan/Dropbox/Public/Analyses/fragexample")
nf <- readOGR(dsn = "shapes", layer = "nfwood")
## OGR data source with driver: ESRI Shapefile
## Source: "shapes", layer: "nfwood"
## with 73 features and 1 fields
## Feature type: wkbPolygon with 2 dimensions
save(nf, file = "nf.rda")
This brings down some example data from my dropbox. The example consists of patches of habitat classified as broadleaf forest in good condition in the New Forest.
load(url("http://dl.dropboxusercontent.com/u/2703650/Analyses/fragexample/nf.rda"))
plot(nf, col = "green")
axis(1)
axis(2)
box()
grid()
In order to measure area and perimeter correctly it is important that the shapefile is projected in an appropriate CRS.
nf@proj4string
## CRS arguments:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy
## +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894
This shows that the CRS is OS grid. IN this case we can get the TOTAL area very simply.
gArea(nf)
## [1] 25409507
If the shapefile had been in geographic coordinates you would receive a warning as the units would not make sense. For example we can transform the polygons into WGS84 non projected coordinates this way.
nf2 <- spTransform(nf, CRS("+init=epsg:4326"))
gArea(nf2)
## Warning: Spatial object is not projected; GEOS expects planar coordinates
## [1] 0.003247
This produces a warning. However if the data are in geographic coordinates it is easy to transform them.
# Retransforming to OSGrid (or maybe UTM in Mexico .. look up appropriate
# EPSG code)
nf3 <- spTransform(nf, CRS("+init=epsg:27700"))
gArea(nf3)
## [1] 25409507
Using the gArea function by default gives the total area. Providing each row corresponds to a polygon we can easily get the area of each patch.
areas <- gArea(nf, byid = T)
Perimeter takes two simple steps using rgeos. We need the boundary first. The measure it. Read the manual in order to understand what happens if you have donut shaped polygons.
edge <- gBoundary(nf, byid = T)
perims <- gLength(edge, byid = T)
d <- data.frame(id = names(areas), areas, perims)
d <- d[order(d$areas, decreasing = T), ]
head(d)
## id areas perims
## 33 33 3966293 26042
## 17 17 1859594 18922
## 58 58 1729279 12483
## 31 31 1466547 13412
## 24 24 1399295 10352
## 8 8 1242838 11943
Now we have enough to get the basic fragstats.
To plot a histogram of areas.
hist(d$areas, col = "grey")
hist(log10(d$areas), col = "grey")
summary(d$perims)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8 1140 2090 3660 4810 26000
summary(d$areas)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 31300 118000 348000 402000 3970000
Areas are in square meters, so we may want to use hectares.
summary(d$areas/10000)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 3.1 11.8 34.8 40.2 397.0
A simple shape metric that takes values from 1 (perfectly compact) to infinty is derived by dividing the perimeter by the perimeter of a circle of the same area.
d$shape <- d$perims/(2 * pi * sqrt(d$areas/pi))
head(d)
## id areas perims shape
## 33 33 3966293 26042 3.689
## 17 17 1859594 18922 3.914
## 58 58 1729279 12483 2.678
## 31 31 1466547 13412 3.124
## 24 24 1399295 10352 2.469
## 8 8 1242838 11943 3.022
The simplest way to do this is to use a negative buffer. For example 100 m.
The smaller patches will have no core areas and so will return NA. RGeos produces a multipolygon for each id, so if there are multiple core areas within a patch you will get a single area. If you want to analyse the number of cores per patch you need some additional steps that I have not shown here.
cores <- gBuffer(nf, width = -100, byid = TRUE)
plot(nf, col = "darkgreen")
plot(cores, col = "red", add = T)
axis(1)
axis(2)
box()
grid()
c_area <- gArea(cores, byid = T)
c_area <- data.frame(id = names(c_area), c_area)
d <- merge(d, c_area, all = T)
d$corepercent <- (d$c_area/d$areas) * 100
d$edgepercent <- 100 - d$corepercent
head(d)
## id areas perims shape c_area corepercent edgepercent
## 1 0 3.129e+04 811.0 1.293 NA NA NA
## 2 1 5.024e+04 1675.1 2.108 NA NA NA
## 3 10 7.626e-01 12.3 3.972 NA NA NA
## 4 11 5.690e+04 1320.6 1.562 NA NA NA
## 5 12 2.266e+04 1014.6 1.901 NA NA NA
## 6 13 2.511e+05 2431.7 1.369 65336 26.02 73.98
This is always very challenging, as there are many decisions to be made regarding the nature of connectivity.
We can get a matrix of all edge to edge distances between patches using rgeos.
To demonstrate the fact that Rgeos calculates edge to edge distance look at this example.
p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
p2 = readWKT("POLYGON((2 0,3 1,4 0,2 0))")
plot(p1, xlim = c(0, 3), ylim = c(0, 4))
plot(p2, add = T)
axis(1)
axis(2)
box()
grid()
gDistance(p1, p2)
## [1] 1
Rgeos distance calculations will run slowly with many complex polygons so you will have to be quite patient if there are more than 100 patches. The reason for the slow calculation is that unless the algorithm is told to do otherwise it calcualates the distance between all points along the edge of each polygon and all points along the edge of all the others. This is brute force, but it does work in the end. PostGIS is orders of magnitude faster as it uses extremely efficient spatial indices. It is good when you have big data sets, but to use it you need to set it up and know how to import data.
mat <- gDistance(nf, byid = TRUE)
The matrix contains the distances between each patch and all the others This can be used in various ways.
If we set the distance for movement between patches to 1km we can look at the network of linked patches.
centroids <- coordinates(nf)
mat[mat > 1000] <- 0
aa <- mat2listw(mat)
str(summary(aa, zero.policy = T))
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 73
## Number of nonzero links: 208
## Percentage nonzero weights: 3.903
## Average number of links: 2.849
## 9 regions with no links:
## 3 11 15 20 34 35 42 50 60
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8
## 9 7 19 16 9 4 5 3 1
## 7 least connected regions:
## 12 41 49 55 66 67 69 with 1 link
## 1 most connected region:
## 8 with 8 links
##
## Weights style: M
## Weights constants summary:
## n nn S0 S1 S2
## M 64 4096 82454 112707947 633784046
## 'data.frame': 1 obs. of 5 variables:
## $ n : num 64
## $ nn: num 4096
## $ S0: num 82454
## $ S1: num 1.13e+08
## $ S2: num 6.34e+08
plot(nf, col = "darkgreen")
axis(1)
axis(2)
box()
grid()
plot(aa, centroids, add = T, pch = 21, cex = 0.2)
If your data come in raster form to start with then you can still use the vector approach shown above. You just need to convert the patches back into vector form.
To demonstrate this we can rasterise the vector layer.
library(raster)
## Attaching package: 'raster'
## The following object is masked from 'package:nlme':
##
## getData
## The following object is masked from 'package:MASS':
##
## area, select
## The following object is masked from 'package:Matrix':
##
## expand
## The following object is masked from 'package:sp':
##
## disaggregate
b <- bbox(nf)
r <- raster(nrow = 500, ncol = 5000, xmn = b[1, 1], xmx = b[1, 2], ymn = b[2,
1], ymx = b[2, 2])
rnf <- rasterize(nf, r)
## Found 73 region(s) and 90 polygon(s)
plot(rnf)
writeRaster(rnf, file = "test/test.tif")
## Error: filename exists; use overwrite=TRUE
You can see the artefacts that rasterisation introduces. It is not a good idea to rasterize a vector layer if you can avoid it. Fragstats used raster data as input because it was designed to work with classified imagery. If your study area is small it would usually be better to digitise habitat yourself with the aid of high resolution images and other sources of information.
There is not a simple “vectorize” function in R raster to convert clumps of pixels to polygons.
However I have found a trick that works very well for moderate sized rasters.
First cluster pixels into clumps. You will need igraph installed to do this.
clumps <- clump(rnf)
## Loading required package: igraph
## Attaching package: 'igraph'
## The following object is masked from 'package:raster':
##
## edge
plot(clumps)
This clumping mechanism could be very useful in itself as it is then easy to calculate fragstat like metrics on the clumps of pixels.
To vectorise the clumps I first extract the coordinates of the pixel centres.
id <- values(rnf)
crds <- coordinates(rnf)
crds <- crds[!is.na(id), ]
id <- id[!is.na(id)]
crds <- data.frame(id, crds)
coordinates(crds) <- ~x + y
Now the next step is a bit clunky computationally, but it works really well in most cases. First buffer out from the pixel centres. Then join the buffers. Finally buffer back in again. The difference between the two buffers should be about half the cell size. This produces quite a good restoration of the original vector, although some of the features will inevitably be lost.
buf1 <- gBuffer(crds, width = 100, byid = T)
buf2 <- gUnaryUnion(buf1, id)
buf3 <- gBuffer(buf2, width = -95, byid = T)
gArea(buf3)
## [1] 25775051
gArea(nf)
## [1] 25409507
We can write this out as shapefile.
buf3@proj4string <- nf@proj4string
n <- length(buf_final@polygons)
bb <- SpatialPolygonsDataFrame(buf3, data = data.frame(id = 1:n), match.ID = FALSE)
writeOGR(bb, "test", "test2", driver = "ESRI Shapefile")
