======================================================
Code snippet from Wetzin'kwa Community Forest Project:
  Tree detection and segmentation
======================================================
library(lidR)
library(terra)
library(sf)
library(tidyverse)
#=====================================
#STEP 1: IMPORT LIDAR DATA
#=====================================
#import lidar data
las_data <- readLAScatalog("Data/LAS_data")       
#Optional: check lidar data and view CRS
las_check(las_data)
st_crs(las_data)
            
#==========================================================
#STEP 2: CONVERT CRS and CLIP LAZ FILES TO AREA OF INTEREST
#==========================================================
#convert KML files to same CRS
stand_kml_converted <- st_transform(stand_union_kml, st_crs(las_data))
#clip .las data to the roi
las_clipped <- clip_roi(las_data, stand_kml_converted)
#=======================================================
#STEP 3: NORMALIZE HEIGHTS AND FILTER OUTLIERS 
#=======================================================
#normalize the data using DEM
las_norm <- normalize_height(las_clipped, wetzinkwa_dem)
#filter undesired data points
range(las_norm$Z, na.rm = TRUE)
las_norm<- filter_poi(las_norm, Z >= 0 & Z <= 20)
#======================================
#STEP 4: CREATE CANOPY HEIGHT MODEL
#======================================
#create simple CHM with a 0.5m resolution using p2r CHM
chm <- rasterize_canopy(las_norm, res = 0.5, algorithm = p2r(0.15))
  plot(chm)
#apply a mean smoothing filter
chm_smoothed <- focal(chm, w = matrix(1,3,3), fun = mean, na.rm = TRUE)
#plot a 3D chm
plot_dtm3d(chm, color="green")
#=======================================
#STEP 5: TREE SEGMENTATION
#=======================================
#define a variable window size  
f <- function(h) {
    pmax(1.0, pmin(6.0, 1.5 + 0.05 * h))}
#locate trees
ttops <- locate_trees(chm, lmf(f))
  plot(chm)
    plot(ttops, add=TRUE, col="red")
  
crowns_dal <- segment_trees(las_norm, dalponte2016(chm_smoothed, ttops, th_tree = 6))
  plot(crowns_dal, color = "treeID")

=====================================================
Code Snippet from Fire Analysis from Landsat Data
=====================================================
#generate a fire polygon for the Prouton fire
prouton_fire_polygon <- fire_polygons %>% 
  filter(FIRE_NO == "C30870") %>% 
  summarise()
#list.dirs for Landsat folders
folders <- list.dirs("data/Landsat 8 OLI_TIRS C2 L2", full.names = TRUE, recursive = FALSE)
#create an output directory
out_dir <- "data/output"
dir.create(out_dir, showWarnings = FALSE)
#create one large for loop to iterate through the folders,
#separate files, crop and mask to QA PIXEL, calculate ndvi and nbr, 
#then save outputs into three separate folders
for (f in folders) {
  prod_id <- basename(f)
  
  #Extract Date Using lubridate
  date_string <- stringr::str_sub(prod_id, 18, 25)
  acq_date <- lubridate::as_date(date_string, format = "%Y%m%d")
  # Create a clean filename base
  clean_name <- format(acq_date, "%Y%m%d")
  files <- list.files(f, full.names = TRUE, pattern = "TIF$", recursive = TRUE)
  # extract sr + qa files
  # Select only SR_B1 ... SR_B7
  sr_files <- files[grepl("SR_B[1-7]", files)]
 #build sr stack
  sr_stack <- terra::rast(sr_files)
  qa_files <- terra::rast(files[grepl("QA_PIXEL", files)])
  
  # Set CRS and crop rasters to fire polygon
  prouton_fire_polygon <- st_transform(prouton_fire_polygon, crs(sr_stack))
  
  # crop by polygon and then create the clear mask
  qa_crop <- terra::crop(qa_files, prouton_fire_polygon, mask = TRUE)
  clear_mask <- qa_crop == 21824
  
  #crop sr files by polygon, then mask using the clear mask
  sr_crop <- terra::crop(sr_stack, prouton_fire_polygon, mask = TRUE)
  sr_mask <- mask(sr_crop, clear_mask, maskvalues = 0)
  
  # calculate ndvi and nbr
  ndvi <- (sr_mask[[5]] - sr_mask[[4]]) / (sr_mask[[5]] + sr_mask[[4]])
  nbr  <- (sr_mask[[5]] - sr_mask[[7]]) / (sr_mask[[5]] + sr_mask[[7]])
  
  # create directories 
  sr_dir   <- file.path(out_dir, "SR")
  ndvi_dir <- file.path(out_dir, "NDVI")
  nbr_dir  <- file.path(out_dir, "NBR")
  dir.create(sr_dir, showWarnings = FALSE)
  dir.create(ndvi_dir, showWarnings = FALSE)
  dir.create(nbr_dir, showWarnings = FALSE)
  
  # write with date
  writeRaster(sr_mask, file.path(sr_dir, paste0(clean_name,"_SR.tif")),
              overwrite = TRUE)
  writeRaster(ndvi, file.path(ndvi_dir, paste0(clean_name,"_NDVI.tif")),
              overwrite = TRUE)
  writeRaster(nbr, file.path(nbr_dir, paste0(clean_name,"_NBR.tif")),
              overwrite = TRUE)