lidR4smartieS

This document has three exercises that will teach you the basics of working in R with lidar data using the lidR package.

You can view this online at https://rpubs.com/jesseast/lidR4smartieS

Exercise 1

In 2021 PJ woodland vegetation data was collected in the form of 30 m x 30 m field plots in Contact, Nevada. For each lidar tile and each plot, generate a digital terrain model (DTM), digital surface model (DSM) and a canopy height model (CHM). For the purpose of this exercise, all the data has been downloaded and processed and is in the folder “lidar4smarties”, but you should follow along with the data acquisition steps when appropriate.

The plot you will be using for this exercise is entitled “UOFU_NVCO_P07”. The plot spatial information is in the form of a polygon feature in the shapefile nvco7shp. The shapefile also contains the plot center coordinates and the coordinates of the four plot corners.

plot_id	plot_y	plot_x	plot_ne_y	plot_ne_x	plot_se_y	plot_se_x	plot_sw_y	plot_sw_x	plot_nw_y	plot_nw_x
UOFU_NVCO_P07	41.6388547	-114.4650434	41.6389866	-114.4648553	41.6387179	-114.4648729	41.6387285	-114.4652269	41.6389905	-114.4652248

Now that we have a shapefile that contains the polygon of the field plot, it is necessary to acquire high quality, recent lidar that intersects the plot polygon.

We do not know what lidar tiles intersect with our plot polygon, therefore, we will go to the national map downloader and upload the zipped shapefile nvco7shp.zip using the “upload shapefile” button under “Advanced Search”. The polygon will appear on the map and you can zoom into the plot on the map. For “Area of Interest” select “Extent”. Select “Elevation Source Data (3DEP) - Lidar, IfSAR” and click the blue “Search Products” button. After searching products, you should see the resulting lidar tiles selected based on the query. In the “Collapse View” box click the “Show All Footprints” and squares should appear on the map indicating the extent of a lidar tile. The result of this search should result should be one lidar tile that is spatially coincident to our plot polygon, entitled NV_NorthWestElko_2020_D20 11TQG110125. The lidar tile can be downloaded by selecting “Download Link (LAZ)” and moving the LAZ file to your R directory. This lidar file is already downloaded for you in the “/las_files” folder.

Now, we will use the lidR package to process our lidar data. There are two options for reading las/laz files in the lidR package: readLAS and readLAScatalog. readLAS is for reading an individual las file and readLAScatalog is used for reading a collection of multiple las files You may need to use readLAScatalog if a plot polygon is not entirely covered by a single lidar tile or if you want to process multiple plot polygons across different lidar tiles. For now, only one plot that encompasses a single las file will be processed, therefore we can use the readLAS function.

First we need to read in our libraries.

library(lidR)
library(sf)
library(terra)
library(mapview)
library(tmaptools)
library(tmap)

Here we read our las file NV_NorthWestElko_2020_D20 11TQG110125 and plot polygon shapefile nvco7shp into R. Next we filter out duplicate points from the point cloud and remove any noise using the filter_duplicates() and filter_noise() lidR functions, respectively.

Here is the step by step of the preprocessing process

las = readLAS("S:/ursa/eastburn/URSA_work/Tutorials/lidR4smartieS/las_files/USGS_LPC_NV_NorthWestElko_2020_D20_11TQG110125.laz") #filter = "-keep_first" #Use readLAS to read your las file.

## Warning: There are 1706954 points flagged 'withheld'.

## Warning: There are 1706954 points flagged 'synthetic'.

polygon = sf::st_read("S:/ursa/eastburn/URSA_work/Tutorials/lidR4smartieS/nvco7shp/plot_level_dat_NVCO7_1983_11.shp") #read the shapefile of your plot polygon

## Reading layer `plot_level_dat_NVCO7_1983_11' from data source 
##   `S:\ursa\eastburn\URSA_work\Tutorials\lidR4smartieS\nvco7shp\plot_level_dat_NVCO7_1983_11.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 30 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 711115.6 ymin: 4612769 xmax: 711146.4 ymax: 4612799
## Projected CRS: NAD83 / UTM zone 11N

las = filter_duplicates(las) #filter out duplicated points if necessary. can use the function las_check() to see if there are any duplicated points in your point cloud

las = filter_poi(las, Classification != LASNOISE) #remove noise from las file

Then we convert the plot polygon shapefile into the correct coordinate system +init=EPSG:26911, if it was not correct already. We then set the coordinate system of the las file to the same coordinate system as the plot polygon shapefile. Lastly, we clip the lidar point cloud to the extent of the plot polygon.

