The 2002 Hayman Fire (Graham 2003; Moriarty & Cheng, 2012) was the largest in Colorado history at the time and it severely impacted the Cheesman Reservoir and watershed. The reservoir, along with ~8,000 ac surrounding, is managed by Denver Water and is adjacent to the Pike San Isabel National Forest. The high severity fire resulted in near total mortality of the ponderosa pine and Douglas-fir forest surrounding the reservoir which necessitated active planting of seedlings through partnerships between Denver Water and the US Forest Service, Colorado State Forest Service, National Resource Conservation Service and the Colorado Forest Restoration Institute through the From Forests to Faucets program (Denver Water).

The objective of this project is to demonstrate the use of UAS remote sensing to evaluate the survival and growth of the reforestation efforts surrounding the Cheesman Reservoir. An exploratory site visit was made on April 03, 2026 by representatives from Denver Water (Jessica Jackman), CSFS (Ethan Bucholz, Andy Whelan, and Spencer Weston), CFRI (Caden Chamberlain), and CSU (George Woolsey). A single, preliminary UAS flight was completed with the DJI Mavic xx over one of the sites that had some of the first plantings. The flight took ~10 minutes with flight settings programmed to capture data for SfM processing to generate aerial point cloud data and a spectral orthomosaic.

Preliminary Data

let’s check out the preliminary data captured from a single UAS flight that was completed on April 03, 2026 with the DJI Mavic xx over one of the sites that had some of the first plantings. The flight took ~10 minutes with flight settings programmed to capture data for SfM processing to generate aerial point cloud data and a spectral orthomosaic.

load the standard libraries we use to do work

# bread-and-butter
library(tidyverse) # the tidyverse
library(viridis) # viridis colors
library(harrypotter) # hp colors
library(RColorBrewer) # brewer colors
library(scales) # work with number and plot scales
library(latex2exp)

# visualization
library(mapview) # interactive html maps
library(kableExtra) # tables
library(patchwork) # combine plots
library(ggnewscale) # new scale
library(ggrepel) # repel labels

# spatial analysis
library(terra) # raster
library(sf) # simple features
library(lidR) # lidR
library(cloud2trees) # cloud2trees

Though not necessary for cloud2trees data processing, let’s quickly check out the location and structure of the data we have

# directory with the downloaded .las|.laz files
point_cld_folder <- "../data/CheesmanLAS/"
# is there data?
list.files(point_cld_folder, pattern = ".*\\.(laz|las)$") %>% length()
## [1] 1
# what files are in here?
list.files(point_cld_folder, pattern = ".*\\.(laz|las)$")[1]
## [1] "cloudR.las"

again, this is not necessary for cloud2trees data processing but we can use lidR to read the point cloud folder as a catalog which doesn’t read in the actual points but just the point cloud header data which includes information on things like the spatial location of the data, the point density, and other point attributes

# read folder as LAScatalog
ctg_temp <- lidR::readLAScatalog(point_cld_folder)
# what information do we get about the point cloud?
ctg_temp
## class       : LAScatalog (v1.2 format 2)
## extent      : 473435.3, 474225.5, 4339062, 4339903 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 13N 
## area        : 0.66 km²
## points      : 38.68 million points
## type        : airborne
## density     : 58.2 points/m²
## num. files  : 1

let’s look at the point cloud extent on a map to orient ourselves in space

ctg_temp %>% 
  cloud2trees:::check_las_ctg_empty() %>% 
  purrr::pluck("data") %>% 
  mapview::mapview(popup = F, layer.name = "point cloud tile")

i told you that we didn’t need to do any of that for cloud2trees data processing and to prove it, we’ll remove the ctg_temp object from our session

remove(list = ls()[grep("_temp",ls())])
gc()
##           used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 3998997 213.6    7749105 413.9  5168294 276.1
## Vcells 5855892  44.7   12255594  93.6  8388599  64.0

we also got an orthomosaic of at least RGB spectral bands. we might also have gotten 4-band imagery if the Mavic xx used for the flight included a sensor that also captures the NIR band?

#### pc extent vs rgb rast
rgb_rast_fnm <- "../data/CheesmanOrtho/result.tif"
rgb_rast <- terra::rast(rgb_rast_fnm)

what is this raster data?

rgb_rast
## class       : SpatRaster 
## size        : 18519, 17021, 4  (nrow, ncol, nlyr)
## resolution  : 0.0306, 0.0306  (x, y)
## extent      : 473599.2, 474120, 4339087, 4339654  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 13N (EPSG:32613) 
## source      : result.tif 
## names       : result_1, result_2, result_3, result_4

it looks like there might actually be a fourth band…

terra::hist(rgb_rast)

…or not. that fourth band is probably the alpha channel. let’s check out the RGB

rgb_rast <- rgb_rast %>% terra::subset(c(1:3))
# get them from the cloud2trees hidden package
# plot
terra::plotRGB(rgb_rast)

# terra::plot(
#   rgb_rast
#   , nc = 3
#   # , nr = 3
#   , mar = c(0.5,0.5,0.5,0.5)
#   , axes = FALSE
#   , legend = F
#   # , col = grDevices::gray.colors(n=111)
# )

looks pretty good and there are certainly shadows on the ground from trees (hopefully). let’s zoom in on the 5000 square meter central area

rgb_rast %>%
  terra::crop(
    terra::ext(rgb_rast[[1]]) %>% 
      sf::st_bbox() %>% 
      sf::st_as_sfc() %>% 
      sf::st_centroid() %>% 
      sf::st_buffer(sqrt(5000)/2) %>% 
      sf::st_bbox() %>% 
      sf::st_as_sfc() %>% 
      sf::st_set_crs(terra::crs(rgb_rast)) %>% 
      terra::vect()
    , mask = T
  ) %>% 
  terra::plotRGB()

green trees!

Point Cloud to Raster: cloud2raster()

Although the cloud2trees::cloud2trees() function combines methods in the cloud2trees package for an all-in-one approach, we’ll instead use the cloud2trees::cloud2raster() function to generate a CHM from the point cloud that we can work with to perform both individual tree detection and slash pile identification (discussed later).

We’ll set the options in the function to generate a CHM which represents a DSM with the ground removed and no other filtering. This high resolution (i.e. fine-grain) CHM will serve as the foundation for tree detection and slash pile detection as we can manipulate it to optimize the processing for both methods.

# outdir
c2t_output_dir <- "../data"
c2t_process_dir <- file.path(c2t_output_dir, "point_cloud_processing_delivery")
c2t_tracking_fnm <- file.path(c2t_process_dir, "processed_tracking_data.csv")
##############################################################
# cloud2trees::cloud2raster
##############################################################
if(
  !file.exists( file.path(c2t_process_dir, "chm_0.1m.tif") )
  || !file.exists( file.path(c2t_process_dir, "dtm_0.25m.tif") )
){
  # time it
  st_temp <- Sys.time()
  # run it
  # cloud2trees
  cloud2raster_ans <- cloud2trees::cloud2raster(
    output_dir = c2t_output_dir
    , input_las_dir = point_cld_folder
    , accuracy_level = 2
    , keep_intrmdt = T
    , dtm_res_m = 0.25
    , chm_res_m = 0.1
    , min_height = 0 # effectively generates a DSM based on non-ground points
  )
  
  # timer
  mins_temp <- difftime(Sys.time(),st_temp,units = "mins") %>% as.numeric()
  # save tracking
  dplyr::tibble(
    timer_cloud2raster_mins = mins_temp
  ) %>% 
    write.csv(
      file = c2t_tracking_fnm
      , row.names = F, append = F
    )
}else{
  dtm_temp <- terra::rast( file.path(c2t_process_dir, "dtm_0.25m.tif") )
  chm_temp <- terra::rast( file.path(c2t_process_dir, "chm_0.1m.tif") )
  
  cloud2raster_ans <- list(
    "dtm_rast" = dtm_temp
    , "chm_rast" = chm_temp
  )
}

let’s see what we got from cloud2trees::cloud2raster()

cloud2raster_ans %>% names()
## [1] "dtm_rast" "chm_rast"

there’s a DTM

# plot to check out the fine-resolution DTM raster
cloud2raster_ans$dtm_rast %>% 
  terra::plot(col = harrypotter::hp(n=100, option = "mischief"), main = "DTM (m)", axes = F)

there’s a CHM

# plot to check out the fine-resolution CHM raster
cloud2raster_ans$chm_rast %>% 
  terra::plot(col = viridis::plasma(n=100), main = "CHM (m)", axes = F)

let’s see some details about the CHM

# what chm?
cloud2raster_ans$chm_rast
## class       : SpatRaster 
## size        : 8411, 7902, 1  (nrow, ncol, nlyr)
## resolution  : 0.1, 0.1  (x, y)
## extent      : 473435.3, 474225.5, 4339062, 4339903  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 13N (EPSG:32613) 
## source      : chm_0.1m.tif 
## name        :          chm 
## min value   : 4.609016e-05 
## max value   : 2.564435e+01
# what data?
cloud2raster_ans$chm_rast %>% terra::summary()
##       chm       
##  Min.   : 0.02  
##  1st Qu.: 0.67  
##  Median : 1.12  
##  Mean   : 1.43  
##  3rd Qu.: 1.74  
##  Max.   :19.94  
##  NA's   :98391

for tree detection, we’ll aggregate the CHM to a lower resolution (i.e. coarser) to smooth out some of the fine detail that can increase the noise in the tree detection processing. This aggregation process also will speed up the processing and reduce the chances of computational memory issues.

