Skip to content

Commit 1281e09

Browse files
committed
Episode 11 updated by replacing calls to raster and rgdal packages with calls to terra package. #368 #363
1 parent 231af26 commit 1281e09

File tree

1 file changed

+83
-74
lines changed

1 file changed

+83
-74
lines changed

episodes/11-vector-raster-integration.Rmd

Lines changed: 83 additions & 74 deletions
Original file line numberDiff line numberDiff line change
@@ -24,22 +24,24 @@ source("setup.R")
2424

2525
```{r load-libraries, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
2626
library(sf)
27-
library(raster)
28-
library(rgdal)
27+
library(terra)
2928
library(ggplot2)
3029
library(dplyr)
3130
```
3231

3332
```{r load-data, echo=FALSE, results="hide"}
3433
# Learners will have this data loaded from earlier episodes
3534
# shapefiles
36-
point_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
37-
lines_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")
38-
aoi_boundary_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")
35+
point_HARV <-
36+
st_read("data/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
37+
lines_HARV <-
38+
st_read("data/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")
39+
aoi_boundary_HARV <-
40+
st_read("data/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")
3941
4042
# CHM
4143
CHM_HARV <-
42-
raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")
44+
rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")
4345
4446
CHM_HARV_df <- as.data.frame(CHM_HARV, xy = TRUE)
4547
@@ -56,8 +58,9 @@ plot_locations_sp_HARV <- st_as_sf(plot_locations_HARV,
5658

5759
## Things You'll Need To Complete This Episode
5860

59-
See the [lesson homepage](.) for detailed information about the software,
60-
data, and other prerequisites you will need to work through the examples in this episode.
61+
See the [lesson homepage](.) for detailed information about the software, data,
62+
and other prerequisites you will need to work through the examples in this
63+
episode.
6164

6265

6366
::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -71,7 +74,7 @@ points.
7174

7275
We often work with spatial layers that have different spatial extents. The
7376
spatial extent of a shapefile or R spatial object represents the geographic
74-
"edge" or location that is the furthest north, south east and west. Thus is
77+
"edge" or location that is the furthest north, south east and west. Thus it
7578
represents the overall geographic coverage of the spatial object.
7679

7780
![](fig/dc-spatial-vector/spatial_extent.png){alt='Extent illustration'} Image Source: National
@@ -92,7 +95,8 @@ CHM_HARV_sp <- st_as_sf(CHM_HARV_df, coords = c("x", "y"), crs = utm18nCRS)
9295
# approximate the boundary box with a random sample of raster points
9396
CHM_rand_sample <- sample_n(CHM_HARV_sp, 10000)
9497
lines_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")
95-
plots_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/PlotLocations_HARV.shp")
98+
plots_HARV <-
99+
st_read("data/NEON-DS-Site-Layout-Files/HARV/PlotLocations_HARV.shp")
96100
```
97101

