Tour of LidarBC https://lidar.gov.bc.ca/

Look at different layers, inspect the index. Where is the data? What are the file formats?

Download the .laz file over the university

From the Point Cloud 1:2,500 layer.
Save to K:/Lab3/

Download the .tif file over the university

From the DEM 1:20,000 layer.

Save to K:/Lab3/

QGIS Plugins

Using Plugins > Manage and Install Plugins, add the folowing plugins in QGIS:

  • Quick Map Service in order to add a basemap
  • Profile Tool to plot cross sections of elevation values

Add layers to Map

Just drag and drop the layers to the map.

Color by Z

Using the Layer Styling Panel color the point cloud with the different point cloud attributes.

3D View

Using the View > 3D Maps View > New 3D Map View, look at your data in 3D.

Cross Section Profile

  1. Add the DEM to the map.
  2. Change the styling to a hillshade.
  3. Open the profile tool.
  4. Add the DEM tif to the cross section profile tool.
  5. Draw cross sections and explore the data.

Bulk Download with R

Tour of `R`

Open R

Install and Load Packages

# INSTALL LIBRARIES
  install.packages("arcpullr")
  install.packages("mapview")
  install.packages("bcmaps")
  install.packages("dplyr")
  install.packages("sf")

# LOAD LIBRARIES
  library(arcpullr)
  library(mapview)
  library(bcmaps)
  library(dplyr)
  library(sf)  

Ouput Directory

# DEST DIR 
  destdir <- "C:/Users/bevington/Downloads/Lab3/"
  dir.create(destdir) 

BC Maps

# AOI 
  aoi <- bcmaps::bc_cities() %>% filter(NAME == "Prince George") %>% 
    st_buffer(5000) %>%
    st_transform(4326)
  mapview(aoi)

Define ArcGIS REST end point

# ARCGIS REST FEATURE SERVER URL
LiDAR_DEM_Index_20000 <- "https://services6.arcgis.com/ubm4tcTYICKBpist/arcgis/rest/services/LiDAR_BC_S3_Public/FeatureServer/6"

Query ArcGIS REST Tiles

# GET LIST OF TILES
  tiles <- arcpullr::get_layer_by_poly(LiDAR_DEM_Index_20000, aoi, sp_rel = "intersects")
  print(tiles)
  mapview(aoi) + mapview(tiles)

Download Tiles

# DOWNLOAD TILES 
  for(i in 1:nrow(tiles)){
    print(i)
    filename <- tiles[i,]$filename
    url <- tiles[i,]$s3Url
    download.file(url, paste0(destdir, filename), mode="wb")
  }

Stich Together Tiles into a Virtual Raster

# MOSAIC TILES 
  sf::gdal_utils(util = "buildvrt", 
                 source = list.files(destdir, full.names = T, pattern = "*.tif$"), 
                 destination = paste0(destdir,"dem.vrt"))

lidR Package

# READ LAS FILE
  # install.packages("lidR")
  library(lidR)
  las <- readLAS(paste0(destdir, "bc_093g086_4_4_4_xyes_8_utm10_2019.laz"))
  print(las)

Canopy Height Model

  # install.packages("terra")
  library(terra)

  # MAKE 2x2m surface model
  dsm = pixel_metrics(las = las, func = ~max(Z), res = 2)
  plot(dsm)

  # MAKE 2x2m terrain model
  dtm = pixel_metrics(las = las, func = ~min(Z), res = 2)
  plot(dtm)
  
  # Make Canopy Height Model
  chm = dsm-dtm
  plot(chm)

  # Export
  terra::writeRaster(dsm, paste0(destdir,"dsm.tif"))
  terra::writeRaster(dtm, paste0(destdir,"dtm.tif"))
  terra::writeRaster(chm, paste0(destdir,"chm.tif"))

Questions:

  • How many points are in the bc_093g086_4_4_4_xyes_8_utm10_2019.laz file? What is the point density?
  • Name 3 strengths and 3 challenges of using lidar data.
  • Look at the canopy height model. Name two artifacts that are visible in this dataset.
  • In QGIS, using the profile tool, make a cross section of the DTM and DSM in a single plot (as shown below), and in a seperate plot make a cross section of the CHM. Make sure you use a line that crosses the forest and over the UNBC buildings and down the cutbanks. Use the following colors. DTM (RED), DSM (GREEN), CHM (BLUE). Roughly how tall are the tallest trees in your CHM cross section?
Categories: GEOG 457Labs