lidR::crs(las) #check crs of the las 

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +vunits=m +no_defs 
## WKT2 2019 representation:
## COMPOUNDCRS["NAD83(2011) / UTM zone 11N + NAVD88 height - Geoid18 (Meters)",
##     PROJCRS["NAD83(2011) / UTM zone 11N",
##         BASEGEOGCRS["NAD83(2011)",
##             DATUM["NAD83 (National Spatial Reference System 2011)",
##                 ELLIPSOID["GRS 1980",6378137,298.257222101,
##                     LENGTHUNIT["metre",1]]],
##             PRIMEM["Greenwich",0,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             ID["EPSG",6318]],
##         CONVERSION["UTM zone 11N",
##             METHOD["Transverse Mercator",
##                 ID["EPSG",9807]],
##             PARAMETER["Latitude of natural origin",0,
##                 ANGLEUNIT["degree",0.0174532925199433],
##                 ID["EPSG",8801]],
##             PARAMETER["Longitude of natural origin",-117,
##                 ANGLEUNIT["degree",0.0174532925199433],
##                 ID["EPSG",8802]],
##             PARAMETER["Scale factor at natural origin",0.9996,
##                 SCALEUNIT["unity",1],
##                 ID["EPSG",8805]],
##             PARAMETER["False easting",500000,
##                 LENGTHUNIT["meter",1],
##                 ID["EPSG",8806]],
##             PARAMETER["False northing",0,
##                 LENGTHUNIT["meter",1],
##                 ID["EPSG",8807]]],
##         CS[Cartesian,2],
##             AXIS["easting",east,
##                 ORDER[1],
##                 LENGTHUNIT["meter",1]],
##             AXIS["northing",north,
##                 ORDER[2],
##                 LENGTHUNIT["meter",1]],
##         ID["EPSG",6340]],
##     VERTCRS["NAVD88 height - Geoid18 (Meters)",
##         VDATUM["North American Vertical Datum 1988"],
##         CS[vertical,1],
##             AXIS["up",up,
##                 LENGTHUNIT["meter",1]],
##         GEOIDMODEL["GEOID18"],
##         ID["EPSG",5703]]]

lidR::crs(polygon) #check crs of the plot polygon

## [1] "PROJCRS[\"NAD83 / UTM zone 11N\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"UTM zone 11N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-117,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    ID[\"EPSG\",26911]]"

#They should be in the same crs, if not, use the project function to put the plot and/or las file in the correct crs. 

#For example:
#polygon = project(polygon, "+init=EPSG:26911") # or use st_transform from the sf package

#Once both are in the same crs, there might be small differences in syntax that can make the las and polygon not compatible, therefore we set the crs equal to each other. Note that this does not project the polygon or las file, it essentially just overwrites the CRS information

lidR::st_crs(las) = lidR::st_crs(polygon)  #ensure the point cloud crs is the same crs as the plot polygon (small differences in syntax can make this whole thing not work even though they are technically in the same CRS)

Next we can clip the las file to our polygon shapefile using the clip_roi() function.

lidar.plot = lidR::clip_roi(las, polygon) #clip the lidar point cloud to the extent of the plot polygon
lidR::plot(lidar.plot) #view the lidar tile and plot

Now that we have done some preprocessing, we can focus on some deliverables! The goal of this exercise, which we will complete now, is to create a DTM, DSM, and CHM, for both the plot polygon area and the lidar tile area.

DTM

The DTM can be calculated using the rasterize_terrain function. To calculate the DTM for the entire lidar tile area, use the “las” variable as your data. To calculate the DTM for the plot polygon area, use the clipped lidar point cloud “lidar.plot” as your data. The rasterize_canopy() function also takes the arguments “res”, which is the resolution of the raster, and the tin() argument, which is a spatial interpolation algorithm that is based on Delaunay triangulation. See the lidR documentation for additional spatial interpolation algorithms that can be used to derive a dtm.

DTM for tile

dtm = rasterize_terrain(las, res = 1, tin()) #las is your las file, res is the resolution of the output, and tin() is a spatial interpolation algorithm

