Skip to content

Commit 8744982

Browse files
authored
Merge pull request #396 from albhasan/m2022q1_e11
Episode 11 updated by replacing calls to raster and rgdal packages
2 parents fa1d1ca + 2cf06e7 commit 8744982

File tree

1 file changed

+93
-82
lines changed

1 file changed

+93
-82
lines changed

episodes/11-vector-raster-integration.Rmd

Lines changed: 93 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -18,28 +18,31 @@ source("setup.R")
1818

1919
:::::::::::::::::::::::::::::::::::::::: questions
2020

21-
- How can I crop raster objects to vector objects, and extract the summary of raster pixels?
21+
- How can I crop raster objects to vector objects, and extract the summary of
22+
raster pixels?
2223

2324
::::::::::::::::::::::::::::::::::::::::::::::::::
2425

2526
```{r load-libraries, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
2627
library(sf)
27-
library(raster)
28-
library(rgdal)
28+
library(terra)
2929
library(ggplot2)
3030
library(dplyr)
3131
```
3232

3333
```{r load-data, echo=FALSE, results="hide"}
3434
# Learners will have this data loaded from earlier episodes
3535
# 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")
36+
point_HARV <-
37+
st_read("data/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
38+
lines_HARV <-
39+
st_read("data/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")
40+
aoi_boundary_HARV <-
41+
st_read("data/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")
3942
4043
# CHM
4144
CHM_HARV <-
42-
raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")
45+
rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")
4346
4447
CHM_HARV_df <- as.data.frame(CHM_HARV, xy = TRUE)
4548
@@ -56,8 +59,9 @@ plot_locations_sp_HARV <- st_as_sf(plot_locations_HARV,
5659

5760
## Things You'll Need To Complete This Episode
5861

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.
62+
See the [lesson homepage](.) for detailed information about the software, data,
63+
and other prerequisites you will need to work through the examples in this
64+
episode.
6165

6266

6367
::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -71,7 +75,7 @@ points.
7175

7276
We often work with spatial layers that have different spatial extents. The
7377
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
78+
"edge" or location that is the furthest north, south east and west. Thus it
7579
represents the overall geographic coverage of the spatial object.
7680

7781
![](fig/dc-spatial-vector/spatial_extent.png){alt='Extent illustration'} Image Source: National
@@ -92,7 +96,8 @@ CHM_HARV_sp <- st_as_sf(CHM_HARV_df, coords = c("x", "y"), crs = utm18nCRS)
9296
# approximate the boundary box with a random sample of raster points
9397
CHM_rand_sample <- sample_n(CHM_HARV_sp, 10000)
9498
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")
99+
plots_HARV <-
100+
st_read("data/NEON-DS-Site-Layout-Files/HARV/PlotLocations_HARV.shp")
96101
```
97102

