Skip to content

Instantly share code, notes, and snippets.

@ptompalski
Created February 2, 2024 03:25
Show Gist options
  • Save ptompalski/94904eca2e1628fb52010c2890431715 to your computer and use it in GitHub Desktop.
Save ptompalski/94904eca2e1628fb52010c2890431715 to your computer and use it in GitHub Desktop.
Code to create an animation explaining how individual tree detection works (using a local maximum filter)
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment