You are probably well aware that, unfortunately, the Brazilian Amazon has been experiencing rapid deforestation for several decades, with the exception of a brief pause during stepped up monitoring and enforcement during roughly 2005 through 2012, corresponding to reforms undertaken under the first Da Silva administration.
In this tutorial, we will getting some practice using some geoprocessing tools in R in the context of analyzing the drivers of Amazonian forest loss. It may be helpful to refer back to the slides on geoprocessing to refresh your memory about the different operations we will be conducting.
Big-picture connection: why overlay matters for global environmental issues
Overlay operations (clip, intersect, union, buffer, erase) are how we turn layers into claims — for example, where roads intersect intact forests, or how much forest loss occurs within Indigenous territories. This is exactly the kind of spatial reasoning used in research that links environmental outcomes to governance, land rights, and resource development.
Indigenous lands & intact forest landscapes: Week 4’s analysis echoes work showing that a substantial share of intact forests overlaps Indigenous Peoples’ lands, and that these lands can be crucial for conserving forest integrity.
Resource development & overlapping risks: A different global overlay example maps thousands of mineral projects and quantifies how often they occur on or near lands of Indigenous and other land-connected peoples, alongside water/food/conflict risk indicators.
The Objectives
These are the main operations we will be working with in this tutorial:
Reading vector layers into R
Making basic plots with vector layers
Reprojecting vector layers
Simplifying vector geometries
Clipping
Unioning
Buffering
Erasing
Computing areas
The Data
Here are the datasets we will be using, along with their sources:
mato_grosso.gpkg: Boundaries of Mato Grosso, Brazil, a primary states in the “Arc of Deforestation” in the Brazilian Amazon, sourced from the Global Administrative Areas dataset.
IFL_Brazil_2000.gpkg and IFL_Brazil_2020.gpkg: Intact forest landscapes (i.e. undisturbed primary forests) in Brazil. Sourced from the Intact Forest Landscapes Team.
As we learned in our introduction to geospatial data in R, the sf library is now the primary workhorse for manipulating vector data. We are also going to use a package called rmapshaper, which has a few geoprocessing tools that work a bit better than their counterparts in the sf package.
Let’s install and activate the packages we’ll be using for this tutorial:
require(tidyverse)
Loading required package: tidyverse
Warning: package 'tidyverse' was built under R version 4.5.2
Warning: package 'ggplot2' was built under R version 4.5.2
Warning: package 'tidyr' was built under R version 4.5.2
Warning: package 'readr' was built under R version 4.5.2
Warning: package 'purrr' was built under R version 4.5.2
Warning: package 'forcats' was built under R version 4.5.2
Warning: package 'lubridate' was built under R version 4.5.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
require(sf)
Loading required package: sf
Warning: package 'sf' was built under R version 4.5.2
Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
require(rmapshaper)
Loading required package: rmapshaper
Warning: package 'rmapshaper' was built under R version 4.5.2
Reading in Vector Data with the sf Package
Just like with the terra package, we have to tell R where to look for the files holding any vector data we want to use. As before, it is easiest to do this by locating our code in the same folder as all of our data, setting our working directory accordingly, and then pointing R to the datasets we want to read in.
First, we can set the working directory. Remember that you can always use Session -> Set Working Directory -> To Source File Location to quickly find the path to the folder where your code and data are hanging out:
setwd("G:/My Drive/Classes/GIS for Global Environment/Labs/R_Tutorials_/Overlay")
Now let’s just read in one vector layer and plot it using the sf package. To read in vector data in the sf package, we use the st_read() function:
mato <-st_read("mato_grosso.gpkg")
Reading layer `arc_of_deforestation__gadm41_bra_1' from data source
`G:\My Drive\Classes\GIS for Global Environment\Labs\R_Tutorials_\Overlay\mato_grosso.gpkg'
using driver `GPKG'
Simple feature collection with 1 feature and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -61.6333 ymin: -18.0416 xmax: -50.2248 ymax: -7.349
Geodetic CRS: WGS 84
plot(mato)
Warning: plotting the first 9 out of 12 attributes; use max.plot = 12 to plot
all
Warning in min.default(structure(numeric(0), class = c("POSIXct", "POSIXt": no
non-missing arguments to min; returning Inf
Warning in max.default(structure(numeric(0), class = c("POSIXct", "POSIXt": no
non-missing arguments to max; returning -Inf
Warning in min.default(structure(numeric(0), class = c("POSIXct", "POSIXt": no
non-missing arguments to min; returning Inf
Warning in max.default(structure(numeric(0), class = c("POSIXct", "POSIXt": no
non-missing arguments to max; returning -Inf
Warning in min.default(structure(numeric(0), class = c("POSIXct", "POSIXt": no
non-missing arguments to min; returning Inf
Warning in max.default(structure(numeric(0), class = c("POSIXct", "POSIXt": no
non-missing arguments to max; returning -Inf
As you may remember from last time, if we just use the basic plot() function with an sf object, it will plot up to the first nine variables in the layer.
Deliverable 1
Now that we know how to read in a layer, let’s try reading in the boundaries of the legal Amazon. Using the code above as an example, read in the Legal_Amazon.gpkg file and assign it to an object called leg_amazon. Then plot() it. Include both the resulting plot and the code you used to generate it in your tutorial report. The result should look similar to the map below.
Reading layer `limites_amazonia_legal_2022' from data source
`G:\My Drive\Classes\GIS for Global Environment\Labs\R_Tutorials_\Overlay\Legal_Amazon.gpkg'
using driver `GPKG'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -73.99045 ymin: -18.0416 xmax: -44 ymax: 5.271841
Geodetic CRS: SIRGAS 2000
Basic Plotting with the sf Package
Now that we have two layers in memory, we can make some basic maps. Later in the tutorial, we’ll look at how to make more polished maps using the ggplot2 and ggspatial packages, but for right now we’ll just use the basic tools, which are helpful for doing quick checks on what your code is doing as you work through a project.
We’ll start by plotting just the boundaries of Mato Grosso. Because we’ll be adding another layer to the plot, we want to make sure that sf doesn’t plot all the variables in the mato object. Here’s a way to do that:
plot(st_geometry(mato))
Now, if we want to plot the outline of the Legal Amazon on top of this, we can’t just use the same plot command, because it’ll just create a new plot of only the Legal Amazon:
plot(st_geometry(leg_amazon))
Remember that this will only work if you already read in the Legal_Amazon.gpkg file and assigned it to an object called leg_amazon.
So we need a way to add multiple layers to a same plot. Here’s how we do that:
By including the add=TRUE argument, we’re telling R that we want to plot the second layer on top of the first.
But now we have the problem that it’s a bit unclear what we’re seeing - for all we know, those black lines could just be a river or something. Fortunately, we have a way to set the color for different layers:
That’s better! You might also have noticed that one cool thing about doing this in RStudio is that it will actually highlight color words with the color that will show up, giving you a bit of a preview of what the plot will look like.
There are lots of ways to reference colors in R, including commonly used RGB and Hex codes. When starting out, though, I tend to like to just use the color words. You can find all the color words for Rat this link.
Deliverable 2
Read in the IFL_Brazil_2020.gpkg file and assign it to an object called ifl_2020. Then read in the National_Roads_Brasil_Amazon.gpkg file and assign it to an object called nat_roads. Then create a plot showing the Legal Amazon in beige, the 2020 Intact Forest Landscapes in dark green, and the national roads in blue (be patient - these files are a bit large and can take a moment to read into memory). Include both your code and the resulting map in your tutorial report. The result should look similar to the map below.
Reading layer `clipped' from data source
`G:\My Drive\Classes\GIS for Global Environment\Labs\R_Tutorials_\Overlay\IFL_Brazil_2020.gpkg'
using driver `GPKG'
Simple feature collection with 206 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.01847 ymin: -25.58273 xmax: -43.16131 ymax: 5.267225
Geodetic CRS: WGS 84
Reading layer `vias_nacional' from data source
`G:\My Drive\Classes\GIS for Global Environment\Labs\R_Tutorials_\Overlay\National_Roads_Brasil_Amazon.gpkg'
using driver `GPKG'
Simple feature collection with 41089 features and 17 fields
Geometry type: MULTILINESTRING
Dimension: XYZ
Bounding box: xmin: -79.54228 ymin: -20.31399 xmax: -43.97888 ymax: 10.01478
z_range: zmin: 0 zmax: 0
Geodetic CRS: SIRGAS 2000
Clipping in the `rmapshaper`` Package
You’ve probably noticed that the nat_roads layer is operating with a different definition of the extent of the Amazon than we are. It covers the whole Amazon, not just the part in Brazil.
For the rest of this tutorial, we’ll focus our analysis just on Mato Grosso. That means we’re going to need to clip some of our other layers to fit.
The function in rmapshaper for clipping is just ms_clip():
nat_roads_mato <-ms_clip(nat_roads, mato)
Error: target and clip do not have identical CRS.
Uh-oh! We’ve been here before! Even though sf will do on-the-fly reprojection to plot objects with different CRSes, it won’t let us do calculations when we have different CRS. Time to learn how to reproject objects in the sf package!
Reprojecting in the sf Package
So what CRSes are we working with? We can see the CRS of any sf object using the st_crs() function:
st_crs(nat_roads)
Coordinate Reference System:
User input: SIRGAS 2000
wkt:
GEOGCRS["SIRGAS 2000",
DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
BBOX[-59.87,-122.19,32.72,-25.28]],
ID["EPSG",4674]]
st_crs(mato)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
The nat_roads object is in a datum called SIRGAS 2000 and is not projected (see the BBOX part of the output above? It’s clearly in decimal degrees). That means we need to reproject it into an appropriate local projection. The same is true of the mato object, which is in our much more familiar WGS 1984 datum.
So what projection should we use?
One good place to look up appropriate projections is https://spatialreference.org/explorer.html, which presents a map that you can use to find projections that would be appropriate in a particular location. You can also search for specific locations at the EPSG.
For this project, we’ll use the Universal Trans-Mercator system. That’s a system that splits the Earth into separate zones with their local projections. Mato Grosso is in UTM Zone 22 South (or UTM Zone 22S).
To reproject sf layers, we use the st_transform() function. This function takes two main arguments: first, the layer that we want to reproject, and, second, the CRS we want to reproject to. There are lots of ways to designate the CRS we want, but one way is to use the EPSG number, for the projection, which you can find at either of the two websites above.
Here’s how to reproject these objects into UTM Zone 22S:
Now should be able to perform our clip (be patient; this can take a bit):
nat_roads_mato <-ms_clip(nat_roads, mato)
Then we can have a look at our plots again:
Now let’s move on to some of our analysis. As we discussed at the outset, we want to look at how the loss of intact forest in Mato Grosso is affected by roads and indigenous territories. That means we need to compare the intact forest landscapes as of 2000 with those in 2020.
Deliverable 3
Read in the IFL_Brazil_2000.gpkg file and assign it to an object called ifl_2000. Plot it in blue Then plot the 2020 on top of it in dark green. To make the plot clearer, you can add in an argument, border = NA to the plot() function to keep R from drawing the borders. Include both your code and the resulting map in your tutorial report. Discuss what you see. Where were the losses of intact forest landscape in the Legal Amazon between 2000 and 2020?
Reading layer `ifl_brazil_2000' from data source
`G:\My Drive\Classes\GIS for Global Environment\Labs\R_Tutorials_\Overlay\IFL_Brazil_2000.gpkg'
using driver `GPKG'
Simple feature collection with 239 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.01847 ymin: -25.58273 xmax: -41.52281 ymax: 5.267225
Geodetic CRS: WGS 84
So it’s those blue areas that we’re really interested in, because those are the areas of intact forest landscape that were lost between 2000 and 2020.
Simplifying polygons in the rmapshaper Package
Before we continue with our analysis, we have to deal with a quick problem.
The intact forest landscape layers feature a common problem we encounter with data that have been produced by vectorizing raster layers - they have a very high vertex to area ratio. In other words, the polygons have really jagged edges.
This not only dramatically increases the file size and computational load when we’re working with these layers, it also can result in creating a bunch of invalid polygons when we use geoprocessing tools on them.
One way to mitigate these problems is to use simplification algorithms. ms_simplify() in rmapshaper, for example, uses a pretty sophisticated algorithm to remove vertices in order to shrink the detail in a layer while maintaining the overall shapes and neighborhood (or topological) relationships among the polygons. While using tools will sacrifice some accuracy, it will reduce our computational load and reduce the risk of invalid polygons later on.
ms_simplify() takes up to three arguments: first, the layer you want to simplify, then, the proportion of the original vertices you want to keep.
Because at the end of the day we want to be able to do some area calculations, we need to address the fact that the intact forest areas layers are not projected. Have a look:
st_crs(ifl_2000)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
st_crs(ifl_2020)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
As you can see, both layers are in WGS 1984.
We can use a neat trick to pull the CRS information from one object to set it as the target CRS for reprojecting another object. Here’s how:
A final data cleaning step that we should take has to do with the fact that, even simplified, the intact forest landscape layer edges are so jagged that comparing them is still likely to generate invalid polygons, just because even the boundary lines of the same forest area might not line up exactly in the two datasets.
One way to mitigate this problem is to slightly buffer the dataset we are using to do the erasing - in this case the 2020 intact forest landscapes layer. A very small buffer can help avoid polygon errors due to slight errors in polygon alignment between the two years. As with simplification, it comes with a cost in accuracy, as it will slightly bias the results downward, but this bias should be tiny at the scale of the entire Legal Amazon.
Buffering in the sf package is handled with the st_buffer() function. It takes two arguments: first, the layer whose features you want to buffer, and, second, the distance of the buffer, in the units used by the layer’s projection. In this case, the units are meters:
Now that these layers are simplified, we’re back on track and continue trying to identify the areas in the Arc of Deforestation that saw intact forest landscapes degraded between 2000 and 2020.
So how do we get them?
Well, you probably already have guessed from the title of this section that the answer is by erasing areas in the 2020 intact forest landscapes layer from the 2000 layer. We can infer that parts of the 2000 intact forest area that are not present in 2020 have been degraded.
Like it’s clipping function, rmapshaper’s erasing function is pretty simple. It’s just ms_erase(). The remove_slivers = TRUE argument tells R to automatically remove “sliver polygons”, which are tiny polygons that show up as a result of small differences in the coordinates of features in two different layers:
ifl_change <-ms_erase(ifl_2000, ifl_2020,remove_slivers =TRUE)plot(st_geometry(ifl_change), col ="red", border =NA)
Computing Areas in the sf Package
Once we have our forest change layer in the correct projection, we can compute its area. Once again, the function is pretty self-explanatory. It’s just st_area(). The only argument you need is the layer the area of whose polygons you want to compute:
Well, that isn’t the most helpful, is it? We really wanted to see the total area that was lost, not the area in each random polygon in the ifl_change object.
One way to get the entire area would be to sum this value up:
st_area(ifl_change) |>sum()
49796427578 [m^2]
Since square meters is a bit tough to wrap our heads around, maybe we should also convert it to square kilometers while we’re at it (there are 1 million square meters in a square kilometer):
st_area(ifl_change) |>sum()/1000000
49796.43 [m^2]
Okay, so that’s almost 50 thousand square kilometers of intact forest landscape in Mato Grosso degraded between 2000 and 2020 - depressing!
We could also do this calculation by unioning the polygons first and then computing the area. We perform unions with the st_union() function.
st_union(ifl_change) |>st_area()/1000000
49796.43 [m^2]
Note that R didn’t automatically convert the units designation when we did the division to convert to kilometers. That’s something to keep an eye on.
Deliverable 4
Write a line of code that would compute the area of intact forest landscapes in Mato Grosso in 2000. Then run the computation. Include both the code and the result in your tutorial report.
Putting it all Together: Intact Forest Loss along National Roads in the Arc of Deforestation
Finally, we’ve introduced almost all the tools that we need to conduct some analysis on the drivers of intact forest degradation in Mato Grosso. Like we’ve already discussed a lot in the class, geospatial analysis is really all about chaining tools together in clever ways. That’s what we’ll be doing here.
It has been well documented that road construction has been a major force driving Amazonian deforestation. Here, we’re going to estimate the share of intact forest landscape degradation in Mato Grosso between 2000 and 2020 within 5 kilometers of a national road.
So here are the steps we need to take to make this calculation: * Create a buffer of 10 kilometers around the national roads in Mato Grosso. * Clip the ifl_change layer with buffer layer. * Compute the area in the resulting clipped layer. * Compute the area of the clipped region as a percentage of the total degraded area between 2000 and 2020.
Since the emphasis in this section is understanding the flow of the process, I’m going to limit my commentary to simply labeling what we’re doing at each step.
Create Buffer around National Roads
nat_road_buff <-st_buffer(nat_roads_mato, dist =10000) |>#remember units are in meters!st_union() #only need one bufferplot(st_geometry(mato), col="beige")plot(st_geometry(nat_road_buff), col="blue", add=TRUE)
Clip Intact Forest Change Area with Buffer
nat_road_ch <-ms_clip(target = ifl_change,clip =st_as_sf(nat_road_buff), #had to add this in because rmapshaper wasn't recognizing the object typeremove_slivers =TRUE)plot(st_geometry(mato), col="beige")plot(st_geometry(nat_road_buff), col="blue", add=TRUE)plot(st_geometry(nat_road_ch), col="red", add=TRUE)
Compute Area of Clipped Region
st_area(nat_road_ch) |>sum()/1000000
1605.53 [m^2]
Compute Percentage of Total Degraded Area
#Compute degraded area around national roadsroad_area <-st_area(nat_road_ch) |>sum()#Compute total degraded areachange_area <-ms_clip(ifl_change, mato) |>st_area() |>sum()#Compute percentageroad_area/change_area*100
3.224186 [1]
So only about 3.2% of intact forest degradation in Mato Grosso between 2000 and 2020 took place within five kilometers of a national road. That seems surprising . . .
Final Tasks
Now it’s your turn to conduct some analysis.
Deliverable 5
Compute the share of total forest degradation that took place in indigenous territories in the Arc of Deforestation between 2000 and 2020. Then compute the share that took place within 5 kilometers of the boundary of indigenous territories (remember that you need to erase the area within indigenous territories in order to get just the buffer zone). Include all the code required to replicate your results.
Congratulations! You have completed this tutorial!