## Delaunay rasterization[===================-------------------------------] 39% (18 threads)
Delaunay rasterization[====================------------------------------] 40% (18 threads)
Delaunay rasterization[====================------------------------------] 41% (18 threads)
Delaunay rasterization[=====================-----------------------------] 42% (18 threads)
Delaunay rasterization[=====================-----------------------------] 43% (18 threads)
Delaunay rasterization[======================----------------------------] 44% (18 threads)
Delaunay rasterization[======================----------------------------] 45% (18 threads)
Delaunay rasterization[=======================---------------------------] 46% (18 threads)
Delaunay rasterization[=======================---------------------------] 47% (18 threads)
Delaunay rasterization[========================--------------------------] 48% (18 threads)
Delaunay rasterization[========================--------------------------] 49% (18 threads)
Delaunay rasterization[=========================-------------------------] 50% (18 threads)
Delaunay rasterization[=========================-------------------------] 51% (18 threads)
Delaunay rasterization[==========================------------------------] 52% (18 threads)
Delaunay rasterization[==========================------------------------] 53% (18 threads)
Delaunay rasterization[===========================-----------------------] 54% (18 threads)
Delaunay rasterization[===========================-----------------------] 55% (18 threads)
Delaunay rasterization[============================----------------------] 56% (18 threads)
Delaunay rasterization[============================----------------------] 57% (18 threads)
Delaunay rasterization[=============================---------------------] 58% (18 threads)
Delaunay rasterization[=============================---------------------] 59% (18 threads)
Delaunay rasterization[==============================--------------------] 60% (18 threads)
Delaunay rasterization[==============================--------------------] 61% (18 threads)
Delaunay rasterization[===============================-------------------] 62% (18 threads)
Delaunay rasterization[===============================-------------------] 63% (18 threads)
Delaunay rasterization[================================------------------] 64% (18 threads)
Delaunay rasterization[================================------------------] 65% (18 threads)
Delaunay rasterization[=================================-----------------] 66% (18 threads)
Delaunay rasterization[=================================-----------------] 67% (18 threads)
Delaunay rasterization[==================================----------------] 68% (18 threads)
Delaunay rasterization[==================================----------------] 69% (18 threads)
Delaunay rasterization[===================================---------------] 70% (18 threads)
Delaunay rasterization[===================================---------------] 71% (18 threads)
Delaunay rasterization[====================================--------------] 72% (18 threads)
Delaunay rasterization[====================================--------------] 73% (18 threads)
Delaunay rasterization[=====================================-------------] 74% (18 threads)
Delaunay rasterization[=====================================-------------] 75% (18 threads)
Delaunay rasterization[======================================------------] 76% (18 threads)
Delaunay rasterization[======================================------------] 77% (18 threads)
Delaunay rasterization[=======================================-----------] 78% (18 threads)
Delaunay rasterization[=======================================-----------] 79% (18 threads)
Delaunay rasterization[========================================----------] 80% (18 threads)
Delaunay rasterization[========================================----------] 81% (18 threads)
Delaunay rasterization[=========================================---------] 82% (18 threads)
Delaunay rasterization[=========================================---------] 83% (18 threads)
Delaunay rasterization[==========================================--------] 84% (18 threads)
Delaunay rasterization[==========================================--------] 85% (18 threads)
Delaunay rasterization[===========================================-------] 86% (18 threads)
Delaunay rasterization[===========================================-------] 87% (18 threads)
Delaunay rasterization[============================================------] 88% (18 threads)
Delaunay rasterization[============================================------] 89% (18 threads)
Delaunay rasterization[=============================================-----] 90% (18 threads)
Delaunay rasterization[=============================================-----] 91% (18 threads)
Delaunay rasterization[==============================================----] 92% (18 threads)
Delaunay rasterization[==============================================----] 93% (18 threads)
Delaunay rasterization[===============================================---] 94% (18 threads)
Delaunay rasterization[===============================================---] 95% (18 threads)
Delaunay rasterization[================================================--] 96% (18 threads)
Delaunay rasterization[================================================--] 97% (18 threads)
Delaunay rasterization[=================================================-] 98% (18 threads)
Delaunay rasterization[=================================================-] 99% (18 threads)
Delaunay rasterization[==================================================] 100% (18 threads)

lidR::plot(dtm) #plot the dtm

DTM for plot

dtm.plot = rasterize_terrain(lidar.plot, res = 1, tin()) #las is your las file, res is the resolution of the output, and tin() is a spatial interpolation algorithm

lidR::plot(dtm.plot) #plot the dtm

DSM

The DSM can be calculated using the rasterize_canopy function. The dsmtin() argument is a spatial interpolation algorithm that is based on Delaunay triangulation. See the lidR documentation for additional spatial interpolation algorithms.

DSM for tile

dsm = rasterize_canopy(las, res = 1, algorithm = dsmtin()) #las is your las file, res is the resolution of the output, and dsmtin() is a spatial interpolation algorithm

plot(dsm) #plot the dsm

DSM for plot

dsm.plot = rasterize_canopy(lidar.plot, res = 1, dsmtin()) #las is your las file, res is the resolution of the output, and dsmtin() is a spatial interpolation algorithm

plot(dsm.plot) #plot the dsm

Height Normalization and the CHM

In order to calculate the CHM we first need to normalize our point cloud. Normalizing a point cloud makes the point cloud ground flat by removing terrain. Then we can use the rasterize_canopy() function on the normalized point cloud, similar to how we calculated the DSM.

Note: When the original (non-normalized) point cloud with absolute elevations is used, the derived layer represents the elevation of the top of the canopy above sea level, and is referred to as DSM. So, to calculate the CHM, we need to get thee DSM of a normalized point cloud! The CHM can also be referred to as a normalized DSM or nDSM.

Height Normalization for tile

nlas = normalize_height(las, res = 1, tin(), dtm = dtm) #normalize the height of the tile using the dtm

