Code to create an animation explaining how individual tree detection works (using a local maximum filter)

ITD_animation.R

library(lidR)
library(tidyverse)
library(glue)
library(ggforce) #geom_circle
library(gifski)

#set your paths here
dir_out <- "C:/Users/ptompals/OneDrive - NRCan RNCan/__WORK/graphics"

#temp dir will be a subfolder of the main dir
dir_temp_frames <- file.path(dir_out, "frames_temp")


#using the built-in lidR dataset in this example
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, filter = "-inside 481280 3812940 481330 3812990")

#generate a CHM
chm <- rasterize_canopy(las, res = 0.5, pitfree(c(0,2,5,10,15), c(0, 1.5)))

#crop the CHM further. This can be done with the readLAS call as well. Advantage 
# was that I could select a nice looking subset for the visualization from the larger CHM
chm <- terra::crop(chm, ext(481280, 481310 ,3812950,  3812970))


#define the variable moving window size function. Using the function provided in lidR locate_trees() example
f <- function(x) { x * 0.07 + 3}

#run tree detection
ttops <- locate_trees(chm, lmf(ws = f))

#check results
# plot(chm, col = col)
# plot(ttops, add=T, col="black")


# link treetops to chm pixels
ttops_extract <- terra::extract(chm, ttops, xy=T, ID=T)

#convert CHM to a data frame (to the plot using ggplot, easier than dealing with rasters)
R <- as.data.frame(chm, xy=T)

#calculate window size (mimicing what lidR::locate_trees does)
R <- R %>% mutate(ws = f(Z)) #ws is in meters (same units as the data)

#calculate window extent
R <- R %>% mutate(ws_x_min = x - ws/2,
                  ws_x_max = x + ws/2,
                  ws_y_min = y - ws/2,
                  ws_y_max = y + ws/2)

#calculate radius for circular window
R <- R %>% mutate(radius = ws / 2)

# create a second df to show tree detection on a subset of the CHM (entire data would take too long)
R1 <- R
R1 <- R1 %>% mutate(id = row_number())
R1 <- R1[450:1450,] #modify this range to change what is shown in the animation


# copy x and y coords of treetops to two different fields (so that the treetops are separated from CHM pixel coords)
ttops_extract <- ttops_extract %>% mutate(x_top = x, y_top = y)

#link R1 with treetops
R1 <- R1 %>% left_join(ttops_extract, by = join_by(x, y))

#cleanup
R1 <- R1 %>% select(-Z.y) %>% rename(Z = Z.x)

buffer <- 2.5 #to extend plot size and show the moving window when it is partially outside the CHM


#create frames one by one




# circular

dir_temp_frames_circular <- file.path(dir_temp_frames, "circular")

if(dir.exists(dir_temp_frames_circular)) unlink(dir_temp_frames_circular)
dir.create(dir_temp_frames_circular, recursive = T)


for (i in 1:nrow(R1)) {
  
  pix <- R1[i,]
  
  current_tops <- R1[1:i,]
  
  
  p <- 
    
    #draw the background CHM first
    ggplot(data=R, aes(x, y, fill=Z)) +
    geom_raster() +
    
    #coord fixed so that pixels are square
    coord_fixed()+
    
    #color palette for the CHM
    scale_fill_viridis_c()+
    
    #add a point showing current pixel...
    geom_point(data=pix, aes(x, y), color="red")+
    
    #... and the moving window aroung it
    # two options here: moving window can be a square or a circle
    
    #if square - use the commented-out code below
    # geom_rect(data=pix,
    #           aes(xmin = ws_x_min , xmax = ws_x_max , ymin = ws_y_min, ymax = ws_y_max),
    #           fill=NA, color="red")+
    
    #if circle:
    geom_circle(data=pix, aes(x0=x, y0=y, r=radius),fill=NA, color="red")+
    
    #show treetops if they were already detected
    geom_point(data=current_tops, aes(x=x_top, y=y_top), color="black", inherit.aes = F, size=4, shape=4)+
    
    #extend plot boundaries to show moving window outside the plot as well
    xlim(min(R$x)-buffer, max(R$x)+buffer )+
    ylim(min(R$y)-buffer, max(R$y)+buffer )+
    
    #remove plot elements
    theme_bw()+
    theme(panel.grid = element_blank(),
          panel.border = element_blank(),
          axis.title = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          legend.position = "none")
  
  fout <- file.path(dir_temp_frames_circular, glue("frame{sprintf('%04d', i)}.png"))
  
  ggsave(plot = p, 
         filename = fout,
         width = 4, height = 3)
}

#make a gif
png_files <- list.files(dir_temp_frames_circular, full.names = T)
gif_file <- file.path(dir_out, "ITD_animation.gif")
gifski(png_files, gif_file, delay = 0.05)