Lab 04 – LiDAR

Introduction

RStudio

https://www.rstudio.com/

lidR

https://github.com/r-lidar/lidR

Open a las file

Open RStudio then create a new Project

Create a New Directory, and New Project. Name it Lab4 and save it in your GEOG457 Folder.

Create a new script to put your work into (The page with the green + is the new file button, choose R Script as type).

• Open the library lidR
• Read a file using readLAS() and save it into a variable named las
• Print information about the file we opened to the console.

library(lidR)
#las <- readLAS("L:/GEOG457/Lab4/PrinceGeorge/pg_2010/pg_2010_093G08644.las")
# Open decimated version instead to prevent computer crash
las <- readLAS("L:/GEOG457/Lab4/1m_las_tile/pg_2010/pg_2010_093G08644.las")
print(las)

To run your code click on the first line and press [Ctrl] + [Enter] or press the run button

We will wait here and get everyone caught up to having information about the file printed to the screen

# Open a directory with LASCatalog
ctg <- readLAScatalog("K:/GEOG457/AlezaLake/las_tile_reproject/")

# For every las file in the catalog perform an operation
for (i in 1:(length(ctg))) {
  tmp <- ctg[i,]
  
  # Computationally determin new file name
  # (Technically you could just overwrite the original file if you live dangerous)
  part_path = strsplit(tmp$filename, "las_tile_reproject/aleza")[[1]]
  new_path = paste("L:/GEOG457/Lab4/AcientForest/las_tile/af", part_path[2], sep = "")
  
  # Open the Las File
  las <- readLAS(tmp$filename)
  
  # Set the projection (This does not modify points only helps 
  # the computer to know what the locations mean!)
  # The projection of 32610 (UTM 10N, WGS 84) is prior knowlege 
  # must be obtained from producer of .las file
  epsg(las) <- 32610 
  
  # Save the file, also lets add a .lax or index file 
  # to make future catalog operations faster
  writeLAS(las, new_path, index = TRUE)
}

Many times LiDAR data is processed as a Local Coordinate System (LCS), this is to say the computer does not care where the points are in the world only in the cloud, you will notice that extent has no units attached only numbers. One easy wat to see this is in an open-source program that we will not be using today called Cloud Compare, when you open a LiDAR file it typically will shift all points on import, allow you to do your work, then apply the reverse shift when you export the file, this is intended to help keep precision by reducing the number of significant digits.

To get a 3D representation of the point cloud you can simply plot it. Have some patience it may take some time to think and may open in the background (watch your taskbar for a new icon to appear)

plot(las)

gc() #Free up RAM used by plot()

Digital Elevation Models

Digital Elevation Models or DEM’s are something you should have interacted with by now and are often generated from LiDAR. However when possible avoid saying DEM there are multiple types, and at the scales, we are working with it matters.

The two most common are Digital Terrain Models (DSM): These are DEM’s where the elevation is the ground without any Vegetation or Structures. On the other hand Digital Surface Models are created from the first return (highest point) and show trees and buildings.

There are more forms of DEM’s you may make for domain speific purposes, for example a Canopy Height Model (CHM) which is the difference between the DSM and DTM, (Technically a CHM should also filture out buildings, however that is extra effort rarely required).

In your assignment please do not refer to speific products as DEM’s they should all be DTM, DSM, or CHM.

Generate a DTM:

dtm <- grid_terrain(las, res = 1, algorithm = tin())
plot(dtm)

Generate a DSM:

thr <- c(0,2,5,10,15)
edg <- c(0, 1.5)
dsm <- grid_canopy(las, res = 1, pitfree(thr, edg))

Generate a CHM

For this we have two options, first is the simple way. Though could potentially be more affected by the effects of raster resolution.

chm <- dsm - dtm

The second way is to make a normalized point cloud then produce the CHM from that. We will revisit this after switching over the Ancient forest where we will be using higher resolution data and more of it.

LasCatalog