## Delaunay rasterization[========------------------------------------------] 17% (18 threads)
Delaunay rasterization[=========-----------------------------------------] 18% (18 threads)
Delaunay rasterization[=========-----------------------------------------] 19% (18 threads)
Delaunay rasterization[==========----------------------------------------] 20% (18 threads)
Delaunay rasterization[==========----------------------------------------] 21% (18 threads)
Delaunay rasterization[===========---------------------------------------] 22% (18 threads)
Delaunay rasterization[===========---------------------------------------] 23% (18 threads)
Delaunay rasterization[============--------------------------------------] 24% (18 threads)
Delaunay rasterization[============--------------------------------------] 25% (18 threads)
Delaunay rasterization[=============-------------------------------------] 26% (18 threads)
Delaunay rasterization[=============-------------------------------------] 27% (18 threads)
Delaunay rasterization[==============------------------------------------] 28% (18 threads)
Delaunay rasterization[==============------------------------------------] 29% (18 threads)
Delaunay rasterization[===============-----------------------------------] 30% (18 threads)
Delaunay rasterization[===============-----------------------------------] 31% (18 threads)
Delaunay rasterization[================----------------------------------] 32% (18 threads)
Delaunay rasterization[================----------------------------------] 33% (18 threads)
Delaunay rasterization[=================---------------------------------] 34% (18 threads)
Delaunay rasterization[=================---------------------------------] 35% (18 threads)
Delaunay rasterization[==================--------------------------------] 36% (18 threads)
Delaunay rasterization[==================--------------------------------] 37% (18 threads)
Delaunay rasterization[===================-------------------------------] 38% (18 threads)
Delaunay rasterization[===================-------------------------------] 39% (18 threads)
Delaunay rasterization[====================------------------------------] 40% (18 threads)
Delaunay rasterization[====================------------------------------] 41% (18 threads)
Delaunay rasterization[=====================-----------------------------] 42% (18 threads)
Delaunay rasterization[=====================-----------------------------] 43% (18 threads)
Delaunay rasterization[======================----------------------------] 44% (18 threads)
Delaunay rasterization[======================----------------------------] 45% (18 threads)
Delaunay rasterization[=======================---------------------------] 46% (18 threads)
Delaunay rasterization[=======================---------------------------] 47% (18 threads)
Delaunay rasterization[========================--------------------------] 48% (18 threads)
Delaunay rasterization[========================--------------------------] 49% (18 threads)
Delaunay rasterization[=========================-------------------------] 50% (18 threads)
Delaunay rasterization[=========================-------------------------] 51% (18 threads)
Delaunay rasterization[==========================------------------------] 52% (18 threads)
Delaunay rasterization[==========================------------------------] 53% (18 threads)
Delaunay rasterization[===========================-----------------------] 54% (18 threads)
Delaunay rasterization[===========================-----------------------] 55% (18 threads)
Delaunay rasterization[============================----------------------] 56% (18 threads)
Delaunay rasterization[============================----------------------] 57% (18 threads)
Delaunay rasterization[=============================---------------------] 58% (18 threads)
Delaunay rasterization[=============================---------------------] 59% (18 threads)
Delaunay rasterization[==============================--------------------] 60% (18 threads)
Delaunay rasterization[==============================--------------------] 61% (18 threads)
Delaunay rasterization[===============================-------------------] 62% (18 threads)
Delaunay rasterization[===============================-------------------] 63% (18 threads)
Delaunay rasterization[================================------------------] 64% (18 threads)
Delaunay rasterization[================================------------------] 65% (18 threads)
Delaunay rasterization[=================================-----------------] 66% (18 threads)
Delaunay rasterization[=================================-----------------] 67% (18 threads)
Delaunay rasterization[==================================----------------] 68% (18 threads)
Delaunay rasterization[==================================----------------] 69% (18 threads)
Delaunay rasterization[===================================---------------] 70% (18 threads)
Delaunay rasterization[===================================---------------] 71% (18 threads)
Delaunay rasterization[====================================--------------] 72% (18 threads)
Delaunay rasterization[====================================--------------] 73% (18 threads)
Delaunay rasterization[=====================================-------------] 74% (18 threads)
Delaunay rasterization[=====================================-------------] 75% (18 threads)
Delaunay rasterization[======================================------------] 76% (18 threads)
Delaunay rasterization[======================================------------] 77% (18 threads)
Delaunay rasterization[=======================================-----------] 78% (18 threads)
Delaunay rasterization[=======================================-----------] 79% (18 threads)
Delaunay rasterization[========================================----------] 80% (18 threads)
Delaunay rasterization[========================================----------] 81% (18 threads)
Delaunay rasterization[=========================================---------] 82% (18 threads)
Delaunay rasterization[=========================================---------] 83% (18 threads)
Delaunay rasterization[==========================================--------] 84% (18 threads)
Delaunay rasterization[==========================================--------] 85% (18 threads)
Delaunay rasterization[===========================================-------] 86% (18 threads)
Delaunay rasterization[===========================================-------] 87% (18 threads)
Delaunay rasterization[============================================------] 88% (18 threads)
Delaunay rasterization[============================================------] 89% (18 threads)
Delaunay rasterization[=============================================-----] 90% (18 threads)
Delaunay rasterization[=============================================-----] 91% (18 threads)
Delaunay rasterization[==============================================----] 92% (18 threads)
Delaunay rasterization[==============================================----] 93% (18 threads)
Delaunay rasterization[===============================================---] 94% (18 threads)
Delaunay rasterization[===============================================---] 95% (18 threads)
Delaunay rasterization[================================================--] 96% (18 threads)
Delaunay rasterization[================================================--] 97% (18 threads)
Delaunay rasterization[=================================================-] 98% (18 threads)
Delaunay rasterization[=================================================-] 99% (18 threads)
Delaunay rasterization[==================================================] 100% (18 threads)

Calculate and plot the CHM of the tile

chm = rasterize_canopy(nlas, res = 1, algorithm = dsmtin()) #use the normalized height tile (nlas) to calculate the CHM

chm[chm<0]=0 #some values could be below zero due to noise, so we remove them here
tm_shape(chm)+
tm_raster(style= "quantile", n=5, palette=get_brewer_pal("Greens", n=5, plot=FALSE))+
tm_layout(legend.outside = TRUE)