##################################################################
# aggregate to make raster more coarse for tree detection
##################################################################
# first, we'll borrow from the `cloud2trees` codebase to get a function to change the resolution of a raster exactly
###___________________________________________###
# adjust the resolution of a raster to be in exactly the target resolution
###___________________________________________###
if(
  !file.exists( file.path(c2t_process_dir, "chm_0.20m.tif") )
){
  agg_chm_rast <- cloud2trees:::adjust_raster_resolution(
    cloud2raster_ans$chm_rast
    , target_resolution = 0.20
    , fun = max
    , resample_method = "max"
    , ofile = file.path(c2t_process_dir, "chm_0.20m.tif")
  )
}else{
  agg_chm_rast <- terra::rast( file.path(c2t_process_dir, "chm_0.20m.tif") )
}
# what chm?
agg_chm_rast
## class       : SpatRaster 
## size        : 4205, 3951, 1  (nrow, ncol, nlyr)
## resolution  : 0.2, 0.2  (x, y)
## extent      : 473435.3, 474225.5, 4339062, 4339903  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 13N (EPSG:32613) 
## source      : chm_0.20m.tif 
## name        :          chm 
## min value   : 1.460902e-04 
## max value   : 2.564435e+01
# what data?
agg_chm_rast %>% terra::summary()
##       chm       
##  Min.   : 0.01  
##  1st Qu.: 0.70  
##  Median : 1.18  
##  Mean   : 1.51  
##  3rd Qu.: 1.85  
##  Max.   :19.94  
##  NA's   :97832

let’s go back to the central 5000 square meter area and look at the CHM over the RGB

first, we’ll plot the 0.1 m resolution CHM overlaid

# aoi_temp
aoi_temp <- terra::ext(rgb_rast[[1]]) %>% 
  sf::st_bbox() %>% 
  sf::st_as_sfc() %>% 
  sf::st_centroid() %>% 
  sf::st_buffer(sqrt(5000)/2) %>% 
  sf::st_bbox() %>% 
  sf::st_as_sfc() %>% 
  sf::st_set_crs(terra::crs(rgb_rast)) %>% 
  terra::vect()
rgb_rast %>%
  terra::crop(aoi_temp, mask = T) %>% 
  terra::plotRGB(
    mar = c(0.1,0.1,1.5,5)
    , main = "0.1 m resolution CHM"
  )
cloud2raster_ans$chm_rast %>%
  # agg_chm_rast %>% 
  terra::crop(aoi_temp, mask = T) %>% 
  terra::plot(
    add = T
    , axes = F
    , col = viridis::plasma(n=100) # harrypotter::hp(n=100, option = "gryffindor")
    , alpha = 0.4
    , plg = list(title = "CHM (m)", title.cex = 0.9, shrink = 0.6)
  )

now, we’ll plot the 0.2 m resolution CHM overlaid

rgb_rast %>%
  terra::crop(aoi_temp, mask = T) %>% 
  terra::plotRGB(
    mar = c(0.1,0.1,1.5,5)
    , main = "0.2 m resolution CHM"
  )
agg_chm_rast %>% 
  terra::crop(aoi_temp, mask = T) %>% 
  terra::plot(
    add = T
    , axes = F
    , col = viridis::plasma(n=100) # harrypotter::hp(n=100, option = "gryffindor")
    , alpha = 0.4
    , plg = list(title = "CHM (m)", title.cex = 0.9, shrink = 0.6)
  )

…interesting

a few initial observations and thoughts:

  • The timing of the data collection (early April) allows the green conifers to stand out against the background terrain and any dormant or senescent grasses and forbs
    • Tree segmentation could start with identifying green objects from the RGB and then confirm with structural data (or vice-versa)
  • There are standing dead trees which are represented in the upper limits of the CHM (>8 m)
    • These will need to be removed in an analysis of the tree plantings via spectral filtering and/or height filtering
    • Could also filter based on some height/area ratio
  • The taller tree plantings (about 1.5-4.5 meters tall) are clearly visible in the CHM
    • What window function to use for ITD from CHM to separate trees with overlapping crowns?
  • Some of the shorter tree plantings (< 1.5 m) are not represented in the CHM (no CHM pixels overlapping green areas)
    • Since they are not in the CHM, it is not likely they are in the point cloud since this CHM is effectively a height-normalized DSM based on non-ground points (cloud2raster(..., min_height = 0))
    • Explore identifying these potential trees first using RGB, then check the raw point cloud in these areas for segmentation and height
    • If these short trees are not represented in the point cloud, try different SfM processing settings?
    • If these short trees are represented in the point cloud, would need to integrate a tree segmentation method that uses point cloud data directly
    • If these short trees are not represented in the point cloud even after using different SfM processing, use the RGB for segmentation and shadows for height?