98103
```{r compare-data-extents, echo=FALSE}
@@ -113,10 +118,10 @@ ggplot() +
113118

114119
Frequent use cases of cropping a raster file include reducing file size and
115120
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.
121+
study area or area of interest. It is often more efficient to crop the raster
122+
to the extent of our study area to reduce file sizes as we process our data.
123+
Cropping a raster can also be useful when creating pretty maps so that the
124+
raster layer matches the extent of the desired vector layers.
120125

121126
## Crop a Raster Using Vector Extent
122127

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

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

130138
```{r crop-by-vector-extent}
131139
ggplot() +
@@ -144,12 +152,12 @@ that falls within the boundaries of the AOI.
144152
CHM_HARV_Cropped <- crop(x = CHM_HARV, y = aoi_boundary_HARV)
145153
```
146154

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:
155+
Now we can plot the cropped CHM data, along with a boundary box showing the
156+
full CHM extent. However, remember, since this is raster data, we need to
157+
convert to a data frame in order to plot using `ggplot`. To get the boundary
158+
box from CHM, the `st_bbox()` will extract the 4 corners of the rectangle that
159+
encompass all the features contained in this object. The `st_as_sfc()` converts
160+
these 4 coordinates into a polygon that we can plot:
153161

154162
```{r show-cropped-area}
155163
CHM_HARV_Cropped_df <- as.data.frame(CHM_HARV_Cropped, xy = TRUE)
@@ -186,9 +194,9 @@ st_bbox(aoi_boundary_HARV)
186194
st_bbox(plot_locations_sp_HARV)
187195
```
188196

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.
197+
Our plot location extent is not the largest but is larger than the AOI
198+
Boundary. It would be nice to see our vegetation plot locations plotted on top
199+
of the Canopy Height Model information.
192200

193201
::::::::::::::::::::::::::::::::::::::: challenge
194202

@@ -208,7 +216,8 @@ CHM_plots_HARVcrop <- crop(x = CHM_HARV, y = plot_locations_sp_HARV)
208216
CHM_plots_HARVcrop_df <- as.data.frame(CHM_plots_HARVcrop, xy = TRUE)
209217
210218
ggplot() +
211-
geom_raster(data = CHM_plots_HARVcrop_df, aes(x = x, y = y, fill = HARV_chmCrop)) +
219+
geom_raster(data = CHM_plots_HARVcrop_df,
220+
aes(x = x, y = y, fill = HARV_chmCrop)) +
212221
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
213222
geom_sf(data = plot_locations_sp_HARV) +
214223
coord_sf()
@@ -219,21 +228,22 @@ ggplot() +
219228
::::::::::::::::::::::::::::::::::::::::::::::::::
220229

221230
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?
231+
(black dots) appear on the Canopy Height Model raster layer except for one. One
232+
is situated on the blank space to the left of the map. Why?
224233

225234
A modification of the first figure in this episode is below, showing the
226235
relative extents of all the spatial objects. Notice that the extent for our
227236
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).
237+
raster (bright green). The `crop()` function will make a raster extent smaller,
238+
it will not expand the extent in areas where there are no data. Thus, the
239+
extent of our vegetation plot layer will still extend further west than the
240+
extent of our (cropped) raster data (dark green).
232241

233242
```{r, echo=FALSE}
234243
# code not shown, demonstration only
235244
# 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)
245+
CHM_plots_HARVcrop_sp <- st_as_sf(CHM_plots_HARVcrop_df, coords = c("x", "y"),
246+
crs = utm18nCRS)
237247
# approximate the boundary box with random sample of raster points
238248
CHM_plots_HARVcrop_sp_rand_sample = sample_n(CHM_plots_HARVcrop_sp, 10000)
239249
```
@@ -244,13 +254,13 @@ CHM_plots_HARVcrop_sp_rand_sample = sample_n(CHM_plots_HARVcrop_sp, 10000)
244254
## Define an Extent
245255

246256
So far, we have used a shapefile to crop the extent of a raster dataset.
247-
Alternatively, we can also the `extent()` function to define an extent to be
257+
Alternatively, we can also the `ext()` function to define an extent to be
248258
used as a cropping boundary. This creates a new object of class extent. Here we
249-
will provide the `extent()` function our xmin, xmax, ymin, and ymax (in that
259+
will provide the `ext()` function our xmin, xmax, ymin, and ymax (in that
250260
order).
251261

252262
```{r}
253-
new_extent <- extent(732161.2, 732238.7, 4713249, 4713333)
263+
new_extent <- ext(732161.2, 732238.7, 4713249, 4713333)
254264
class(new_extent)
255265
```
256266

@@ -259,8 +269,8 @@ class(new_extent)
259269
## Data Tip
260270

261271
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`).
272+
a list. For more details see the `ext()` function help file
273+
(`?terra::ext`).
264274

265275

266276
::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -278,7 +288,8 @@ To plot this data using `ggplot()` we need to convert it to a dataframe.
278288
CHM_HARV_manual_cropped_df <- as.data.frame(CHM_HARV_manual_cropped, xy = TRUE)
279289
```
280290

281-
Now we can plot this cropped data. We will show the AOI boundary on the same plot for scale.
291+
Now we can plot this cropped data. We will show the AOI boundary on the same
292+
plot for scale.
282293

283294
```{r show-manual-crop-area}
284295
ggplot() +
@@ -292,38 +303,41 @@ ggplot() +
292303
## Extract Raster Pixels Values Using Vector Polygons
293304

294305
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).
306+
for example, plot locations that we are sampling on the ground. We can extract
307+
all pixel values within 20m of our x,y point of interest. These can then be
308+
summarized into some value of interest (e.g. mean, maximum, total).
296309