## stars object downsampled to 1000 by 1000 cells. See tm_shape manual (argument raster.downsample)

Height Normalization for the plot

nlas.plot = normalize_height(lidar.plot, res = 1, tin(), dtm = dtm.plot) #Normalize the plot using the clipped plot and the dtm of the plot

Calculate and plot the CHM of the plot

chm.plot = rasterize_canopy(nlas.plot, res = 1, algorithm = dsmtin()) #Make CHM using normalized las, and dsmtin() algorithm

plot(chm.plot)

Exercise 2

Situations could arise where you want to process multiple lidar tiles. For this exercise, we want to retrieve the DTM, DSM, and CHM of multiple lidar tiles using the readLAScatalog function.

Now we need to read in the downloaded las files. When working with more than one las file, it is necessary to use the readLAScatalog function.

las.cat = readLAScatalog("S:/ursa/eastburn/URSA_work/Tutorials/lidR4smartieS/las_files/") #read the folder where multiple las files are located using the readLAScatalog function

Use the following function to filter duplicate points for a las catalog. This function can also be used for a las file, however the function filter_duplicates() is also available to filter single las files (as seen in Exercise 1).

filter_dups = function(las)
{
  if (is(las, "LAS"))
  {
      las = filter_duplicates(las)
      return(las)
  }
  
  if (is(las, "LAScatalog"))
  {
     options = list(
       need_output_file = TRUE,    # Throw an error if no output template is provided
       need_buffer = FALSE)         # Throw an error if buffer is 0
     res <- catalog_map(las, filter_dups, .options = options)
     return(res)
  }
}

opt_chunk_buffer(las.cat) <- 50 #chunk_buffer and chunk_size options help with processing by breaking up the lidar tile into smaller...chunks!
opt_chunk_size(las.cat) <- 500 
opt_output_files(las.cat) <- file.path(tempdir(), "temp") #writes raster chunks to a directory so you do not need to store them in memory 
las.cat = filter_dups(las.cat) #name variable las.cat so we are working with the filtered las catalog in the following steps 

## Warning: There are 3099 points flagged 'withheld'.

## Warning: There are 3099 points flagged 'synthetic'.

## Chunk 1 of 36 (2.8%): state <U+26A0>

## Warning: There are 398012 points flagged 'withheld'.

## Warning: There are 398012 points flagged 'synthetic'.

## Chunk 2 of 36 (5.6%): state <U+26A0>

## Warning: There are 9116 points flagged 'withheld'.

## Warning: There are 9116 points flagged 'synthetic'.

## Chunk 3 of 36 (8.3%): state <U+26A0>

## Warning: There are 453695 points flagged 'withheld'.

## Warning: There are 453695 points flagged 'synthetic'.

## Chunk 4 of 36 (11.1%): state <U+26A0>

## Warning: There are 381932 points flagged 'withheld'.

## Warning: There are 381932 points flagged 'synthetic'.

## Chunk 5 of 36 (13.9%): state <U+26A0>

## Warning: There are 303146 points flagged 'withheld'.

## Warning: There are 303146 points flagged 'synthetic'.

## Chunk 6 of 36 (16.7%): state <U+26A0>

## Warning: There are 4871 points flagged 'withheld'.

## Warning: There are 4871 points flagged 'synthetic'.

## Chunk 7 of 36 (19.4%): state <U+26A0>

## Warning: There are 429006 points flagged 'withheld'.

## Warning: There are 429006 points flagged 'synthetic'.

## Chunk 8 of 36 (22.2%): state <U+26A0>

## Warning: There are 93005 points flagged 'withheld'.

## Warning: There are 93005 points flagged 'synthetic'.

## Chunk 9 of 36 (25%): state <U+26A0>

## Warning: There are 503504 points flagged 'withheld'.

## Warning: There are 503504 points flagged 'synthetic'.

## Chunk 10 of 36 (27.8%): state <U+26A0>

## Warning: There are 434933 points flagged 'withheld'.

## Warning: There are 434933 points flagged 'synthetic'.

## Chunk 11 of 36 (30.6%): state <U+26A0>

## Warning: There are 236694 points flagged 'withheld'.

## Warning: There are 236694 points flagged 'synthetic'.

## Chunk 12 of 36 (33.3%): state <U+26A0>

## Warning: There are 1803 points flagged 'withheld'.

## Warning: There are 1803 points flagged 'synthetic'.

## Chunk 13 of 36 (36.1%): state <U+26A0>

## Warning: There are 442906 points flagged 'withheld'.

## Warning: There are 442906 points flagged 'synthetic'.

## Chunk 14 of 36 (38.9%): state <U+26A0>

## Warning: There are 251965 points flagged 'withheld'.

## Warning: There are 251965 points flagged 'synthetic'.

## Chunk 15 of 36 (41.7%): state <U+26A0>

## Warning: There are 499007 points flagged 'withheld'.

## Warning: There are 499007 points flagged 'synthetic'.

## Chunk 16 of 36 (44.4%): state <U+26A0>

## Warning: There are 435023 points flagged 'withheld'.

## Warning: There are 435023 points flagged 'synthetic'.

## Chunk 17 of 36 (47.2%): state <U+26A0>

## Warning: There are 144547 points flagged 'withheld'.