LAS files are very large (so much so we have concerns about crashing the lab computers, however, we also want to do tree detection today so want the full resolution.

As a rough overview of what this does, we will open a folder instead of a file, and then the computer can work on the entire area one piece at a time so we hopefully have enough ram. Additionally, many options can be run in parallel on a multicore computer in a catalogue. Below is an example of what processing looks like.

Open the Acient Forest data

Be very careful here attempting to work with the point data of the entire region will crash your computer (que reminder to save your script!)

# Enample Multicore processing using future package
library(future)
plan(multisession)

#open data
acient_forest_ctg <- readLAScatalog("L:/GEOG457/Lab4/AcientForest/las_tile/")
plot(acient_forest_ctg)

Would you like a map with that?

plot(acient_forest_ctg, map=TRUE)

You can use the catalog like a las file in most lidR functions. And often time with better results:

https://jean-romain.github.io/lidRbook/engine.html#engine-rationale

It is also convient, not sure exactly what tile you need? Just extract an area from the catalog, lets take a look at UNBC

pg_2010_ctg <- readLAScatalog("L:/GEOG457/Lab4/1m_las_tile/pg_2010/")
plot(clip_circle(pg_2010_ctg, x = 512179, y = 5971394, radius = 500), size = 2)
gc()

Tree Detection

acient_forest_ctg <- readLAScatalog("L:/GEOG457/Lab4/AcientForest/las_tile/")
plot(acient_forest_ctg)

# Subset
acient_forest_subset <- clip_rectangle(acient_forest_ctg, xleft = 616950, xright = 617500, ytop = 5958786, ybottom = 5958100)
# Generate Products for later
dtm <- grid_terrain(acient_forest_subset, algorithm = tin())

thr <- c(0,2,5,10,15)
edg <- c(0, 1.5)
dsm <- grid_canopy(acient_forest_subset, 1, pitfree(thr, edg))


# Normalize Data
af_normalized <- normalize_height(acient_forest_subset, dtm)
plot(af_normalized)

# Better Canopy height model
chm <- grid_canopy(af_normalized, 1, pitfree(thr, edg))
plot(chm, col = height.colors(50))



#Find Trees
ttops <- find_trees(af_normalized, lmf(ws = 10))
plot(chm, col = height.colors(50))
plot(ttops, add = TRUE, col = 'Red')


#Segement Trees
af_trees <- segment_trees(af_normalized, li2012(R = 3, speed_up = 5))
col <- random.colors(200)
tree_segments <- plot(af_trees, color = "treeID")
add_treetops3d(tree_segments, ttops)

# Deliniate Crowns
crowns <- delineate_crowns(af_trees)
plot(crowns)

Assignment (5 points)

For this assignment create an R script, each question should start with a comment of the question number. There is no need to submit any plots, just the R file, please put written answers in comments using the # symbol.

# Boilerplate ####
library(lidR)
library(future)
plan(multisession)

# Lab4 ####
## Question 1 ####
Code......

## Question 2 ####
Code....

## Question 3 ####
Code......


## Question 4 ####
Code......

Tip: The #### at the end of a comment denotes it as part of an outline, the number of leading #’s sets the heading level. You can click on the outline to automatically scroll to the code in question.

Questions

• (0.5 Points) There is a bridge over the Fraser River at the coordinates 53°54’36″N 122°43’24″W, what is the file name of the full resolution Prince George 2018 las file that contains this location?
  Expected Computation Time: less than 1 minute
• (1 Point) Create a DEM of Prince George (yes all of it). Your output should be at a 25m resolution and should be created from the 2014 data from the 10m LiDAR, use the algorithm knnidw(), tin() will run but takes about 3x as long though with a few less holes in the river. Note: you can safely ignore the warnings about points outside of the convex hull.
  
  Export your image as a GeoTIFF (See line 2 of code at bottom of page). No need to submit image so long as code is working.
  Expected Computation Time: 5-6 minutes

• (1.5 Points) Make a subset of the 1 metre resolution 2018 Prince George data in Lheidli T’enneh Memorial Park, roughly the same area as the blue rectangle in the picture below. If you are not familiar with the area it is 53°54’23″N 122°44’03″W. TIPs: At some point in this question you will want this equitation: lmf(ws = 15, hmin = 20), when running the tree_segmentation() use dalponte2016 it runs substantially faster than li2012; if you do not know how to use dalponte2016 as this is the first time seeing it, fear not we have the book.
  Finally I would advise against plotting the point clouds, and only plotting the raster products.
  
  • Make a Canopy Height Model (from a normalized point cloud)
    You may get an excessive amount of warnings, again these again be ignored for our purposes.
  • Make a plot with markings over the tree tops
  • Make a 2nd plot with the crowns outlined

Expected Computation Time: 10-15 minutes

• (1 Point) Create a Plot showing the change between 2018 and 2010 of the Prince George Cutbanks 53°55’45″N 122°45’34″W
  Expected Computation Time: 2-3 Minutes.

• (0.5 Points) Why do we prefer LASCatalog over individual large LAS files? And why do we prefer it over simply applying operations to each file individually.
• (0.5 Point) Well organized code; part of being a skilled professional is making your work understandable to others (or yourself in a few months), clean code is also easier to debug. Before submitting your assignment make sure that the R file is well organized, and variables have logical easy to understand names.

Not part of the lab but possibly useful commands

# Save a geoTIFF
writeRaster(dtm, "dtm.tif", options=c('TFW=YES'))


# Generate Hill shade
slope <- terrain(dtm, opt="slope")
aspect <- terrain(dtm, opt="aspect")
hill_shade <- hillShade(slope, aspect) # hill_shade variable has underscore, don't name variables same as function
plot(hill_shade, col = grey.colors(256))