297-
![Extract raster information using a polygon boundary. From https://www.neonscience.org/sites/default/files/images/spatialData/BufferSquare.png](fig//BufferSquare.png)
310+
![](fig//BufferSquare.png){alt='Image shows raster information extraction using 20m polygon boundary.'}
311+
Image Source: National Ecological Observatory Network (NEON)
298312

299313
To do this in R, we use the `extract()` function. The `extract()` function
300314
requires:
301315

302316
- The raster that we wish to extract values from,
303317
- The vector layer containing the polygons that we wish to use as a boundary or
304318
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.) .
319+
- we can tell it to store the output values in a data frame using
320+
`raw = FALSE` (this is optional).
307321

308322
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.
323+
`aoi_boundary_HARV` polygon which surrounds the tower located at the NEON
324+
Harvard Forest field site.
311325

312326
```{r extract-from-raster}
313-
tree_height <- extract(x = CHM_HARV, y = aoi_boundary_HARV, df = TRUE)
327+
tree_height <- extract(x = CHM_HARV, y = aoi_boundary_HARV, raw = FALSE)
314328
315329
str(tree_height)
316330
```
317331

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
332+
When we use the `extract()` function, R extracts the value for each pixel
333+
located within the boundary of the polygon being used to perform the extraction
334+
- in this case the `aoi_boundary_HARV` object (a single polygon). Here, the
321335
function extracted values from 18,450 pixels.
322336

323337
We can create a histogram of tree height values within the boundary to better
324338
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.
339+
use the column `HARV_chmCrop` from our data frame as our x values, as this
340+
column represents the tree heights for each pixel.
327341

328342
```{r view-extract-histogram}
329343
ggplot() +
@@ -333,10 +347,9 @@ ggplot() +
333347
ylab("Frequency of Pixels")
334348
```
335349

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.
350+
We can also use the `summary()` function to view descriptive statistics
351+
including min, max, and mean height values. These values help us better
352+
understand vegetation at our field site.
340353

341354
```{r}
342355
summary(tree_height$HARV_chmCrop)
@@ -345,12 +358,12 @@ summary(tree_height$HARV_chmCrop)
345358
## Summarize Extracted Raster Values
346359

347360
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.
361+
of summary statistic we are interested in using the `fun =` argument. Let's
362+
extract a mean height value for our AOI.
351363

352364
```{r summarize-extract}
353-
mean_tree_height_AOI <- extract(x = CHM_HARV, y = aoi_boundary_HARV, fun = mean)
365+
mean_tree_height_AOI <- extract(x = CHM_HARV, y = aoi_boundary_HARV,
366+
fun = mean)
354367
355368
mean_tree_height_AOI
356369
```
@@ -361,23 +374,22 @@ canopy height model is 22.43 meters.
361374
## Extract Data using x,y Locations
362375

363376
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.
368-
369-
![Extract raster information using a buffer region. From: https://www.neonscience.org/sites/default/files/images/spatialData/BufferCircular.png](fig/BufferCircular.png)
377+
surrounding individual point locations using the `st_buffer()` function. To do
378+
this we define the summary argument (`fun = mean`) and the buffer distance
379+
(`dist = 20`) which represents the radius of a circular region around each
380+
point. By default, the units of the buffer are the same units as the data's
381+
CRS. All pixels that are touched by the buffer region are included in the
382+
extract.
370383

371-
Source: National Ecological Observatory Network (NEON).
384+
![](fig/BufferCircular.png){alt='Image shows raster information extraction using 20m buffer region.'}
385+
Image Source: National Ecological Observatory Network (NEON)
372386

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.
387+
Let's put this into practice by figuring out the mean tree height in the 20m
388+
around the tower location (`point_HARV`).
376389

377390
```{r extract-point-to-buffer}
378391
mean_tree_height_tower <- extract(x = CHM_HARV,
379-
y = point_HARV,
380-
buffer = 20,
392+
y = st_buffer(point_HARV, dist = 20),
381393
fun = mean)
382394
383395
mean_tree_height_tower
@@ -387,11 +399,10 @@ mean_tree_height_tower
387399

388400
## Challenge: Extract Raster Height Values For Plot Locations
389401

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.
402+
1) Use the plot locations object (`plot_locations_sp_HARV`) to extract an
403+
average tree height for the area within 20m of each vegetation plot location
404+
in the study area. Because there are multiple plot locations, there will be
405+
multiple averages returned.
395406

396407
2) Create a plot showing the mean tree height of each area.
397408

@@ -402,10 +413,9 @@ mean_tree_height_tower
402413
```{r hist-tree-height-veg-plot}
403414
# extract data at each plot location
404415
mean_tree_height_plots_HARV <- extract(x = CHM_HARV,
405-
y = plot_locations_sp_HARV,
406-
buffer = 20,
407-
fun = mean,
408-
df = TRUE)
416+
y = st_buffer(plot_locations_sp_HARV,
417+
dist = 20),
418+
fun = mean)
409419
410420
# view data
411421
mean_tree_height_plots_HARV
@@ -427,8 +437,9 @@ ggplot(data = mean_tree_height_plots_HARV, aes(ID, HARV_chmCrop)) +
427437
:::::::::::::::::::::::::::::::::::::::: keypoints
428438

429439
- 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.
431-
- Use the `extent()` function to define an extent.
440+
- Use the `extract()` function to extract pixels from a raster object that fall
441+
within a particular extent boundary.
442+
- Use the `ext()` function to define an extent.
432443

433444
::::::::::::::::::::::::::::::::::::::::::::::::::
434445

0 commit comments

Comments
 (0)