## Warning: There are 144547 points flagged 'synthetic'.

## Chunk 18 of 36 (50%): state <U+26A0>

## Warning: There are 9930 points flagged 'withheld'.

## Warning: There are 9930 points flagged 'synthetic'.

## Chunk 19 of 36 (52.8%): state <U+26A0>

## Warning: There are 440618 points flagged 'withheld'.

## Warning: There are 440618 points flagged 'synthetic'.

## Chunk 20 of 36 (55.6%): state <U+26A0>

## Warning: There are 424764 points flagged 'withheld'.

## Warning: There are 424764 points flagged 'synthetic'.

## Chunk 21 of 36 (58.3%): state <U+26A0>

## Warning: There are 514485 points flagged 'withheld'.

## Warning: There are 514485 points flagged 'synthetic'.

## Chunk 22 of 36 (61.1%): state <U+26A0>

## Warning: There are 434560 points flagged 'withheld'.

## Warning: There are 434560 points flagged 'synthetic'.

## Chunk 23 of 36 (63.9%): state <U+26A0>

## Warning: There are 69834 points flagged 'withheld'.

## Warning: There are 69834 points flagged 'synthetic'.

## Chunk 24 of 36 (66.7%): state <U+26A0>

## Warning: There are 126375 points flagged 'withheld'.

## Warning: There are 126375 points flagged 'synthetic'.

## Chunk 25 of 36 (69.4%): state <U+26A0>

## Warning: There are 465089 points flagged 'withheld'.

## Warning: There are 465089 points flagged 'synthetic'.

## Chunk 26 of 36 (72.2%): state <U+26A0>

## Warning: There are 506882 points flagged 'withheld'.

## Warning: There are 506882 points flagged 'synthetic'.

## Chunk 27 of 36 (75%): state <U+26A0>

## Warning: There are 487111 points flagged 'withheld'.

## Warning: There are 487111 points flagged 'synthetic'.

## Chunk 28 of 36 (77.8%): state <U+26A0>

## Warning: There are 443068 points flagged 'withheld'.

## Warning: There are 443068 points flagged 'synthetic'.

## Chunk 29 of 36 (80.6%): state <U+26A0>

## Warning: There are 29494 points flagged 'withheld'.

## Warning: There are 29494 points flagged 'synthetic'.

## Chunk 30 of 36 (83.3%): state <U+26A0>

## Warning: There are 207421 points flagged 'withheld'.

## Warning: There are 207421 points flagged 'synthetic'.

## Chunk 31 of 36 (86.1%): state <U+26A0>

## Warning: There are 396228 points flagged 'withheld'.

## Warning: There are 396228 points flagged 'synthetic'.

## Chunk 32 of 36 (88.9%): state <U+26A0>

## Warning: There are 447558 points flagged 'withheld'.

## Warning: There are 447558 points flagged 'synthetic'.

## Chunk 33 of 36 (91.7%): state <U+26A0>

## Warning: There are 300975 points flagged 'withheld'.

## Warning: There are 300975 points flagged 'synthetic'.

## Chunk 34 of 36 (94.4%): state <U+26A0>

## Warning: There are 398740 points flagged 'withheld'.

## Warning: There are 398740 points flagged 'synthetic'.

## Chunk 35 of 36 (97.2%): state <U+26A0>

## Warning: There are 35673 points flagged 'withheld'.

## Warning: There are 35673 points flagged 'synthetic'.

## Chunk 36 of 36 (100%): state <U+26A0>

#try renaming the last las.cat something else

Use the following function to denoise the lidar point clouds in a las catalog. This function can also be used for a las file, however the function filter_noise() is also available to filter single las files (as seen in Exercise 1).

denoise = function(las)
{
  if (is(las, "LAS"))
  {
      las = readLAS(chunk)
      las = filter_poi(las, Classification != LASNOISE) #denoise las.cat
      return(las)
  }
  
  if (is(las, "LAScatalog"))
  {
     options = list(
       need_output_file = TRUE,    # Throw an error if no output template is provided
       need_buffer = FALSE)         # Throw an error if buffer is 0
     res <- catalog_map(las, denoise, .options = options)
     return(res)
  }
}

opt_chunk_buffer(las.cat) <- 50 
opt_chunk_size(las.cat) <- 500 
opt_output_files(las.cat) = file.path(tempdir(), "temp8") #writes raster chunks to a directory so you do not need to store them in memory 
las.cat = denoise(las.cat) #name variable las.cat so we are working with the filtered and denoised las catalog in the following steps 

## Warning: There are 1190736 points flagged 'withheld'.

## Warning: There are 1190736 points flagged 'synthetic'.

## Warning: There are 1190736 points flagged 'withheld'.

## Warning: There are 1190736 points flagged 'synthetic'.

## Chunk 1 of 1 (100%): state <U+26A0>

DTM

Make DTM for LASCatalog object

opt_chunk_buffer(las.cat) <- 50 
opt_chunk_size(las.cat) <- 500
opt_output_files(las.cat) = file.path(tempdir(), "temp.00") #set directory for output. Can use a temporary directory as well using tempdir() 
dtm.cat.2 = rasterize_terrain(las.cat, res = 5, tin())