98102
```{r compare-data-extents, echo=FALSE}
@@ -113,10 +117,10 @@ ggplot() +
113117

114118
Frequent use cases of cropping a raster file include reducing file size and
115119
creating maps. Sometimes we have a raster file that is much larger than our
116-
study area or area of interest. It is often more efficient to crop the
117-
raster to the extent of our study area to reduce file sizes as we process
118-
our data. Cropping a raster can also be useful when creating pretty maps so
119-
that the raster layer matches the extent of the desired vector layers.
120+
study area or area of interest. It is often more efficient to crop the raster
121+
to the extent of our study area to reduce file sizes as we process our data.
122+
Cropping a raster can also be useful when creating pretty maps so that the
123+
raster layer matches the extent of the desired vector layers.
120124

121125
## Crop a Raster Using Vector Extent
122126

@@ -125,7 +129,10 @@ spatial object. To do this, we need to specify the raster to be cropped and the
125129
spatial object that will be used to crop the raster. R will use the `extent` of
126130
the spatial object as the cropping boundary.
127131

128-
To illustrate this, we will crop the Canopy Height Model (CHM) to only include the area of interest (AOI). Let's start by plotting the full extent of the CHM data and overlay where the AOI falls within it. The boundaries of the AOI will be colored blue, and we use `fill = NA` to make the area transparent.
132+
To illustrate this, we will crop the Canopy Height Model (CHM) to only include
133+
the area of interest (AOI). Let's start by plotting the full extent of the CHM
134+
data and overlay where the AOI falls within it. The boundaries of the AOI will
135+
be colored blue, and we use `fill = NA` to make the area transparent.
129136

130137
```{r crop-by-vector-extent}
131138
ggplot() +
@@ -144,12 +151,12 @@ that falls within the boundaries of the AOI.
144151
CHM_HARV_Cropped <- crop(x = CHM_HARV, y = aoi_boundary_HARV)
145152
```
146153

147-
Now we can plot the cropped CHM data, along with a boundary box showing the full
148-
CHM extent. However, remember, since this is raster data, we need to convert to
149-
a data frame in order to plot using `ggplot`. To get the boundary box from CHM,
150-
the `st_bbox()` will extract the 4 corners of the rectangle that encompass all
151-
the features contained in this object. The `st_as_sfc()` converts these 4
152-
coordinates into a polygon that we can plot:
154+
Now we can plot the cropped CHM data, along with a boundary box showing the
155+
full CHM extent. However, remember, since this is raster data, we need to
156+
convert to a data frame in order to plot using `ggplot`. To get the boundary
157+
box from CHM, the `st_bbox()` will extract the 4 corners of the rectangle that
158+
encompass all the features contained in this object. The `st_as_sfc()` converts
159+
these 4 coordinates into a polygon that we can plot:
153160

154161
```{r show-cropped-area}
155162
CHM_HARV_Cropped_df <- as.data.frame(CHM_HARV_Cropped, xy = TRUE)
@@ -186,9 +193,9 @@ st_bbox(aoi_boundary_HARV)
186193
st_bbox(plot_locations_sp_HARV)
187194
```
188195

189-
Our plot location extent is not the largest but is larger than the AOI Boundary.
190-
It would be nice to see our vegetation plot locations plotted on top of the
191-
Canopy Height Model information.
196+
Our plot location extent is not the largest but is larger than the AOI
197+
Boundary. It would be nice to see our vegetation plot locations plotted on top
198+
of the Canopy Height Model information.
192199

193200
::::::::::::::::::::::::::::::::::::::: challenge
194201

@@ -219,21 +226,22 @@ ggplot() +
219226
::::::::::::::::::::::::::::::::::::::::::::::::::
220227

221228
In the plot above, created in the challenge, all the vegetation plot locations
222-
(black dots) appear on the Canopy Height Model raster layer except for one. One is
223-
situated on the blank space to the left of the map. Why?
229+
(black dots) appear on the Canopy Height Model raster layer except for one. One
230+
is situated on the blank space to the left of the map. Why?
224231

225232
A modification of the first figure in this episode is below, showing the
226233
relative extents of all the spatial objects. Notice that the extent for our
227234
vegetation plot layer (black) extends further west than the extent of our CHM
228-
raster (bright green). The `crop()` function will make a raster extent smaller, it
229-
will not expand the extent in areas where there are no data. Thus, the extent of our
230-
vegetation plot layer will still extend further west than the extent of our
231-
(cropped) raster data (dark green).
235+
raster (bright green). The `crop()` function will make a raster extent smaller,
236+
it will not expand the extent in areas where there are no data. Thus, the
237+
extent of our vegetation plot layer will still extend further west than the
238+
extent of our (cropped) raster data (dark green).
232239

233240
```{r, echo=FALSE}
234241
# code not shown, demonstration only
235242
# create CHM_plots_HARVcrop as a shape file
236-
CHM_plots_HARVcrop_sp <- st_as_sf(CHM_plots_HARVcrop_df, coords = c("x", "y"), crs = utm18nCRS)
243+
CHM_plots_HARVcrop_sp <- st_as_sf(CHM_plots_HARVcrop_df, coords = c("x", "y"),
244+
crs = utm18nCRS)
237245
# approximate the boundary box with random sample of raster points
238246
CHM_plots_HARVcrop_sp_rand_sample = sample_n(CHM_plots_HARVcrop_sp, 10000)
239247
```
@@ -250,7 +258,7 @@ will provide the `extent()` function our xmin, xmax, ymin, and ymax (in that
250258
order).
251259

252260
```{r}
253-
new_extent <- extent(732161.2, 732238.7, 4713249, 4713333)
261+
new_extent <- ext(732161.2, 732238.7, 4713249, 4713333)
254262
class(new_extent)
255263
```
256264

@@ -259,8 +267,8 @@ class(new_extent)
259267
## Data Tip
260268

261269
The extent can be created from a numeric vector (as shown above), a matrix, or
262-
a list. For more details see the `extent()` function help file
263-
(`?raster::extent`).
270+
a list. For more details see the `ext()` function help file
271+
(`?terra::ext`).
264272

265273

266274
::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -278,7 +286,8 @@ To plot this data using `ggplot()` we need to convert it to a dataframe.
278286
CHM_HARV_manual_cropped_df <- as.data.frame(CHM_HARV_manual_cropped, xy = TRUE)
279287
```
280288

