Skip to content

Commit

Permalink
Merge pull request #396 from albhasan/m2022q1_e11
Browse files Browse the repository at this point in the history 
Episode 11 updated by replacing calls to raster and rgdal packages
  • Loading branch information
@drakeasberry
drakeasberry committed Mar 23, 2023
2 parents fa1d1ca + 2cf06e7 commit 8744982
Showing 1 changed file with 93 additions and 82 deletions.
175 changes: 93 additions & 82 deletions episodes/11-vector-raster-integration.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,28 +18,31 @@ source("setup.R")

:::::::::::::::::::::::::::::::::::::::: questions

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

::::::::::::::::::::::::::::::::::::::::::::::::::

```{r load-libraries, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
library(sf)
library(raster)
library(rgdal)
library(terra)
library(ggplot2)
library(dplyr)
```

```{r load-data, echo=FALSE, results="hide"}
# Learners will have this data loaded from earlier episodes
# shapefiles
point_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
lines_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")
aoi_boundary_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")
point_HARV <-
st_read("data/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
lines_HARV <-
st_read("data/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")
aoi_boundary_HARV <-
st_read("data/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")
# CHM
CHM_HARV <-
raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")
rast("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")
CHM_HARV_df <- as.data.frame(CHM_HARV, xy = TRUE)
Expand All @@ -56,8 +59,9 @@ plot_locations_sp_HARV <- st_as_sf(plot_locations_HARV,

## Things You'll Need To Complete This Episode

See the [lesson homepage](.) for detailed information about the software,
data, and other prerequisites you will need to work through the examples in this episode.
See the [lesson homepage](.) for detailed information about the software, data,
and other prerequisites you will need to work through the examples in this
episode.


::::::::::::::::::::::::::::::::::::::::::::::::::
Expand All @@ -71,7 +75,7 @@ points.

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

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

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

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

## Crop a Raster Using Vector Extent

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

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.
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.

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

Now we can plot the cropped CHM data, along with a boundary box showing the full
CHM extent. However, remember, since this is raster data, we need to convert to
a data frame in order to plot using `ggplot`. To get the boundary box from CHM,
the `st_bbox()` will extract the 4 corners of the rectangle that encompass all
the features contained in this object. The `st_as_sfc()` converts these 4
coordinates into a polygon that we can plot:
Now we can plot the cropped CHM data, along with a boundary box showing the
full CHM extent. However, remember, since this is raster data, we need to
convert to a data frame in order to plot using `ggplot`. To get the boundary
box from CHM, the `st_bbox()` will extract the 4 corners of the rectangle that
encompass all the features contained in this object. The `st_as_sfc()` converts
these 4 coordinates into a polygon that we can plot:

```{r show-cropped-area}
CHM_HARV_Cropped_df <- as.data.frame(CHM_HARV_Cropped, xy = TRUE)
Expand Down Expand Up @@ -186,9 +194,9 @@ st_bbox(aoi_boundary_HARV)
st_bbox(plot_locations_sp_HARV)
```

Our plot location extent is not the largest but is larger than the AOI Boundary.
It would be nice to see our vegetation plot locations plotted on top of the
Canopy Height Model information.
Our plot location extent is not the largest but is larger than the AOI
Boundary. It would be nice to see our vegetation plot locations plotted on top
of the Canopy Height Model information.

::::::::::::::::::::::::::::::::::::::: challenge

Expand All @@ -208,7 +216,8 @@ CHM_plots_HARVcrop <- crop(x = CHM_HARV, y = plot_locations_sp_HARV)
CHM_plots_HARVcrop_df <- as.data.frame(CHM_plots_HARVcrop, xy = TRUE)
ggplot() +
geom_raster(data = CHM_plots_HARVcrop_df, aes(x = x, y = y, fill = HARV_chmCrop)) +
geom_raster(data = CHM_plots_HARVcrop_df,
aes(x = x, y = y, fill = HARV_chmCrop)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
geom_sf(data = plot_locations_sp_HARV) +
coord_sf()
Expand All @@ -219,21 +228,22 @@ ggplot() +
::::::::::::::::::::::::::::::::::::::::::::::::::

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

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

```{r, echo=FALSE}
# code not shown, demonstration only
# create CHM_plots_HARVcrop as a shape file
CHM_plots_HARVcrop_sp <- st_as_sf(CHM_plots_HARVcrop_df, coords = c("x", "y"), crs = utm18nCRS)
CHM_plots_HARVcrop_sp <- st_as_sf(CHM_plots_HARVcrop_df, coords = c("x", "y"),
crs = utm18nCRS)
# approximate the boundary box with random sample of raster points
CHM_plots_HARVcrop_sp_rand_sample = sample_n(CHM_plots_HARVcrop_sp, 10000)
```
Expand All @@ -244,13 +254,13 @@ CHM_plots_HARVcrop_sp_rand_sample = sample_n(CHM_plots_HARVcrop_sp, 10000)
## Define an Extent

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

```{r}
new_extent <- extent(732161.2, 732238.7, 4713249, 4713333)
new_extent <- ext(732161.2, 732238.7, 4713249, 4713333)
class(new_extent)
```

Expand All @@ -259,8 +269,8 @@ class(new_extent)
## Data Tip

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


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

Now we can plot this cropped data. We will show the AOI boundary on the same plot for scale.
Now we can plot this cropped data. We will show the AOI boundary on the same
plot for scale.

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

Often we want to extract values from a raster layer for particular locations -
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).
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).

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

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

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

We will begin by extracting all canopy height pixel values located within our
`aoi_boundary_HARV` polygon which surrounds the tower located at the NEON Harvard
Forest field site.
`aoi_boundary_HARV` polygon which surrounds the tower located at the NEON
Harvard Forest field site.

```{r extract-from-raster}
tree_height <- extract(x = CHM_HARV, y = aoi_boundary_HARV, df = TRUE)
tree_height <- extract(x = CHM_HARV, y = aoi_boundary_HARV, raw = FALSE)
str(tree_height)
```

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

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

```{r view-extract-histogram}
ggplot() +
Expand All @@ -333,10 +347,9 @@ ggplot() +
ylab("Frequency of Pixels")
```

We can also use the
`summary()` function to view descriptive statistics including min, max, and mean
height values. These values help us better understand vegetation at our field
site.
We can also use the `summary()` function to view descriptive statistics
including min, max, and mean height values. These values help us better
understand vegetation at our field site.

```{r}
summary(tree_height$HARV_chmCrop)
Expand All @@ -345,12 +358,12 @@ summary(tree_height$HARV_chmCrop)
## Summarize Extracted Raster Values

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