## Chunk 1 of 1 (100%): state <U+2713>

#visualize DTM
tm_shape(dtm.cat.2)+
tm_raster(style= "quantile", palette=get_brewer_pal("Greens", plot=FALSE))+
tm_layout(legend.outside = TRUE)

DSM

Make DSM for LASCatalog object

opt_chunk_buffer(las.cat) <- 50 
opt_chunk_size(las.cat) <- 500 
opt_output_files(las.cat) = file.path(tempdir(), "temp22") #set directory for output. Can use a temporary directory as well using tempdir() 
dsm.cat.2 = rasterize_canopy(las.cat, res = 5, algorithm = p2r())

## Chunk 1 of 1 (100%): state <U+26A0>

#visualize dsm
tm_shape(dsm.cat.2)+
tm_raster(style= "quantile", n=7, palette=get_brewer_pal("Greens", n=7, plot=FALSE))+
tm_layout(legend.outside = TRUE)

Height Normalization

Normalize LASCatalog object using normalize_height() function

opt_chunk_size(las.cat) <- 500
opt_output_files(las.cat) = file.path(tempdir(), "temp3") #set directory for output. Can use a temporary directory
las.cat.norm = normalize_height(las.cat, res = 5, tin(), dtm = dtm.cat.2) 

## Warning: There are 1190052 points flagged 'withheld'.

## Warning: There are 1190052 points flagged 'synthetic'.

## Chunk 1 of 1 (100%): state <U+26A0>

CHM

Make the CHM of the lidar tile and plot. Make sure to use normalized las catalog for the CHM.

opt_chunk_buffer(las.cat.norm) <- 50 
opt_chunk_size(las.cat.norm) <- 500 
opt_output_files(las.cat.norm) <- file.path(tempdir(), "temp44")
chm.cat1 = rasterize_canopy(las.cat.norm, res = 5, algorithm = dsmtin()) #use the normalized height lidar tile to calculate the CHM. Can decrease the res argument for finer resolution.

## Warning: There are 1190052 points flagged 'withheld'.

## Warning: There are 1190052 points flagged 'synthetic'.

## Chunk 1 of 1 (100%): state <U+26A0>

#visualize the chm
chm.cat1[chm.cat1<0]=0
tm_shape(chm.cat1)+
tm_raster(style= "quantile", n=5, palette=get_brewer_pal("Greens", n=5, plot=FALSE))+
tm_layout(legend.outside = TRUE)

If we want to write any raster to a tiff, say the CHM, we can use the function writeRaster and a filepath.

f <- file.path(tempdir(), "chmcat1.tif") #make sure to specify format in file name (.tif)
terra::writeRaster(x = chm.cat1, filename = f, overwrite = TRUE) #Save the raster file with the function writeRaster(‘variable name’, filename = ‘save path’)

Exercise 3

We want to determine the canopy density and the canopy cover for the UOFU_NVCO_P07 plot. User derived metrics can be calculated pretty easily for both las and lascatalog objects. For Exercise 3 we will create a user defined function to calculate canopy density and canopy cover for both LAS and LAScatalog objects.

Canopy density is defined as the number of all returns at or above breast height (1.37 m) divided by total number of all returns. Meanwhile, canopy cover is defined as the number of first returns at or above breast height (1.37 m) divided by total number of first returns.

First, we will create our own function to calculate canopy metrics.

met.fun = function(Z){
  n.pts = length(Z) #number of points in point cloud
  n.pts.above.1.37m = length(Z[Z>=1.37]) #number of points above 1.37m in height
  ratio.above.1.37m = n.pts.above.1.37m / n.pts #number of points above 1.37m in height divided by number of points
  met.list = list(
    number.of.points = n.pts,
    number.of.points.above.1.37.meter = n.pts.above.1.37m,
    ratio.of.points.above.1.37.meter = ratio.above.1.37m)
  return(met.list)
}

Canopy density

#Single plot and LAS object
mets1 = cloud_metrics(nlas.plot, ~met.fun(Z)) #use normalized point cloud file to calculate canopy density for a single plot
mets1$ratio.of.points.above.1.37.meter

## [1] 0.09106884

mets2 = cloud_metrics(nlas, ~met.fun(Z)) #use normalized point cloud file to calculate canopy density for a tile
mets2$ratio.of.points.above.1.37.meter

## [1] 0.02119882

Canopy Cover Use filter to filter normalized point cloud to first returns to calculate canopy cover

#Single plot and LAS object
mets3 = cloud_metrics(nlas.plot, ~met.fun(Z), filter = ~ReturnNumber == 1
) #use normalized point cloud file to calculate canopy density for a single plot
mets3$ratio.of.points.above.1.37.meter

## [1] 0.09363666

mets4 = cloud_metrics(nlas, ~met.fun(Z), filter = ~ReturnNumber == 1
) #use normalized point cloud file to calculate canopy density for a single plot
mets4$ratio.of.points.above.1.37.meter

## [1] 0.02037484

Bonus

Voxelize points

Sometimes you want to have fun with voxels! To derive metrics from voxels, simply use the voxel_metrics function which converts a las file into voxels at your chosen resolution then calculates your function for each voxels, in this case we are using the canopy density function from exercise 3. To voxelize a lascatalog, you can adapt the function we used to denoise and filter duplicates by replacing the denoise() or filter_duplicates() with voxelize_metrics().