281-
Now we can plot this cropped data. We will show the AOI boundary on the same plot for scale.
289+
Now we can plot this cropped data. We will show the AOI boundary on the same
290+
plot for scale.
282291

283292
```{r show-manual-crop-area}
284293
ggplot() +
@@ -292,38 +301,40 @@ ggplot() +
292301
## Extract Raster Pixels Values Using Vector Polygons
293302

294303
Often we want to extract values from a raster layer for particular locations -
295-
for example, plot locations that we are sampling on the ground. We can extract all pixel values within 20m of our x,y point of interest. These can then be summarized into some value of interest (e.g. mean, maximum, total).
304+
for example, plot locations that we are sampling on the ground. We can extract
305+
all pixel values within 20m of our x,y point of interest. These can then be
306+
summarized into some value of interest (e.g. mean, maximum, total).
296307

297-
![Extract raster information using a polygon boundary. From https://www.neonscience.org/sites/default/files/images/spatialData/BufferSquare.png](fig//BufferSquare.png)
308+
![Extract raster information using a polygon boundary. From https://www.neonscience.org/sites/default/files/images/spatialData/BufferSquare.png](fig//BufferSquare.png){alt = "Extract raster information using a polygon boundary. From https://www.neonscience.org/sites/default/files/images/spatialData/BufferSquare.png"}
298309

299310
To do this in R, we use the `extract()` function. The `extract()` function
300311
requires:
301312

302313
- The raster that we wish to extract values from,
303314
- The vector layer containing the polygons that we wish to use as a boundary or
304315
boundaries,
305-
- we can tell it to store the output values in a data frame using
306-
`df = TRUE`. (This is optional, the default is to return a list, NOT a data frame.) .
316+
- we can tell it to store the output values in a data frame using
317+
`raw = FALSE` (this is optional).
307318

308319
We will begin by extracting all canopy height pixel values located within our
309-
`aoi_boundary_HARV` polygon which surrounds the tower located at the NEON Harvard
310-
Forest field site.
320+
`aoi_boundary_HARV` polygon which surrounds the tower located at the NEON
321+
Harvard Forest field site.
311322

312323
```{r extract-from-raster}
313324
tree_height <- extract(x = CHM_HARV, y = aoi_boundary_HARV, df = TRUE)
314325
315326
str(tree_height)
316327
```
317328

318-
When we use the `extract()` function, R extracts the value for each pixel located
319-
within the boundary of the polygon being used to perform the extraction - in
320-
this case the `aoi_boundary_HARV` object (a single polygon). Here, the
329+
When we use the `extract()` function, R extracts the value for each pixel
330+
located within the boundary of the polygon being used to perform the extraction
331+
- in this case the `aoi_boundary_HARV` object (a single polygon). Here, the
321332
function extracted values from 18,450 pixels.
322333

323334
We can create a histogram of tree height values within the boundary to better
324335
understand the structure or height distribution of trees at our site. We will
325-
use the column `layer` from our data frame as our x values, as this column
326-
represents the tree heights for each pixel.
336+
use the column `HARV_chmCrop` from our data frame as our x values, as this
337+
column represents the tree heights for each pixel.
327338

328339
```{r view-extract-histogram}
329340
ggplot() +
@@ -333,10 +344,9 @@ ggplot() +
333344
ylab("Frequency of Pixels")
334345
```
335346

336-
We can also use the
337-
`summary()` function to view descriptive statistics including min, max, and mean
338-
height values. These values help us better understand vegetation at our field
339-
site.
347+
We can also use the `summary()` function to view descriptive statistics
348+
including min, max, and mean height values. These values help us better
349+
understand vegetation at our field site.
340350

341351
```{r}
342352
summary(tree_height$HARV_chmCrop)
@@ -345,12 +355,12 @@ summary(tree_height$HARV_chmCrop)
345355
## Summarize Extracted Raster Values
346356

347357
We often want to extract summary values from a raster. We can tell R the type
348-
of summary statistic we are interested in using the `fun =` argument. Let's extract
349-
a mean height value for our AOI. Because we are extracting only a single number, we will
350-
not use the `df = TRUE` argument.
358+
of summary statistic we are interested in using the `fun =` argument. Let's
359+
extract a mean height value for our AOI.
351360

