Chapter 4 - Elevation Data

2025-04-03

Chapter 4: Digital Elevation Data

Authors: Dimitrios Markou, Danielle Ethier

In Chapter 3, you processed vector and raster climate data, combined them with NatureCounts observations, and visualized them using plots and spatiotemporal maps. In this tutorial, you will process LiDAR-derived Digital Terrain Models, apply cropping and masking procedures, and extract elevation values to combine with NatureCounts data. Your focus will be on the NatureCounts and spatial data within La Mauricie National Park.

La Mauricie National Park is situated in the Laurentian mountains and covers 536 km2 within the Eastern Canadian Temperate-Boreal Forest transition ecoregion. The environment is characterized by mixed forests, lakes, rivers, and hills that range from from 150 m to over 500 m in elevation. The park provides suitable habitat for a variety of wildlife including at least 215 bird species.

Light Detection and Ranging (LiDAR) is an active remote sensing technology. It is performed using laser scanners that emit pulses of light and determine the position of target 3D objects by measuring the amount of time between pulses being emitted and received. It is a revolutionary technology that helps in the acquisition of extremely accurate digital elevation data over wide spatial and temporal scales. For more info and examples of LiDAR data resources specific to your province, see LiDAR Data Resources at the end of this Chapter.

4.0 Learning Objectives

By the end of Chapter 4 - Elevation Data, users will know how to:

• Distinguish between different types of digital elevation datasets: Elevation Models
• Process Digital Terrain Models (DTMs), including cropping and masking: Spatial Extents - Crop and Mask
• Extract elevation data (raster values) at observation sites (vector points): Map & Extract Digital Elevation Data
• Know how to access LiDAR-derived data products: LiDAR Data Resources and Manual Data Download

Load the necessary packages:

library(naturecounts)
library(sf)
library(stars)  # Another package for working with raster data
library(tidyverse)

This tutorial uses the following spatial data.

1. Places administered by Parks Canada - Boundary shapefiles

2. Quebec Breeding Bird Atlas (2010 - 2014) - NatureCounts bird observations

3. Forêt ouverte (lidar-derived products) - Digital Terrain Models (DTMs)

4.1 Data Setup

To read the National Park polygons into R, navigate to Places administered by Parks Canada, click “Go to Resource” to download. If you haven’t already, return to the Chapter 2 Learning Objectives to view (getwd()) and set (setwd()) your working directory and create a folder within this directory called data. This is where you will store all data used in this tutorial series.

Create a subdirectory called nationalparks to neatly store the parks spatial data.

dir.create("data/nationalparks", recursive = TRUE)

Move the files from your Downloads to the nationalparks subdirectory before applying the st_read() function.

# specify the path to your shapefile
national_parks <- st_read("path/to/your/shapefile.shp")

Filter the national_parks dataset for La Mauricie National Park.

View(national_parks) # to find the correct object ID

mauricie_boundary <- national_parks %>%
  filter(OBJECTID == "21")

# Drop the Z-dimension (3D component) to make it 2D
mauricie_boundary <- st_zm(mauricie_boundary, drop = TRUE, what = "ZM")

OPTIONAL: To save the boundary as a shapefile to your working directory, use the st_write() function. The boundary files are used again in subsequent chapters. We have uploaded them to the Google Drive data folder for your convenience.

Create a subdirectory called boundary to neatly store the boundary files.

dir.create("data/boundary", recursive = TRUE)

To execute this code chunk, remove the #

# st_write(mauricie_boundary, "data/boundary/mauricie_boundary.shp") 

# st_write(mauricie_boundary, "data/boundary/mauricie_boundary.kml", driver="KML")

To assess the species distribution within the National Park, download data from the Quebec Breeding Bird Atlas (2010 - 2014) which is a 5 year census that collects data on the distribution and abundance of all species breeding in the province.

Don’t forget to replace testuser with your NatureCounts username. You will be prompted for your password.

Download NatureCounts data:

quebec_atlas <- nc_data_dl(
  collections = "QCATLAS2PC",  
  username = "testuser", info = "spatial_data_tutorial", timeout = 1000)

Read in the list of species names from NatureCounts which can be linked to the species id:

species_names <- search_species()

To create date and doy columns and ensure that the ObservationCount column is in the correct numeric format we can apply the format_dates() and mutate() functions. We will also filter the dataset to exclude rows with missing coordinates.