vox_met <- voxel_metrics(las, ~met.fun(Z), res = 4) # calculate voxel metrics. 
plot(vox_met, color="Z", size = 4, bg = "white", voxel = TRUE)

Classify ground points

If you are working with point cloud data downloaded on the national map explorer, the ground points are usually already classified. If you are working with hovermap data, ground points may not be classified. We can classify ground points using the lidR package and the function classify_ground().

las.ground.class <- classify_ground(las, algorithm = pmf(ws = 5, th = 3)) #refer to lidR book for additional algorithms

## Morphological filter: 1% (1 threads)
Morphological filter: 2% (1 threads)
Morphological filter: 3% (1 threads)
Morphological filter: 4% (1 threads)
Morphological filter: 5% (1 threads)
Morphological filter: 6% (1 threads)
Morphological filter: 7% (1 threads)
Morphological filter: 8% (1 threads)
Morphological filter: 9% (1 threads)
Morphological filter: 10% (1 threads)
Morphological filter: 11% (1 threads)
Morphological filter: 12% (1 threads)
Morphological filter: 13% (1 threads)
Morphological filter: 14% (1 threads)
Morphological filter: 15% (1 threads)
Morphological filter: 16% (1 threads)
Morphological filter: 17% (1 threads)
Morphological filter: 18% (1 threads)
Morphological filter: 19% (1 threads)
Morphological filter: 20% (1 threads)
Morphological filter: 21% (1 threads)
Morphological filter: 22% (1 threads)
Morphological filter: 23% (1 threads)
Morphological filter: 24% (1 threads)
Morphological filter: 25% (1 threads)
Morphological filter: 26% (1 threads)
Morphological filter: 27% (1 threads)
Morphological filter: 28% (1 threads)
Morphological filter: 29% (1 threads)
Morphological filter: 30% (1 threads)
Morphological filter: 31% (1 threads)
Morphological filter: 32% (1 threads)
Morphological filter: 33% (1 threads)
Morphological filter: 34% (1 threads)
Morphological filter: 35% (1 threads)
Morphological filter: 36% (1 threads)
Morphological filter: 37% (1 threads)
Morphological filter: 38% (1 threads)
Morphological filter: 39% (1 threads)
Morphological filter: 40% (1 threads)
Morphological filter: 41% (1 threads)
Morphological filter: 42% (1 threads)
Morphological filter: 43% (1 threads)
Morphological filter: 44% (1 threads)
Morphological filter: 45% (1 threads)
Morphological filter: 46% (1 threads)
Morphological filter: 47% (1 threads)
Morphological filter: 48% (1 threads)
Morphological filter: 49% (1 threads)
Morphological filter: 50% (1 threads)
Morphological filter: 51% (1 threads)
Morphological filter: 52% (1 threads)
Morphological filter: 53% (1 threads)
Morphological filter: 54% (1 threads)
Morphological filter: 55% (1 threads)
Morphological filter: 56% (1 threads)
Morphological filter: 57% (1 threads)
Morphological filter: 58% (1 threads)
Morphological filter: 59% (1 threads)
Morphological filter: 60% (1 threads)
Morphological filter: 61% (1 threads)
Morphological filter: 62% (1 threads)
Morphological filter: 63% (1 threads)
Morphological filter: 64% (1 threads)
Morphological filter: 65% (1 threads)
Morphological filter: 66% (1 threads)
Morphological filter: 67% (1 threads)
Morphological filter: 68% (1 threads)
Morphological filter: 69% (1 threads)
Morphological filter: 70% (1 threads)
Morphological filter: 71% (1 threads)
Morphological filter: 72% (1 threads)
Morphological filter: 73% (1 threads)
Morphological filter: 74% (1 threads)
Morphological filter: 75% (1 threads)
Morphological filter: 76% (1 threads)
Morphological filter: 77% (1 threads)
Morphological filter: 78% (1 threads)
Morphological filter: 79% (1 threads)
Morphological filter: 80% (1 threads)
Morphological filter: 81% (1 threads)
Morphological filter: 82% (1 threads)
Morphological filter: 83% (1 threads)
Morphological filter: 84% (1 threads)
Morphological filter: 85% (1 threads)
Morphological filter: 86% (1 threads)
Morphological filter: 87% (1 threads)
Morphological filter: 88% (1 threads)
Morphological filter: 89% (1 threads)
Morphological filter: 90% (1 threads)
Morphological filter: 91% (1 threads)
Morphological filter: 92% (1 threads)
Morphological filter: 93% (1 threads)
Morphological filter: 94% (1 threads)
Morphological filter: 95% (1 threads)
Morphological filter: 96% (1 threads)
Morphological filter: 97% (1 threads)
Morphological filter: 98% (1 threads)
Morphological filter: 99% (1 threads)
Morphological filter: 100% (1 threads)

## Original dataset already contains 18556307 ground points. These points were reclassified as 'unclassified' before performing a new ground classification.

Points are already classified as ground so the classify_ground() function will reclassify points. We can visualize the result of the ground classification by coloring points by their classification.

plot(las.ground.class, color = "Classification", size = 3)