ESPM 137, Lab 6: Landscape Metrics
#! This lab is worth 10 pts and is due at 5pm on Friday, 10/9 !#
Overview and Goals
In this lab, we will calculate landscape metrics for a forested landscape in Humboldt County. Like many parts of California, this landscape is a ‘working landscape’ where there are several competing land uses. One of these is the longstanding timber harvest industry that has been thoroughly regulated in California, including Humboldt County. Timber harvesters must file a timber harvest plan with the county and records of timber extracted from different parcels. From this information, we can reconstruct data layers for the areas of timber harvest, which we’ll work with today. Another emerging land use is cannabis production, which has accelerated in this region of the state since 2000. Cannabis grows have for a long time, at least until the legalization of recreational cannabis use in California, been a more clandestine land use activity. Nevertheless, we have been able to identify areas of cannabis production from satellite photography, from which we can construct GIS raster layers that we will also work with today. We’ll use the SDMTools package in R to calculate a variety of landscape metrics to analyze patterns of forest cover change on this landscape. The SDMTools package is based on the FragStats program that was originally developed to examine changes in forested landscapes and has a long history in landscape ecology.
Goals– 1: Calculate landscape metrics for forested landscapes in Humboldt County, CA. 2: Examine the effects of land use conversion on the study landscape. 3: Quantify the effects of cannabis grows and timber harvest on the landscape forest structure in the study area.
Set up the R session
Load required packages…
Set markdown preferences…
Part 1: Raster Data and Proportions of Habitat
First, let’s load the data by reading in rasters for forest cover, timber harvested, and cannabis grows.
forest2000 <- raster("forest021.tif")
forest2000## class : RasterLayer
## dimensions : 1314, 2320, 3048480 (nrow, ncol, ncell)
## resolution : 0.0005, 0.0005 (x, y)
## extent : -124.4476, -123.2876, 39.91097, 40.56797 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : /cloud/project/forest021.tif
## names : forest021
## values : 0, 1 (min, max)timber <- raster("timber_harvest.tif")
timber## class : RasterLayer
## dimensions : 1314, 2320, 3048480 (nrow, ncol, ncell)
## resolution : 0.0005, 0.0005 (x, y)
## extent : -124.4476, -123.2876, 39.91097, 40.56797 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : /cloud/project/timber_harvest.tif
## names : timber_harvest
## values : 0, 1 (min, max)grows <- raster("grows634.tif")
grows## class : RasterLayer
## dimensions : 1314, 2320, 3048480 (nrow, ncol, ncell)
## resolution : 0.0005, 0.0005 (x, y)
## extent : -124.4476, -123.2876, 39.91097, 40.56797 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : /cloud/project/grows634.tif
## names : grows634
## values : 0, 1 (min, max)Let’s plot the forest raster…
plot(forest2000, main="Forest Cover")
And then the rasters of cannabis grows and timber harvest…
par(mfrow=c(1,2))
plot(grows, main="Cannabis Grows")
plot(timber, main="Timber Harvest")
We can calculate the proportion of the habitat (p) covered by forest, using the prop() function.
prop(x) x: a raster layer on which to calculate the proportion of habitat type 1
forest.p <- prop(forest2000)
forest.p## [1] 0.7515825Question 1 (1 pt)
1a) In the box below, calculate the proportions of cannabis grows and timber harvest on the landscape.
cannabis.p <- prop(grows)
timber.p <- prop(timber)
cannabis.p## [1] 0.002671212timber.p ## [1] 0.018723321b) Relatively speaking, how much of the landscape is covered by areas of timber harvest compared to cannabis grows?
Cannabis grow corresponds to less than one percent (0.27%) and timber harves to 1.8%. Lanscape covered by cannabis growth is much less than areas of timber harvest.
Part 2: Calculating Landscape Metrics for Forest Cover
Next, let’s calculate some landscape metrics for our original landscape - that means for the forest021.tif raster (our forest object). That raster is a snapshot of forest cover for this study landscape in 2000. We can use the ClassStat() function to calculate a huge list of standard landscape metrics.
ClassStat(map, bkgd = NA, latlon = FALSE) map: a matrix or raster layer with categorical values representing different habitat types bkgd: the background value for which statistics will not be calculated latlon: boolean value representing if the data is geographic
The background value (the non-habitat) is coded as 0 on this raster, so we set bkgd=0.Because we have geographic data with real spatial locations, we set latlon=TRUE.
# Calculate landscape metrics
forest2000.metrics <- ClassStat(forest2000, bkgd=0, latlon=TRUE)
# Display results
showMetrics(forest2000.metrics)## [,1]
## class 1.0000
## n.patches 1657.0000
## total.area 2047345376.5329
## prop.landscape 0.7516
## patch.density 0.0000
## total.edge 8991831.5500
## edge.density 0.0033
## landscape.shape.index 50.0940
## largest.patch.index 0.6486
## mean.patch.area 1235573.5525
## sd.patch.area 43558804.7891
## min.patch.area 2353.4891
## max.patch.area 1766857948.9376
## perimeter.area.frac.dim 0.0088
## mean.perim.area.ratio 0.0547
## sd.perim.area.ratio 0.0216
## min.perim.area.ratio 0.0036
## max.perim.area.ratio 0.0831
## mean.shape.index 1.2551
## sd.shape.index 1.1059
## min.shape.index 1.0000
## max.shape.index 38.5442
## mean.frac.dim.index 1.0371
## sd.frac.dim.index 0.0456
## min.frac.dim.index 1.0001
## max.frac.dim.index 1.3423
## total.core.area 1618081797.0903
## prop.landscape.core 0.5940
## mean.patch.core.area 976512.8528
## sd.patch.core.area 35443527.3469
## min.patch.core.area 0.0000
## max.patch.core.area 1438110759.9901
## prop.like.adjacencies 0.8978
## aggregation.index 94.7167
## lanscape.division.index 0.5763
## splitting.index 2.3598
## effective.mesh.size 1154366727.2123
## patch.cohesion.index 9.9778We can retrieve any of the individual metrics by calling them with their names after the $ sign. Let’s look at some examples…
forest2000.metrics$total.edge # Show us the value calculated for the total edge metric## [1] 8991832forest2000.metrics$prop.landscape # Show us the total proportion of forest cover## [1] 0.7515764forest2000.metrics$prop.landscape.core # Show us the total proportion of forest core habitat## [1] 0.5939946Question 2 (0.5 pts)
What is the proportion of forest edge habitat on our study landscape? >>Because the total proportion of forest cover is 0.7515764 and the total proportion of forest core habitat is 0.5939946, the proportion of forest edge habitat is 0.1575818. That is about 16% of our study landscape.
0.7515764-0.5939946## [1] 0.1575818Part 3: Creating New Raster of Forest Cover
Our goal is to analyze how the metrics of forest cover have changed due to deforestation from cannabis grows and timber harvest. The forest2000 raster represents the forest cover in 2000, and the grows and timber rasters represent areas of cannabis grows and timber harvest observed in 2013. On the forest raster, values of 1 represent forest cover, and values of 0 represent non-forest. On the grows and timber rasters, values of 1 represent area that is now a cannabis grow or that was harvested for timber, and values of 0 represent areas that weren’t modified by these activities. You can also refer to the following table to see what the raster values represent.
Raster 0 1
forest non-forest forest cover grows non-grow cannabis grow timber not harvested timber harvested
To measure how the landscape metrics have changed due to cannabis grows and timber harvest, we need to create new rasters that remove those areas from the raster of forest cover in 2000.
Question 3 (2 pts)
In the chunk below, create a new raster named grows.removed by removing areas of cannabis grows from the forest cover represented on the forest2000 raster. And create a new raster named timber.removed by removing areas of harvested timber from the forest cover represented on the forest2000 raster. We haven’t shown you how to ‘remove’ habitat from a raster, per se, but we think you can figure out how to do so. You don’t need to use a fancy new function or anything we haven’t shown you before.
grows.removed= forest2000 - grows
timber.removed= forest2000 - timberQuestion 4 (1.5 pts)
4a) Check your results by calculating the proportion of forest cover on the new rasters in the chunk below.
cannabis.pro <- prop(grows.removed)
timber.pro <- prop(timber.removed)
cannabis.pro## [1] 0.7489113timber.pro## [1] 0.73267664b) Is this result what you expected? Why or why not? >>Yes, the proportion of forest cover in cannabis grow is grater than forest cover of harvest timber. Timber harvesting causes more deforestation than cannabis grow. The result was what I was expecting because the landscape proportion of cannabis grow that we calculated in question 1 was much smaller than timber harvest.
Part 4: Calculating New Landscape Metrics for Forest Cover
Now that we have rasters for forest cover with cannabis grows removed and with timber harvest areas removed, we can calculate the landscape metrics for forest cover following these landscape modifications.
Question 5 (1.5 pt)
In the chunk below, calculate the landscape metrics using the ClassStat() function for the new rasters you created. Store the results as objects named grows.metrics and timber.metrics so that we can compare them later.
#Original
forest2000.metrics$prop.landscape## [1] 0.7515764# Calculate grow & harvest metrics
grows.metrics <- ClassStat(grows.removed, bkgd=0, latlon=TRUE)
timber.metrics <- ClassStat(timber.removed, bkgd=0, latlon=TRUE)
# Grow & harvest results
grows.metrics$prop.landscape## [1] 0.7489025timber.metrics$prop.landscape## [1] 0.7327153Now let’s view the results…
metrics <- rbind(forest2000=forest2000.metrics, cannabis=grows.metrics, timber=timber.metrics)
showMetrics(metrics)## forest2000 cannabis timber
## class 1.0000 1.0000 1.0000
## n.patches 1657.0000 1669.0000 1772.0000
## total.area 2047345376.5329 2040061493.4407 1995966225.1023
## prop.landscape 0.7516 0.7489 0.7327
## patch.density 0.0000 0.0000 0.0000
## total.edge 8991831.5500 9339725.4700 10263366.0000
## edge.density 0.0033 0.0034 0.0038
## landscape.shape.index 50.0940 52.1292 57.8983
## largest.patch.index 0.6486 0.6460 0.6316
## mean.patch.area 1235573.5525 1222325.6402 1126391.7749
## sd.patch.area 43558804.7891 43228513.1998 41022107.4865
## min.patch.area 2353.4891 2353.4891 2353.4891
## max.patch.area 1766857948.9376 1759746997.9152 1720576259.2370
## perimeter.area.frac.dim 0.0088 0.0092 0.0103
## mean.perim.area.ratio 0.0547 0.0548 0.0559
## sd.perim.area.ratio 0.0216 0.0216 0.0217
## min.perim.area.ratio 0.0036 0.0038 0.0044
## max.perim.area.ratio 0.0831 0.0831 0.0831
## mean.shape.index 1.2551 1.2570 1.2529
## sd.shape.index 1.1059 1.1455 1.2415
## min.shape.index 1.0000 1.0000 1.0000
## max.shape.index 38.5442 40.6315 45.8219
## mean.frac.dim.index 1.0371 1.0372 1.0361
## sd.frac.dim.index 0.0456 0.0457 0.0455
## min.frac.dim.index 1.0001 1.0001 1.0001
## max.frac.dim.index 1.3423 1.3473 1.3590
## total.core.area 1618081797.0903 1589780093.0850 1507692200.3362
## prop.landscape.core 0.5940 0.5836 0.5535
## mean.patch.core.area 976512.8528 952534.5075 850842.0995
## sd.patch.core.area 35443527.3469 34635344.0963 31898746.7874
## min.patch.core.area 0.0000 0.0000 0.0000
## max.patch.core.area 1438110759.9901 1410218483.2856 1338143596.4697
## prop.like.adjacencies 0.8978 0.8937 0.8814
## aggregation.index 94.7167 94.4899 93.7997
## lanscape.division.index 0.5763 0.5797 0.5982
## splitting.index 2.3598 2.3788 2.4880
## effective.mesh.size 1154366727.2123 1145159536.1656 1094873757.4449
## patch.cohesion.index 9.9778 9.9777 9.9776We can also calculate the changes between the original forest2000 metrics and the metrics calculated on your new rasters to quantify how much change in the forest metrics each of these land-uses has caused.
diff.metrics <- rbind(ChangeFromCannabis=grows.metrics-forest2000.metrics, ChangeFromTimber=timber.metrics-forest2000.metrics)
showMetrics(diff.metrics[,2:38])## ChangeFromCannabis ChangeFromTimber
## n.patches 12.0000 115.0000
## total.area -7283883.0922 -51379151.4306
## prop.landscape -0.0027 -0.0189
## patch.density 0.0000 0.0000
## total.edge 347893.9200 1271534.4500
## edge.density 0.0001 0.0005
## landscape.shape.index 2.0352 7.8043
## largest.patch.index -0.0026 -0.0170
## mean.patch.area -13247.9124 -109181.7776
## sd.patch.area -330291.5893 -2536697.3026
## min.patch.area 0.0000 0.0000
## max.patch.area -7110951.0224 -46281689.7006
## perimeter.area.frac.dim 0.0004 0.0015
## mean.perim.area.ratio 0.0000 0.0011
## sd.perim.area.ratio 0.0000 0.0001
## min.perim.area.ratio 0.0002 0.0007
## max.perim.area.ratio 0.0000 0.0000
## mean.shape.index 0.0019 -0.0022
## sd.shape.index 0.0396 0.1356
## min.shape.index 0.0000 0.0000
## max.shape.index 2.0873 7.2777
## mean.frac.dim.index 0.0002 -0.0009
## sd.frac.dim.index 0.0001 0.0000
## min.frac.dim.index 0.0000 0.0000
## max.frac.dim.index 0.0051 0.0168
## total.core.area -28301704.0053 -110389596.7542
## prop.landscape.core -0.0104 -0.0405
## mean.patch.core.area -23978.3453 -125670.7533
## sd.patch.core.area -808183.2506 -3544780.5596
## min.patch.core.area 0.0000 0.0000
## max.patch.core.area -27892276.7044 -99967163.5204
## prop.like.adjacencies -0.0041 -0.0164
## aggregation.index -0.2268 -0.9170
## lanscape.division.index 0.0034 0.0219
## splitting.index 0.0190 0.1282
## effective.mesh.size -9207191.0467 -59492969.7674
## patch.cohesion.index 0.0000 -0.0002Question 6 (2 pts)
6a) Consider the values we estimated for the number of patches, mean patch area, and total edge on each of our rasters. Also look at the landscape division index - this index is the probability than any two cells (drawn at random) are in different patches. In general terms, explain the effects that cannabis production and timber harvesting each have on the forest cover in this study landscape >>The change in the number of forest cover patches was greater in the timber raster, therefore, there was an increase in the total forest edge. The mean patch area of forest cover decreased more on the timber raster than on the cannabis grow raster. The landscape division index was also greater on the timber raster. Both cannabis growth and timber harvesting have an impact on forest cover, however, timber harvesting causes much more deforestation(smaller patches, greater forest edge habitat and less total forest cover) and fragmentation(landscape division index) than cannabis grow.
#forest, cannabis, timber
metrics$n.patches## [1] 1657 1669 1772metrics$mean.patch.area ## [1] 1235574 1222326 1126392metrics$total.edge## [1] 8991832 9339725 10263366metrics$lanscape.division.index## [1] 0.5763346 0.5797108 0.5982241###changes between the original forest2000 metrics and the metrics calculated (cannabis vs timber)
diff.metrics$n.patches## [1] 12 115diff.metrics$mean.patch.area ## [1] -13247.91 -109181.78diff.metrics$total.edge## [1] 347893.9 1271534.4diff.metrics$lanscape.division.index## [1] 0.003376189 0.0218895146b) There is one metric where cannabis grows and timber harvest cause changes in opposite directions: mean shape index. Recall from lecture what patch shape index means. Do cannabis grows result in more regular or more irregular forest patches? Does timber harvest result in more regular or irregular forest patches?
Low shape index values indicate simple, more regular geometric shapes and higher values indicate complex, irregular shapes. Timber harvest result in more regular forest patches. Cannabis grow, because it has a lower effect of forest cover, results in less regular (more natural) patches. In this case, more regular forest patches are an indicator of grater human activity and landscape disturbance.
#mean.shape.index (griow vs timber)
diff.metrics$mean.shape.index ## [1] 0.001855851 -0.002158386We can see that cannabis grows cause a loss of 28,301,704 square meters of patch core area, which is much less than the 110,389,597 square meters of core area lost to timber harvest. However, let’s not forget that the total areas of cannabis grows and timber harvest are very different. So, let’s see what happens if we scale the core area lost relative to the total size of cannabis grows and timber harvest areas.
writeLines("Change in forest patch core area / total area of cannabis grows:")## Change in forest patch core area / total area of cannabis grows:diff.metrics$total.core.area[1] / ClassStat(grows, bkgd=0, latlon=T)$total.area## [1] -3.885524writeLines("Change in forest patch core area / total area of timber harvest:")## Change in forest patch core area / total area of timber harvest:diff.metrics$total.core.area[2] / ClassStat(timber, bkg=0, latlon=T)$total.area## [1] -2.146458Question 7 (1.5 pt)
What does this tell us about the placement of cannabis grows and areas of timber harvest in forest patches? i.e. Which is more likely to be found in the center of a forest patch? >>This is an example of the importance of scale. Even though there was a greater total patch core area lost due to timber harvest, if the core area lost is scaled to the total size of cannabis grows and timber harvest areas, it is easy to see that there is a grater change in forest patch core area caused by cannabis grow. Cannabis grow is more likely to be found in the center of a forest patch.
The End
That’s all for this week. Don’t forget to save your work. And when you’re done, knit it and submit it.