quebec_atlas <- quebec_atlas %>%
  format_dates() %>%  # create the date and doy columns 
  mutate(ObservationCount = as.numeric(ObservationCount)) %>%  # convert to numeric format
  filter(!is.na(longitude) & !is.na(latitude))  # remove rows with missing coordinates

To convert the NatureCounts data to a spatial object and transform its coordinate reference system (CRS) to match the National Park boundary we can use the st_as_sf() and st_transform() functions, respectively.

quebec_atlas_sf <- sf::st_as_sf(quebec_atlas,
                        coords = c("longitude", "latitude"), crs = 4326) # converts the quebec_atlas data to an sf object, crs = 4326 is WGS84 which corresponds to the NatureCounts data

mauricie_boundary <- st_transform(mauricie_boundary, crs = st_crs(quebec_atlas_sf)) # match the CRS between the polygon and point data

Clip the NatureCounts point data to the National Park boundary using st_intersection().

mauricie_birds_sf <- sf::st_intersection(quebec_atlas_sf, mauricie_boundary)

Append the species names to the clipped NatureCounts dataset based on species_id code.

mauricie_birds_sf <- mauricie_birds_sf %>%
  left_join(species_names, by = "species_id")

Tidyverse functions can help us summarize our data in a variety of ways. For example, we can determine the annual total count of birds for each year across all sites. To do this we use the group_by() function to group the observations by year, and summarise() to calculate and create the annual_count column.

mauricie_birds_summary <- mauricie_birds_sf %>%
  group_by(survey_year) %>%
  summarise(annual_count = sum(ObservationCount, na.rm = TRUE)) %>% # calculates the annual_count
  filter(!is.na(survey_year))  # remove rows with missing year

mauricie_birds_summary
#> Simple feature collection with 4 features and 2 fields
#> Geometry type: MULTIPOINT
#> Dimension:     XY
#> Bounding box:  xmin: -73.12685 ymin: 46.65015 xmax: -72.77968 ymax: 46.82687
#> Geodetic CRS:  WGS 84
#> # A tibble: 4 × 3
#>   survey_year annual_count                                                                                  geometry
#> *       <int>        <dbl>                                                                          <MULTIPOINT [°]>
#> 1        2011          909 ((-72.83548 46.76367), (-72.83824 46.76044), (-72.84852 46.75341), (-72.85742 46.74331...
#> 2        2012          578 ((-72.77968 46.74193), (-72.80576 46.75797), (-72.80941 46.75672), (-72.83613 46.77124...
#> 3        2013          517 ((-72.80942 46.75671), (-72.83614 46.77124), (-72.84853 46.75341), (-72.82786 46.74054...
#> 4        2014         1396 ((-72.77968 46.74193), (-72.80576 46.75797), (-72.80941 46.75672), (-72.83613 46.77124...

If you wanted to summarize total count for each species at each Atlas Block you could adjust the pipe using pivot_wider().

mauricie_species_summary <- mauricie_birds_sf %>%
  st_drop_geometry() %>%  # drop the geometry column
  group_by(english_name, Locality) %>%
  summarise(total_count = sum(ObservationCount, na.rm = TRUE)) %>% # calculates the total_count column
  pivot_wider(names_from = english_name, # populates the column names with each species common name
              values_from = total_count, # populates each cell with total_count
              values_fill = list(total_count = 0)) # missing values are zero-filled
#> `summarise()` has grouped output by 'english_name'. You can override using the `.groups` argument.
 
View(mauricie_species_summary)

Drop the geometry column and convert the filtered NatureCounts data back to a regular dataframe. Select key attributes to reduce the dataframe.

mauricie_birds_df <- mauricie_birds_sf %>%
  st_drop_geometry() %>%  # Drops the geometry column
  bind_cols(
    st_coordinates(mauricie_birds_sf) %>%  # Extract coordinates
      as.data.frame()  # Convert the coordinates to a data.frame
  ) %>%
  rename(longitude = X, latitude = Y)

mauricie_birds_df <- mauricie_birds_df %>% dplyr::select(record_id, species_id, english_name, SiteCode, survey_year, survey_month, survey_day, ObservationCount, longitude, latitude) # reduce the dataframe to the columns needed for your analysis