```{r summarize-extract}
mean_tree_height_AOI <- extract(x = CHM_HARV, y = aoi_boundary_HARV, fun = mean)
mean_tree_height_AOI <- extract(x = CHM_HARV, y = aoi_boundary_HARV,
fun = mean)
mean_tree_height_AOI
```
Expand All @@ -361,23 +374,22 @@ canopy height model is 22.43 meters.
## Extract Data using x,y Locations

We can also extract pixel values from a raster by defining a buffer or area
surrounding individual point locations using the `extract()` function. To do this
we define the summary argument (`fun = mean`) and the buffer distance (`buffer = 20`)
which represents the radius of a circular region around each point. By default, the units of the
buffer are the same units as the data's CRS. All pixels that are touched by the buffer region are included in the extract.

![Extract raster information using a buffer region. From: https://www.neonscience.org/sites/default/files/images/spatialData/BufferCircular.png](fig/BufferCircular.png)
surrounding individual point locations using the `st_buffer()` function. To do
this we define the summary argument (`fun = mean`) and the buffer distance
(`dist = 20`) which represents the radius of a circular region around each
point. By default, the units of the buffer are the same units as the data's
CRS. All pixels that are touched by the buffer region are included in the
extract.

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

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

```{r extract-point-to-buffer}
mean_tree_height_tower <- extract(x = CHM_HARV,
y = point_HARV,
buffer = 20,
y = st_buffer(point_HARV, dist = 20),
fun = mean)
mean_tree_height_tower
Expand All @@ -387,11 +399,10 @@ mean_tree_height_tower

## Challenge: Extract Raster Height Values For Plot Locations

1) Use the plot locations object (`plot_locations_sp_HARV`)
to extract an average tree height for the
area within 20m of each vegetation plot location in the study area. Because there are
multiple plot locations, there will be multiple averages returned, so the `df = TRUE`
argument should be used.
1) Use the plot locations object (`plot_locations_sp_HARV`) to extract an
average tree height for the area within 20m of each vegetation plot location
in the study area. Because there are multiple plot locations, there will be
multiple averages returned.

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

Expand All @@ -402,10 +413,9 @@ mean_tree_height_tower
```{r hist-tree-height-veg-plot}
# extract data at each plot location
mean_tree_height_plots_HARV <- extract(x = CHM_HARV,
y = plot_locations_sp_HARV,
buffer = 20,
fun = mean,
df = TRUE)
y = st_buffer(plot_locations_sp_HARV,
dist = 20),
fun = mean)
# view data
mean_tree_height_plots_HARV
Expand All @@ -427,8 +437,9 @@ ggplot(data = mean_tree_height_plots_HARV, aes(ID, HARV_chmCrop)) +
:::::::::::::::::::::::::::::::::::::::: keypoints

- Use the `crop()` function to crop a raster object.
- Use the `extract()` function to extract pixels from a raster object that fall within a particular extent boundary.
- Use the `extent()` function to define an extent.
- Use the `extract()` function to extract pixels from a raster object that fall
within a particular extent boundary.
- Use the `ext()` function to define an extent.

::::::::::::::::::::::::::::::::::::::::::::::::::

Expand Down

0 comments on commit 8744982

Please sign in to comment.