352361
```{r summarize-extract}
353-
mean_tree_height_AOI <- extract(x = CHM_HARV, y = aoi_boundary_HARV, fun = mean)
362+
mean_tree_height_AOI <- extract(x = CHM_HARV, y = aoi_boundary_HARV,
363+
fun = mean)
354364
355365
mean_tree_height_AOI
356366
```
@@ -361,23 +371,23 @@ canopy height model is 22.43 meters.
361371
## Extract Data using x,y Locations
362372

363373
We can also extract pixel values from a raster by defining a buffer or area
364-
surrounding individual point locations using the `extract()` function. To do this
365-
we define the summary argument (`fun = mean`) and the buffer distance (`buffer = 20`)
366-
which represents the radius of a circular region around each point. By default, the units of the
367-
buffer are the same units as the data's CRS. All pixels that are touched by the buffer region are included in the extract.
374+
surrounding individual point locations using the `st_buffer()` function. To do
375+
this we define the summary argument (`fun = mean`) and the buffer distance
376+
(`dist = 20`) which represents the radius of a circular region around each
377+
point. By default, the units of the buffer are the same units as the data's
378+
CRS. All pixels that are touched by the buffer region are included in the
379+
extract.
368380

369-
![Extract raster information using a buffer region. From: https://www.neonscience.org/sites/default/files/images/spatialData/BufferCircular.png](fig/BufferCircular.png)
381+
![Extract raster information using a buffer region. From: https://www.neonscience.org/sites/default/files/images/spatialData/BufferCircular.png](fig/BufferCircular.png){alt = "Extract raster information using a buffer region. From: https://www.neonscience.org/sites/default/files/images/spatialData/BufferCircular.png"}
370382

371383
Source: National Ecological Observatory Network (NEON).
372384

373-
Let's put this into practice by figuring out the mean tree height in the
374-
20m around the tower location (`point_HARV`). Because we are extracting only a single number, we
375-
will not use the `df = TRUE` argument.
385+
Let's put this into practice by figuring out the mean tree height in the 20m
386+
around the tower location (`point_HARV`).
376387

377388
```{r extract-point-to-buffer}
378389
mean_tree_height_tower <- extract(x = CHM_HARV,
379-
y = point_HARV,
380-
buffer = 20,
390+
y = st_buffer(point_HARV, dist = 20),
381391
fun = mean)
382392
383393
mean_tree_height_tower
@@ -387,11 +397,10 @@ mean_tree_height_tower
387397

388398
## Challenge: Extract Raster Height Values For Plot Locations
389399

390-
1) Use the plot locations object (`plot_locations_sp_HARV`)
391-
to extract an average tree height for the
392-
area within 20m of each vegetation plot location in the study area. Because there are
393-
multiple plot locations, there will be multiple averages returned, so the `df = TRUE`
394-
argument should be used.
400+
1) Use the plot locations object (`plot_locations_sp_HARV`) to extract an
401+
average tree height for the area within 20m of each vegetation plot location
402+
in the study area. Because there are multiple plot locations, there will be
403+
multiple averages returned.
395404

396405
2) Create a plot showing the mean tree height of each area.
397406

@@ -402,10 +411,9 @@ mean_tree_height_tower
402411
```{r hist-tree-height-veg-plot}
403412
# extract data at each plot location
404413
mean_tree_height_plots_HARV <- extract(x = CHM_HARV,
405-
y = plot_locations_sp_HARV,
406-
buffer = 20,
407-
fun = mean,
408-
df = TRUE)
414+
y = st_buffer(plot_locations_sp_HARV,
415+
dist = 20),
416+
fun = mean)
409417
410418
# view data
411419
mean_tree_height_plots_HARV
@@ -427,7 +435,8 @@ ggplot(data = mean_tree_height_plots_HARV, aes(ID, HARV_chmCrop)) +
427435
:::::::::::::::::::::::::::::::::::::::: keypoints
428436

429437
- Use the `crop()` function to crop a raster object.
430-
- Use the `extract()` function to extract pixels from a raster object that fall within a particular extent boundary.
438+
- Use the `extract()` function to extract pixels from a raster object that fall
439+
within a particular extent boundary.
431440
- Use the `extent()` function to define an extent.
432441

433442
::::::::::::::::::::::::::::::::::::::::::::::::::

0 commit comments

Comments
 (0)