OPTIONAL: To save the NatureCounts data as a .csv to your disk, use the write.csv() function, specify the name of your .csv file, and use the row.names = FALSE argument to exclude row numbers from the output. The filtered NatureCounts dataset will be used again in subsequent chapters. We have uploaded this file to the Google Drive data folder for your convenience.

To execute this code chunk, remove the #

# write.csv(mauricie_birds_df, "data/mauricie_birds_df.csv", row.names = FALSE)

4.2 Elevation Models

Digital elevation datasets store topographic information like elevation or slope and are often used in landscape ecology. These datasets, i.e., Digital Elevation Models (DEMs), Digital Surface Models (DSMs), and Digital Terrain Models (DTMs) are derived through a variety of remote sensing and spatial interpolation techniques and all help describe land features.

Digital Elevation Model - represents the bare-Earth surface and excludes all terrain vector features (i.e., streams, breaklines, and ridges), and all ground objects (power lines, buildings, trees, and vegetation).

Digital Surface Model - represents the heights of the Earth’s surface and includes all natural and artificial features or ground objects.

Digital Terrain Model - represents the bare-Earth surface topography and includes all terrain vector features. It does not include natural or artificial ground objects. In other words, it is a DEM that is augmented by the presence of streams, breaklines, and ridges.

For quick-access to the elevation data, download the .tif files from their URL using the hyperlinks below. If you wish to gain experience in how-to access LiDAR data yourself, follow the steps at the end of this Chapter.

Click the following hyperlinks to download the LiDAR-derived elevation data for tiles 31I14NE, 31I15NO, 31I14SE, 31I15SO, 31I11NE, and 31I10NO, respectively and then move them to their own folder named dtm stored in your data subdirectory.

Create a subdirectory called dtm to neatly store your downloaded TIF files. Move them from your Downloads to the dtm folder.

dir.create("data/dtm", recursive = TRUE)

Set the path to your TIF file subdirectory.

Create a mosaic of the adjacent DTM rasters.

Reminder: You should store the DTM’s in their own folder after download. The code below will produce an error if there are other .tif files in your subdirectory (like the ones you downloaded in Chapter 3).

Here we use the stars package which allows us to create ‘proxy’ objects, meaning that the files aren’t loaded into memory until the absolutely must be. This allows us to work much more quickly with the files, especially for plotting.

# list all the TIFF files in your directory
dtm_files <- list.files(output_dir, pattern = "\\.tif$", full.names = TRUE) 

# Combine into a single stars object (use proxy to only load what we need)
dtm_mosaic <- st_mosaic(dtm_files)
dtm_mosaic <- read_stars(dtm_mosaic, proxy = TRUE)

Let’s check if the DTM and National Park boundary have the same crs by using the st_crs() function and equality operator (==) which will generate either TRUE or FALSE.

st_crs(dtm_mosaic) == st_crs(mauricie_boundary)
#> [1] FALSE

To reproject the spatial data with the same CRS, we can use the st_transform() function.

mauricie_boundary <- st_transform(mauricie_boundary, crs = st_crs(dtm_mosaic))

4.3 Spatial Extents: Crop & Mask

• Cropping reduces the extent of a raster to the extent of another raster or vector.
• Masking assigns NA values to cells of a raster not covered by a vector extent.

The stars package function st_crop() performs both of these steps in one.

Here, we crop the extent of the raster to the National Park Boundary.

dtm_crop <- st_crop(dtm_mosaic, mauricie_boundary)

OPTIONAL: To write the masked raster to your disk, you can use the write_stars() function from stars (this will take a while as it is a large file). This raster will be used in Chapter 7: Summary Tools. We have uploaded this file to the Google Drive data folder for your convenience.

To execute this code chunk, remove the #

# output_dir <- "path/to/your/subdirectory/data/dtm_mask.tif" # include the masked raster name and '.tif' argument 
# write_stars(dtm_crop, output_dir)

Note: Make sure to include the ‘.tif’ argument when specifying the output_dir and name of your masked raster file. This specifies the file type while using write_stars() which will, otherwise, produce an error.

4.4 Map and Extract Digital Elevation Data

Visualize the regional and masked DTM’s with a two-panel plot.

# Plot dtm_mosaic
plot(dtm_mosaic, main = "regional DTM", 
     col = viridisLite::viridis(n = 100), reset = FALSE) # reset = FALSE allows us to stack plots
#> downsample set to 37

# Overlay the National Park boundary on the first plot
plot(st_as_sfc(mauricie_boundary), add = TRUE, border = "red", lwd = 2)

# Plot dtm_mask
plot(dtm_crop, main = "masked DTM", 
     col = viridisLite::viridis(n = 100), reset = FALSE)
#> downsample set to 28

# Overlay the National Park boundary on the second plot
#plot(st_as_sfc(mauricie_boundary), add = TRUE, border = "red", lwd = 2, col = NA, reset = FALSE)

Summarize the NatureCounts data for mapping. Here we will determine the number of species observed at each site.

# Group by SiteCode and summarize mean annual count across all species
mauricie_site_summary <- mauricie_birds_sf %>%
  group_by(SiteCode, geometry) %>% 
  summarise(n_species = n_distinct(english_name)) # count the number of unique species
#> `summarise()` has grouped output by 'SiteCode'. You can override using the `.groups` argument.

Now we can map the NatureCounts summary output with the DTM data.

# Match the CRS
mauricie_site_summary <- st_transform(mauricie_site_summary, crs = st_crs(dtm_crop))

#scale the n_species column to range from 0-5 for plotting
mauricie_site_summary$n_species_scale <- scales::rescale(mauricie_site_summary$n_species, to = c(0, 2))

# Plot the dtm_mask raster
plot(dtm_crop, 
     main = "Elevation & Species Diversity Summarized across La Mauricie",
     col = viridisLite::viridis(n = 100), reset = FALSE)
#> downsample set to 28

# Overlay mauricie_site_summary points
plot(st_geometry(mauricie_site_summary), 
     pch = 19,        # Solid circle
     col = "black",   # Point color
     cex = mauricie_site_summary$n_species_scale, # Point size based on total
     reset = FALSE, add = TRUE)      

# Define legend values for the point size scale
legend_values <- c(min(mauricie_site_summary$n_species, na.rm = TRUE), 
                   median(mauricie_site_summary$n_species, na.rm = TRUE), 
                   max(mauricie_site_summary$n_species, na.rm = TRUE))
legend_sizes <- scales::rescale(legend_values, to = c(0.5, 2))

# Add legend for the point size scale beside elevation scale
legend("topright",                     # Position
       inset = c(0.05, 0.03),            # Moves the legend outside the map extent
       legend = legend_values, 
       pch = 19, col = "black", 
       pt.cex = legend_sizes, 
       title = "Species Richness",
       title.cex = 0.7,                # Adjust legend title size
       bty = "y",                      # Place box around legend
       xpd = TRUE)                     # Allows placing legend outside the plot

Rename the elevation raster layer.

names(dtm_crop) <- "elevation"

Extract elevation values for each bird observation record and append it to mauricie_birds. First, make sure the CRS of both spatial data match and then use st_extract().

# Match the CRS
mauricie_birds_sf <- st_transform(mauricie_birds_sf, crs = st_crs(dtm_crop))

# Extract the elevation values for each site
elevation_values_sf <- stars::st_extract(dtm_crop, mauricie_birds_sf)

Add in the data columns we want

elevation_values_sf <- bind_cols(elevation_values_sf, 
                                 st_drop_geometry(mauricie_birds_sf)) %>%
  dplyr::select(SiteCode, elevation)

Assign a point id identifier to each location based on its unique geometry and return the result as a dataframe.

elevation_values_df <- elevation_values_sf %>% 
  group_by(SiteCode, geometry) %>%
  mutate(point_id = cur_group_id()) %>% 
  select(point_id, elevation) %>% 
  distinct() %>%
  ungroup() %>%
  sf::st_drop_geometry() %>%  # drop the geometry column
  as.data.frame()  # convert to a regular dataframe


# Note that usually the SiteCode is the unique identifier for each observation site. However, this was not the case for the Quebec Breeding Bird Atlas data. We therefore use the geometry to create a unique identifier for each point.

OPTIONAL: To save the elevation data as a .csv to your disk, use the write.csv() function, specify the name of your .csv file, and use the row.names = FALSE argument to exclude row numbers from the output. The extracted elevation values will be required to complete Chapter 7: Summary Tools. We have uploaded this file to the Google Drive data folder for your convenience.

