http://earthenginepartners.appspot.com/science-2013-global-forest
Packages needed.
library(rgdal)
library(rgeos)
library(raster)
The screenshot was very quickly georefefereced using only two points. In order to load the screenshot use the raster package to produce a “brick”.
https://dl.dropboxusercontent.com/u/2703650/Weblog/Defor/peten.png
r <- stack("peten.png")
The code below loads the results directly from my dropbox.
load(url("http://dl.dropboxusercontent.com/u/2703650/Weblog/Defor/defor.rda"))
The raster brick is called r. The total area can be calculated from the extent.
totalarea <- (r@extent@xmax - r@extent@xmin) * (r@extent@ymax - r@extent@ymin)/10000
totalarea
## [1] 166510
plotRGB(r, r = 1, b = 3, g = 2)
Now to look at the values in the green band.
plot(r[[2]])
So, we can clump all values over 50 to find forest patches.
cl <- clump(r[[2]] > 50)
plot(cl)
To produce a vector layer use the dilation and erosion trick. First extract the coordinates for the raster centroids and check for NA values.
id <- values(cl)
crds <- coordinates(cl)
crds <- crds[!is.na(id), ]
id <- id[!is.na(id)]
crds <- data.frame(id, crds)
coordinates(crds) <- ~x + y
Now use a buffer to join the points. After eroding the buffer we get polygons.
buf1 <- gBuffer(crds, width = 50, byid = T)
buf2 <- gUnaryUnion(buf1, id)
buf3 <- gBuffer(buf2, width = -50, byid = T)
We may want to remove fragments below 5 ha as there are a lot of very small patches. This should be considered when intepreting the statistics.
frags <- buf3[gArea(buf3, byid = T) > 50000]
plot(frags, col = "darkgreen")
axis(1)
axis(2)
box()
grid()
Percent of area forested
forestarea <- gArea(frags)/10000
100 * forestarea/totalarea
## [1] 16.48
area <- gArea(frags, byid = T)
edge <- gBoundary(frags, byid = T)
perims <- gLength(edge, byid = T)
d <- data.frame(id = names(frags), area, perims)
d$shape <- d$perims/(2 * pi * sqrt(d$area/pi))
d$area <- d$area/10000
head(d)
## id area perims shape
## 1 1 22.806 4863 2.872
## 3 3 5082.178 560335 22.173
## 7 7 97.820 17788 5.074
## 9 9 6.386 2382 2.659
## 11 11 17.785 4751 3.178
## 12 12 75.351 14931 4.852
d <- d[order(d$area, decreasing = T), ]
head(d)
## id area perims shape
## 3 3 5082.2 560335 22.17
## 2295 2295 3012.4 209040 10.74
## 691 691 1998.7 265413 16.75
## 707 707 1680.2 232236 15.98
## 1995 1995 770.1 123244 12.53
## 1035 1035 745.3 101270 10.46
hist(d$area, col = "grey")
hist(log10(d$area), col = "grey")
summary(d$area)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5 9 15 79 37 5080
cores <- gBuffer(frags, width = -100, byid = TRUE)
plot(frags, col = "darkgreen")
plot(cores, col = "red", add = T)
axis(1)
axis(2)
box()
grid()
Percent of total area that is core
100 * gArea(cores)/10000/totalarea
## [1] 4.277
c_area <- gArea(cores, byid = T)/10000
c_area <- data.frame(id = names(c_area), c_area)
d <- merge(d, c_area, all = T)
d$corepercent <- (d$c_area/d$area) * 100
d$edgepercent <- 100 - d$corepercent
d <- d[order(d$area, decreasing = T), ]
head(d)
## id area perims shape c_area corepercent edgepercent
## 245 3 5082.2 560335 22.17 1547.5 30.45 69.55
## 170 2295 3012.4 209040 10.74 1586.7 52.67 47.33
## 316 691 1998.7 265413 16.75 409.5 20.49 79.51
## 318 707 1680.2 232236 15.98 366.0 21.78 78.22
## 121 1995 770.1 123244 12.53 154.5 20.06 79.94
## 4 1035 745.3 101270 10.46 155.4 20.85 79.15
frags <- SpatialPolygonsDataFrame(frags, data = d, match.ID = FALSE)
writeOGR(frags, "shapefiles", "frags", driver = "ESRI Shapefile")
## Error: layer exists, use a new layer name