Create a subdirectory called env_covariates to neatly store the output environmental data.

dir.create("data/env_covariates", recursive = TRUE)

To execute this code chunk, remove the #

# write.csv(elevation_values_df, "data/env_covariates/elevation_values_df.csv", row.names = FALSE)

Now, let’s combine the elevation dataframe with the NatureCounts summary data you previously created in order to assess how elevation influences your metric of species diversity.

mauricie_site_summary <- st_transform(mauricie_site_summary, crs = st_crs(dtm_crop))

# Extract the elevation values for each site
elevation_sitesum_sf <- stars::st_extract(dtm_crop, mauricie_site_summary) %>%
  bind_cols(st_drop_geometry(mauricie_site_summary))

# Plot the relationship between elevation and species diversity 
plot(elevation_sitesum_sf$elevation, elevation_sitesum_sf$n_species, 
     xlab = "Elevation (m)", ylab = "Number of Species", 
     main = "Species Diversity vs. Elevation")

There does not appear to be a clear relationship between species diversity and elevation in La Mauricie National Park.

---

Congratulations! You completed Chapter 4 - Digital Elevation Data. In this chapter, you successfully 1) processed raster DTM’s 2) performed cropping and masking procedures and 3) extracted elevation data over vector data points. In Chapter 5, you will explore how to extract unique pixel values from landcover data and link it to NatureCounts data.

Below are some additional resources you may find useful to further your understanding of LiDAR data.

4.5 LiDAR Data Resources

LiDAR systems help produce extremely accurate measurements of surface elevation to digital terrain and elevation models across landscapes. These systems use instruments (sensors) that are commonly flown on aircraft and incorporate laser ranging for distance, GPS for geographic position and height, and aircraft altitude for the orientation of the instrument. The accuracy of LiDAR measurements depends on the strength of point returns, instrument precision & accuracy, and terrain characteristics, among other factors.

Conceptual figure depicts an airborne LiDAR system (‘LiDAR-i lend’, Marek9134 (2012) CC-BY-SA-3.0).

As of June 2022, Natural Resources Canada has produced their largest release of LiDAR-derived data made publicly available on Open Maps. This data is complementary to previous CanElevation Series products: the HRDEM and HRDEM Mosaic. The HRDEM download procedure can be found here. Some notable open LiDAR data resources on the provincial-level include LidarBC - British Columbia and Forêt ouverte - Québec.

The USGS provides LiDAR and its derived products through the LidarExplorer online application. The National Map also represents the primary repository for geospatial data in the U.S.

4.6 Manual Data Download

If you wish to access LiDAR data for La Mauricie National Park yourself to gain experience with the process, follow the steps below. The data are available on the Forêt ouverte platform.

Step 1: Navigate to the Lidar data site. Under Données et ressources > Carte interactive > Forêt ouverte select (Explorer > Accéder) next to the Téléchargement option to explore the Lidar tiles on an interactive map .

Step 2: Using the search bar, search for Parc national de la Mauricie. Under Layers, toggle on Parc nationaux du Canada to visualize the park boundary.

Alternatively, you can import your own .shp or .kml file using the Import and Export tab.

Specify the coordinate reference system (WGS84), click Import, and select your .shp or .kml file to display on the map.

Step 3: Download the Modèle num. terrain (MNT) (Résolution spatiale 1 m) for each of the six Lidar tiles that intersect with the park (14NE, 15NO, 14SE, 15SO, 11NE, 10NO).

First, navigate to the Map tab and toggle on the visibility for the boundary and Lidar MNT layers by clicking on the eye symbol, if necessary. Hovering over this icon will display either the Show Layer or Hide layer option.

Your interactive map should now look something like this:

One by one, click and download the terrain data for each of the 6 Lidar tiles that intersect with the park boundary. Click your target tile, expand the Téléchargement window and click Accéder next to the MNT (1 m spatial resolution) data to download. Repeat for each tile.

This will open a webpage directory displaying a list of spatial data files available for the selected tile. Choose the MNT (.tif) file that includes the tile code in its name. For instance, for the first tile covering northern tip of the national park), select “MNT_31I14NE.tif.” Repeat this process for each respective tile.

Save these files in your working directory for easy access during this tutorial.

getwd()