# The r-universe version (recommended for the latest version)
install.packages('ICESat2VegR', repos = c('https://carlos-alberto-silva.r-universe.dev', 'https://cloud.r-project.org'))

# The CRAN version
#install.packages("ICESat2VegR")

# Required additional libraries
need_pkgs <- c(
  "reticulate",   # Python <-> R interface
  "leaflet",      # interactive maps
  "sf",           # spatial vectors
  "terra",        # rasters & vectors
  "data.table",   # fast tables
  "dplyr",        # tidy helpers
  "mapview",      # interactive map viewer
  "caret",        # modeling / ML
  "hdf5r"         # reading/writing local HDF5 (.h5) files
)
missing <- need_pkgs[!need_pkgs %in% rownames(installed.packages())]
if (length(missing)) {
  message("Installing missing R packages: ", paste(missing, collapse = ", "))
  install.packages(missing, repos = repos, dependencies = TRUE)
}

suppressPackageStartupMessages({
  library(ICESat2VegR)
  library(reticulate)
  library(leaflet)
  library(sf)
  library(terra)
  library(data.table)
  library(dplyr)
  library(mapview)
  library(caret)
})
cat("Loading libraries... done.\n")                                                    # Confirm loading

ICESat2VegR_configure()

# Verify configuration status

# Helper safely() runs an expression and returns a default value if an error occurs
safely <- function(expr, default = NA) tryCatch(expr, error = function(e) default)

# Check if 'reticulate' is installed
have_reticulate <- requireNamespace("reticulate", quietly = TRUE)

# Attempt to discover which Python executable is being used
py_path <- if (have_reticulate) safely(reticulate::py_discover_config()$python, NA) else NA

# Verify if Python is properly initialized and available
py_ready <- have_reticulate && !inherits(try(reticulate::py_config(), silent = TRUE), "try-error")

# Build a list summarizing the environment and module status
status <- list(
  ICESat2VegR_loaded = "package:ICESat2VegR" %in% search(),  # Is package loaded in current session?
  python_used        = py_path,                              # Path to active Python binary
  h5py               = if (py_ready) reticulate::py_module_available("h5py")        else FALSE,
  earthaccess        = if (py_ready) reticulate::py_module_available("earthaccess") else FALSE,
  ee                 = if (py_ready) reticulate::py_module_available("ee")          else FALSE
)

# Print the collected environment information to the console
print(status)

# Load the ICESat2VegR package
library(ICESat2VegR)

# Set output directory
outdir <- tempdir()

# Download example dataset
url <- "https://github.com/carlos-alberto-silva/ICESat2VegR/releases/download/example_datasets/Study_Site.zip"
zip_file <- file.path(outdir, "Study_Site.zip")
download.file(url, zip_file, mode = "wb")

# Unzip the example dataset
unzip(zip_file, exdir = outdir)

# Specifying bounding box coordinates
lower_left_lon <- -96.0
lower_left_lat <- 40.0
upper_right_lon <- -100
upper_right_lat <- 42.0


# Specifying the date range
daterange <- c("2021-10-02", "2021-10-03")

atl03_granules_local <- ATLAS_dataFinder(
short_name = "ATL03",
lower_left_lon,
lower_left_lat,
upper_right_lon,
upper_right_lat,
version = "007",
daterange = daterange,
persist = TRUE,
cloud_computing = FALSE
)

head(atl03_granules_local)

atl08_granules_local <- ATLAS_dataFinder(
  short_name = "ATL08",
  lower_left_lon,
  lower_left_lat,
  upper_right_lon,
  upper_right_lat,
  version = "007",
  daterange = daterange,
  persist = TRUE,
  cloud_computing = FALSE
)

head(atl08_granules_local)

# Download all granules
ATLAS_dataDownload(atl03_granules_local[1:3], outdir)
ATLAS_dataDownload(atl08_granules_local[c(2,4,5)],  outdir)

## ATL03
# Read the granules
atl03_files <- list.files(outdir, "ATL03.*h5", full.names = TRUE)
atl03_h5 <- lapply(atl03_files, ATL03_read)


## ATL08
# Read the granules
atl08_files <- list.files(outdir, "ATL08.*h5", full.names = TRUE)
atl08_h5 <- lapply(atl08_files, ATL08_read)

# List groups within first file of atl08_h5
atl08_h5[[1]]$ls()

# Set NASA Earthdata credentials in Python environment (required for cloud access)
# Replace "your_username" and "your_password" with your NASA Earthdata Login credentials
# Register at: https://urs.earthdata.nasa.gov
reticulate::py_run_string("
import os
os.environ['EARTHDATA_USERNAME'] = 'your_username'
os.environ['EARTHDATA_PASSWORD'] = 'your_password'
")

# Create or update .netrc file with credentials (run once per session)
earthdata_login()

# Search ATL03 granules in cloud (cloud_computing = TRUE streams data without downloading)
atl03_granules_cloud <- ATLAS_dataFinder(
  short_name      = "ATL03",
  lower_left_lon,
  lower_left_lat,
  upper_right_lon,
  upper_right_lat,
  version         = "007",   # latest version
  daterange       = daterange,
  persist         = TRUE,
  cloud_computing = TRUE
)

# Read the granule (the ATL03_read can only read one granule per read)
atl03_h5_cloud <- ATL03_read(atl03_granules_cloud[1])

# List groups within the h5 in cloud
atl03_h5_cloud$beams

## gt1l gt1r gt2l gt2r gt3l gt3r

close(atl03_h5_cloud)

# Extract the Reference Ground Track as a line object
# (rgt_extract() takes a single ATL03_read/ATL08_read product, not the list)
rgt_line <- rgt_extract(h5 = atl03_h5[[1]], line = TRUE)

# Plot the ICESat-2 orbit animation
plot_icesat2_orbit_animation(
  rgt = rgt_line,
  launch = TRUE
)

# Plot the default ICESat-2 orbit animation
plot_icesat2_orbit_animation()

atl03_photons_dt <- lapply(atl03_h5,ATL03_photons_attributes_dt)
atl03_photons_dt <- rbindlist2(atl03_photons_dt)

head(atl03_photons_dt)

# ATL03 seg attributes
atl03_seg_att_ls <- lapply(
  atl03_h5,
  ATL03_seg_metadata_dt,
  attributes = c("delta_time", "solar_elevation", "pitch", "h_ph", "ref_elev")
)
atl03_seg_dt <- rbindlist2(atl03_seg_att_ls)

# Remove segments above 20km
atl03_seg_dt <- atl03_seg_dt[h_ph < 20000]

head(atl03_seg_dt)

# ATL08 seg attributes
atl08_seg_att_ls <- lapply(
atl08_h5,
ATL08_seg_attributes_dt,
attributes = c("h_canopy", "h_te_mean", "terrain_slope", "canopy_openness", "night_flag")
)


atl08_seg_dt <- rbindlist2(atl08_seg_att_ls)


# Consider only segment with h_canopy < 100 and terrain height < 20000
atl08_seg_dt <- atl08_seg_dt[h_canopy < 100 & h_te_mean < 20000]

head(atl08_seg_dt)

layout(t(1:2))

# ATL03 height histogram
hist(atl03_seg_dt$h_ph, col = "#bd8421", xlab = "Elevation (m)", main = "ATL03 h_ph")
hist(atl08_seg_dt$h_canopy, col = "green", xlab = "Height (m)", main = "ATL08 h_canopy")

library(terra)

blueYellowRed <- function(n) grDevices::hcl.colors(n, "RdYlBu")

set.seed(123)
mask <- base::sample(seq_len(nrow(atl03_seg_dt)), 50)
atl03_seg_vect <- to_vect(atl03_seg_dt)

# Plot with mapview
mapview::mapview(
  atl03_seg_vect[mask],
  zcol = "h_ph",
  layer.name = "h_ph",
  breaks = 3,
  col.regions = blueYellowRed,
  map.types = c("Esri.WorldImagery")
)

# Extract vector from atl08_seg_dt
class(atl08_seg_dt)
atl08_seg_vect <- to_vect(atl08_seg_dt)

# Palette function
greenYellowRed <- function(n) {
grDevices::hcl.colors(n, "RdYlGn")
}


# Plot with mapview
leaflet_available <- require("leaflet")
if (!leaflet_available) stop("leaflet not found!")

map_vect <- mapview::mapView(
atl08_seg_vect,
layer.name = "h_canopy",
zcol = "h_canopy",
col.regions = greenYellowRed,
map.types = c("Esri.WorldImagery")
)

map_vect

terra::writeVector(atl03_seg_vect, file.path(outdir, "atl03_seg.gpkg"))
terra::writeVector(atl08_seg_vect, file.path(outdir, "atl08_seg.gpkg"))

redYellowGreen <- function(n) grDevices::hcl.colors(n, "RdYlGn")
max_h_canopy <- ATL08_seg_attributes_dt_gridStat(atl08_seg_dt, func = max(h_canopy), res = 0.02)

mapview::mapView(
max_h_canopy,
map = map_vect,
col.regions = redYellowGreen
)

multiple_attributes <- ATL08_seg_attributes_dt_gridStat(atl08_seg_dt, func = list(
max_h_canopy = max(h_canopy),
min_h_canopy = min(h_canopy),
mean_canopy_openness = mean(canopy_openness),
mean_h_te_mean = mean(h_te_mean)
), res = 0.01)

map_vect_openness <- mapview::mapView(
atl08_seg_vect,
zcol = "canopy_openness",
layer.name = "canopy_openness",
col.regions = redYellowGreen,
map.types = c("Esri.WorldImagery")
)

blueYellowRed <- function(n) grDevices::hcl.colors(n, "RdYlBu", rev = TRUE)

map_vect_terrain <- mapview::mapView(
atl08_seg_vect,
zcol = "h_te_mean",
layer.name = "h_te_mean",
col.regions = blueYellowRed,
map.types = c("Esri.WorldImagery")
)


m1 <- mapview::mapView(multiple_attributes[[1]], layer.name = "Max h_canopy", map = map_vect, col.regions = redYellowGreen)
m2 <- mapview::mapView(multiple_attributes[[2]], layer.name = "Min h_canopy", map = map_vect, col.regions = redYellowGreen)
m3 <- mapview::mapView(multiple_attributes[[3]], layer.name = "Mean canopy openness", map = map_vect_openness, col.regions = redYellowGreen)
m4 <- mapview::mapView(multiple_attributes[[4]], layer.name = "Mean h_te_mean", col.regions = blueYellowRed, map = map_vect_terrain)

leafsync::sync(m1, m2, m3, m4)

# Define bbox
clip_region <- terra::ext(-83.2, -83.14, 32.12, 32.18)
aoi <- file.path(outdir, "example_aoi.gpkg")
aoi_vect <- terra::vect(aoi)

# Define hdf5 output file
output <- tempfile(pattern = "alt08_h5_clip_", fileext = ".h5")

# Clip the data for only the first ATL08 file
atl08_clipped <- ATL08_h5_clipBox(atl08_h5[[1]], output, clip_obj = clip_region)

##atl08_clippeds <- ATL08_h5_clipGeometry(atl08_h5[[2]], output, clip_obj = aoi_vect, split_by="id")

atl08_seg_dt <- ATL08_seg_attributes_dt(atl08_h5[[1]], attributes = c("h_canopy"))
atl08_seg_dt_clip <- ATL08_seg_attributes_dt(atl08_clipped, attributes = c("h_canopy"))

# Display location of clipped data
atl08_seg_vect <- to_vect(atl08_seg_dt)
atl08_seg_clip_vect <- to_vect(atl08_seg_dt_clip)

bbox <- terra::vect(terra::ext(atl08_seg_clip_vect), crs = "epsg:4326")
centroid <- terra::geom(terra::centroids(bbox))

map1 <- mapview::mapview(
  atl08_seg_clip_vect,
  col.regions = "yellow",
  alpha.regions = 1,
  lwd = 5,
  map.types = c("Esri.WorldImagery"),
  alpha = 0,
  cex = 2,
  legend = FALSE
)

map1

dplyr_available <- require("dplyr")
if (!dplyr_available) stop("dplyr not found!")

# Final map
final_map <- map1@map %>%
  leaflet::addCircleMarkers(data = atl08_seg_vect, radius = 2) %>%
  leaflet::addPolygons(
    data = bbox, fillOpacity = 0, weight = 3, color = "white",
    opacity = 1, dashArray = "5, 1, 0"
  ) %>%
  leaflet::addLegend(
    position = "topright",
    colors = c("blue", "yellow", "white"),
    labels = c("atl08_segment", "atl08_segment_clipped", "bbox"),
    opacity = 1
  ) %>%
  leaflet::setView(centroid[, "x"][[1]], centroid[, "y"][[1]], zoom = 13)

final_map

aoi <- file.path(outdir, "example_aoi.gpkg")
aoi_vect <- terra::vect(aoi)

centroid <- terra::geom(terra::centroids(aoi_vect))

# Extract the h_canopy attribute from the first ATL08 file
atl08_seg_dt <- lapply(atl08_h5, ATL08_seg_attributes_dt, attributes = c("h_canopy"))
atl08_seg_dt <- rbindlist2(atl08_seg_dt)
atl08_seg_vect <- to_vect(atl08_seg_dt)

# Clip using geometry
atl08_seg_dt_clip <- ATL08_seg_attributes_dt_clipGeometry(
atl08_seg_dt, aoi_vect, split_by = "id"
)

atl08_seg_clip_vect <- to_vect(atl08_seg_dt_clip)

colors <- c("#00FF00", "#FF00FF")

map1 <- mapview::mapview(
atl08_seg_clip_vect,
alpha = 0,
col.regions = colors,
alpha.regions = 1,
zcol = "poly_id",
lwd = 5,
map.types = c("Esri.WorldImagery"),
cex = 2,
legend = FALSE
)

# Final map
final_map <- map1@map %>%
leaflet::addCircleMarkers(data = atl08_seg_vect, color = "blue", radius = 2) %>%
leaflet::addPolygons(
  data = aoi_vect, fillOpacity = 0, weight = 3, color = colors,
  opacity = 1, dashArray = "5, 1, 0"
) %>%
leaflet::addLegend(
  position = "topright",
  colors = c("blue", "yellow", colors),
  labels = c("atl08_segment", "atl08_segment_clipped", "aoi_1", "aoi_2"),
  opacity = 1
) %>%
leaflet::setView(centroid[, "x"][[1]], centroid[, "y"][[1]], zoom = 13)

final_map

clipped <- clip(data_object, clip_obj = aoi)

# Extract ATL08 segment attributes
atl08_seg <- ATL08_seg_attributes_dt(atl08_h5[[1]], attributes = "h_canopy")

# Define area of interest
aoi <- file.path(outdir, "example_aoi.gpkg")
aoi_vect <- terra::vect(aoi)

# Clip using the generic function
atl08_seg_clip <- clip(atl08_seg, clip_obj = aoi_vect)

# Convert to vector for visualization
atl08_seg_clip_vect <- to_vect(atl08_seg_clip)

# ATL03 HDF5 object
atl03 <- ATL03_read(atl03_files[[1]])

# Define output path
output <- tempfile(pattern = "atl03_clip_", fileext = ".h5")

# Define bounding box with YOUR coordinates
bbox <- terra::ext(c(-84.71409, -84.6435, 30.12047, 30.17786))

# Clip using generic clip()
atl03_clipped <- clip(atl03, output, clip_obj = bbox)

# import the ATL03 h5 files
atl03_files <- list.files(outdir, "ATL03.*h5", full.names = TRUE)
atl03_h5 <- lapply(atl03_files, ATL03_read)

# Herein as we are working with list of h5 files we will need
# to loop over each file and extract the attributes and then
# concatenate them with rbindlist2
atl03_atl08_dts <- lapply(
seq_along(atl03_h5),
function(ii) {
  ATL03_ATL08_photons_attributes_dt_join(
    atl03_h5[[ii]],
    atl08_h5[[ii]]
  )
}
)

atl03_atl08_dt <- rbindlist2(atl03_atl08_dts)

head(atl03_atl08_dt)

oldpar <- par(no.readonly = TRUE)
par(oma = c(0, 0, 0, 0))
par(mar = c(2, 3, 1, 1))
layout(matrix(c(1, 2), ncol = 1))
plot(
atl03_atl08_dt[orbit_number == 3208],
y = "h_ph",
colors = c("gray", "#bd8421", "forestgreen", "green"),
xlim = c(25500, 28500),
beam = "gt2r",
cex = 0.5,
pch = 16
)

par(mar = c(3, 3, 1, 1))

plot(
atl03_atl08_dt[orbit_number == 3208],
y = "ph_h",
colors = c("gray", "#bd8421", "forestgreen", "green"),
xlim = c(25500, 28500),
beam = "gt2r",
cex = 0.5,
pch = 16,
legend = FALSE
)

par(
oldpar
)

h_canopy <- ATL03_ATL08_photons_attributes_dt_gridStat(
atl03_atl08_dt[ph_h < 50 & ph_h > 0],
func = list(
  h_canopy = quantile(ph_h, 0.98),
  count = .N
),
res = 0.01
)

plot(h_canopy,
   col = viridis::inferno(100),
   xlab = "Langitude (degree)",
   ylab = "Latitude (degree)",
   ylim = c(32.1, 32.4)
)

# Herein as we are working with list of h5 files we will need
# to loop over each file and extract the attributes and then
# concatenate them with rbindlist2

stopifnot(length(atl03_h5) == length(atl08_h5))

atl03_atl08_dts <- lapply(
  seq_along(atl03_h5),
  function(ii) {
    ATL03_ATL08_photons_attributes_dt_join(
      atl03_h5[[ii]],
      atl08_h5[[ii]]
    )
  }
)

atl03_atl08_dt <- rbindlist2(atl03_atl08_dts)
head(atl03_atl08_dt)

atl03_atl08_photons_grouped_dt <- ATL03_ATL08_segment_create(atl03_atl08_dt,
                                                             segment_length = 30,
                                                             centroid = "mean",
                                                             output = NA,
                                                             overwrite = FALSE
)

atl03_atl08_seg_dt <- ATL03_ATL08_compute_seg_attributes_dt_segStat(
  atl03_atl08_photons_grouped_dt,
  list(
    h_canopy_ge0 = quantile(ph_h, 0.98),
    h_canopy_gt0 = quantile(ph_h[ph_h > 0], 0.98),
    n_ground = sum(classed_pc_flag == 1),
    n_mid_canopy = sum(classed_pc_flag == 2),
    n_top_canopy = sum(classed_pc_flag == 3),
    n_canopy_total = sum(classed_pc_flag >= 2)
  ),
  ph_class = c(1, 2, 3)
)

head(atl03_atl08_seg_dt)

atl03_atl08_vect <- to_vect(atl03_atl08_seg_dt[h_canopy_gt0 <= 31])

# Get coordinates of segment 1358
seg_1358 <- atl03_atl08_seg_dt[segment_id == 1358]
seg_1358$longitude
seg_1358$latitude

# Use those coordinates as centroid
map_output <- mapview::mapview(
  atl03_atl08_vect,
  zcol        = "h_canopy_gt0",
  col.regions = grDevices::hcl.colors(9, "RdYlGn"),
  alpha       = 0,
  layer.name  = "h_canopy",
  map.types   = c("Esri.WorldImagery"),
  cex         = 4
)@map %>%
  leaflet::setView(
    lng  = seg_1358$longitude, 
    lat  = seg_1358$latitude, 
    zoom = 10
  )

map_output

# For the sake of the example, we will train and test the model with the same data
library(randomForest)

h_canopy <- c(
  13.9, 3.1, 2.2, 4.6, 21.6,
  7.2, 5, 7.7, 0.8, 9.7,
  11, 11.3, 15.5, 5.1, 10.4,
  0.6, 14.6, 13.3, 9.8, 14.7
)

agbd <- c(
  144.8, 27.5, 51.6, 60.5, 232.3,
  102.8, 33.1, 91.3, 23, 120.1,
  125.7, 127.2, 147.4, 48.8, 103.3,
  55.9, 181.8, 139.9, 120.1, 162.8
)

set.seed(172783946)
model <- randomForest::randomForest(data.frame(h_canopy = h_canopy), agbd)

out_h5 <- tempfile(fileext = ".h5")

for (atl08_h5_item in atl08_h5) {
  atl08_seg_dt <- ATL08_seg_attributes_dt(atl08_h5_item, attributes = c("h_canopy"))
  atl08_seg_dt[h_canopy > 100, h_canopy := NA_real_]
  atl08_seg_dt <- na.omit(atl08_seg_dt)
  predicted_h5 <- predict_h5(model, atl08_seg_dt, out_h5)
}

# Fix PROJ path
Sys.setenv(PROJ_DATA = system.file("proj", package = "ICESat2VegR"))
Sys.setenv(PROJ_LIB  = system.file("proj", package = "ICESat2VegR"))

# Verify
system.file("proj", package = "ICESat2VegR")  # should return a path

# Get the longitude/latitude of the predicted points
x <- predicted_h5[["longitude"]][]
y <- predicted_h5[["latitude"]][]

# Use only your AOI coordinates
x_clip <- x[x >= -84.72 & x <= -84.64]
y_clip <- y[y >= 30.12  & y <= 30.18]

bbox <- terra::ext(
  min(x_clip), max(x_clip),
  min(y_clip), max(y_clip)
)
bbox

# Rasterize using the clipped AOI
output_raster <- tempfile(fileext = ".tif")
rasterize_h5(predicted_h5, output_raster, bbox = bbox, res = 0.005)

# Or rasterize the full extent of the predicted points instead
output_raster <- tempfile(fileext = ".tif")
bbox <- terra::ext(min(x), max(x), min(y), max(y))

# Creates the raster with statistics
res <- 0.02
rasterize_h5(predicted_h5, output_raster, bbox = bbox, res = res)

# Open the raster by file path
forest_height_palette <- c("#ffffff", "#4d994d", "#004d00")

# Open band 2 only (mean AGBD)
library(leaflet)

stars_rast <- stars::read_stars(output_raster, RasterIO = list(bands = 2))
res_map <- mapview::mapview(
  stars_rast,
  layer.name = "AGBD mean",
  col.regions = forest_height_palette,
  na.alpha = 0.1,
  map = leaflet::leaflet() %>% leaflet::addProviderTiles("Esri.WorldImagery")
)
res_map

# The r-universe version (recommended for the latest version)
install.packages("ICESat2VegR", , repos = c("https://caiohamamura.r-universe.dev", "https://cloud.r-project.org"))
#install.packages('ICESat2VegR', repos = c('https://carlos-alberto-silva.r-universe.dev', 'https://cloud.r-project.org'))

# The CRAN version
#install.packages("ICESat2VegR")

# Required additional libraries
need_pkgs <- c(
  "reticulate",   # Python <-> R interface
  "leaflet",      # interactive maps
  "sf",           # spatial vectors
  "terra",        # rasters & vectors
  "data.table",   # fast tables
  "dplyr",        # tidy helpers
  "mapview",      # interactive map viewer
  "caret",        # modeling / ML
  "hdf5r"         # reading/writing local HDF5 (.h5) files
)
missing <- need_pkgs[!need_pkgs %in% rownames(installed.packages())]
if (length(missing)) {
  message("Installing missing R packages: ", paste(missing, collapse = ", "))
  install.packages(missing, repos = repos, dependencies = TRUE)
}

suppressPackageStartupMessages({
  library(ICESat2VegR)
  library(reticulate)
  library(leaflet)
  library(sf)
  library(terra)
  library(data.table)
  library(dplyr)
  library(mapview)
  library(caret)
})
cat("Loading libraries... done.\n")                                                    # Confirm loading

aoi_path <- system.file("extdata", "aoi_4326.geojson", package = "ICESat2VegR")
boundary <- sf::st_read(aoi_path, quiet = TRUE)
box      <- sf::st_bbox(boundary)

lower_left_lon  <- box["xmin"];  lower_left_lat  <- box["ymin"]
upper_right_lon <- box["xmax"];  upper_right_lat <- box["ymax"]

daterange <- c("2025-08-01", "2025-08-31")

atl03_granules_cloud <- ICESat2VegR::ATLAS_dataFinder(
  short_name = "ATL03",
  lower_left_lon, lower_left_lat, upper_right_lon, upper_right_lat,
  version = "007", daterange = daterange, persist = TRUE, cloud_computing = FALSE
)

atl08_granules_cloud <- ICESat2VegR::ATLAS_dataFinder(
  short_name = "ATL08",
  lower_left_lon, lower_left_lat, upper_right_lon, upper_right_lat,
  version = "007", daterange = daterange, persist = TRUE, cloud_computing = FALSE
)

outdir <- tempdir()
ATLAS_dataDownload(atl03_granules_cloud, outdir)
ATLAS_dataDownload(atl08_granules_cloud, outdir)

atl03_files <- list.files(outdir, pattern = "ATL03.*h5", full.names = TRUE)
atl08_files <- list.files(outdir, pattern = "ATL08.*h5", full.names = TRUE)
get_timestamp <- function(f) sub(".*_(\\d{14})_.*", "\\1", basename(f))
timestamps_comuni <- intersect(sapply(atl08_files, get_timestamp),
                               sapply(atl03_files, get_timestamp))

ts    <- timestamps_comuni[1]
atl03 <- ICESat2VegR::ATL03_read(atl03_files[grepl(ts, atl03_files)])
atl08 <- ICESat2VegR::ATL08_read(atl08_files[grepl(ts, atl08_files)])

dt <- ICESat2VegR::ATL03_ATL08_photons_attributes_dt_join(atl03, atl08)
head(dt)

seg20 <- ATL03_ATL08_segment_create(dt, 20, centroid = "mean", output = NA, overwrite = FALSE)

stats20 <- ATL03_ATL08_compute_seg_attributes_dt_segStat(
  seg20,
  list(
    rh98           = quantile(ph_h[ph_h > 0 & classed_pc_flag %in% c(2,3)], 0.98),
    n_ground       = sum(classed_pc_flag == 1),
    n_top_canopy   = sum(classed_pc_flag == 3),
    n_canopy_total = sum(classed_pc_flag >= 2),
    mean_solar     = mean(solar_elevation),
    night_flag2    = as.integer(mean(night_flag) > 0.5)
  ),
  ph_class = c(2, 3)
)
# Convert to SpatVector first
stats20_vect <- to_vect(stats20)

centroid <- stats20[, .(x = mean(longitude), y = mean(latitude))]

map_output <- mapview::mapview(
  stats20_vect,
  zcol        = "rh98",
  col.regions = grDevices::hcl.colors(9, "RdYlGn"),
  alpha       = 0,
  layer.name  = "rh98 (m)",
  map.types   = c("Esri.WorldImagery"),
  cex         = 4
)@map %>%
  leaflet::setView(lng = centroid$x, lat = centroid$y, zoom = 10)

map_output

stats20_clip <- ATL03_ATL08_seg_attributes_dt_clipGeometry(stats20, boundary)
stats20_clip$year <- as.integer(substr(ts, 1, 4))
stats20_clip <- stats20_clip[stats20_clip$rh98 <= 50]

stats20_clip_sf <- st_as_sf(
as.data.frame(stats20_clip),
coords = c("longitude", "latitude"), crs = 4326, remove = FALSE
)
st_write(stats20_clip_sf,
       file.path(outdir, paste0("ATL03_ATL08_", ts, "_20.geojson")),
       delete_dsn = TRUE, quiet = TRUE)

Sys.setenv(EE_PROJECT = "your-ee-project")
ICESat2VegR::ee_initialize()

boundary_simple <- sf::st_simplify(
  sf::st_union(sf::st_make_valid(boundary)),
  dTolerance = 0.001
)
region <- .as_ee_geom(boundary_simple)

seg_path <- system.file("extdata", "ATL03_ATL08_example_segments.geojson", package = "ICESat2VegR")
data_raw <- sf::read_sf(seg_path)
df_dt    <- data.table::as.data.table(sf::st_drop_geometry(data_raw))
class(df_dt) <- c("icesat2.atl03_atl08_seg_dt", "data.table", "data.frame")

set.seed(1)
df_sampled <- ICESat2VegR::sample(df_dt, method = randomSampling(500))

df_vect <- terra::vect(
  as.data.frame(df_sampled), geom = c("longitude", "latitude"), crs = "EPSG:4326"
)

yr          <- 2025
predictors2 <- ee_build_AlphaEarth_embedding_terrain_stack(region, yr, yr)
sampling_df <- ICESat2VegR::seg_ancillary_extract(predictors2, df_vect, 20)
head(sampling_df[, 1:12])

predictors2_rgb <- predictors2$select(c("A00", "A20", "A40"))
#
centroid_lon <- mean(terra::ext(terra::vect(boundary))[c(1,2)])
centroid_lat <- mean(terra::ext(terra::vect(boundary))[c(3,4)])
#
leaflet::leaflet() |>
addEEImage(
  predictors2_rgb,
  bands = c("A00", "A20", "A40"),
  group = "RGB Embedding",
  min = c(-0.06),
  max = c( 0.12)
) |>
leaflet::addControl(
  html     = "<b>RGB Embedding</b><br>Bands: A00 / A20 / A40",
  position = "bottomleft"
) |>
leaflet::setView(lng = centroid_lon, lat = centroid_lat, zoom = 16) |>
leaflet::addLayersControl(
  overlayGroups = "RGB Embedding",
  options = leaflet::layersControlOptions(collapsed = FALSE)
)

x <- sampling_df %>%
  dplyr::select(starts_with("A"), "slope", "aspect", "elevation") %>%
  data.frame()
y <- sampling_df %>% dplyr::select("rh98") %>% data.frame()

sel_rf_rfe   <- varSel(x, y$rh98)
best_metrics <- sel_rf_rfe$selvars
print(best_metrics)

# Build importance tables from sel_rf_rfe
imp_desc <- sel_rf_rfe$importance %>%
  dplyr::mutate(selected = parameter %in% sel_rf_rfe$selvars) %>%
  dplyr::arrange(importance)

best_imp_rfe <- imp_desc %>%
  dplyr::filter(selected == TRUE) %>%
  dplyr::arrange(importance)

par(mfrow = c(1, 2), mar = c(4, 6, 2, 1))

# Left — all RFE evaluated, colored by selection
barplot(
  rev(imp_desc$importance),
  names.arg = rev(as.character(imp_desc$parameter)),
  horiz     = TRUE,
  las       = 1,
  main      = "RFE importance\n(green = selected)",
  xlab      = "Importance",
  col       = ifelse(rev(imp_desc$selected), "#1B7837", "grey70"),
  cex.names = 0.7
)
abline(v = 0.2, lty = 2, col = "red")

# Right — selected only
barplot(
  best_imp_rfe$importance,
  names.arg = as.character(best_imp_rfe$parameter),
  horiz     = TRUE,
  las       = 1,
  main      = "Selected predictors",
  xlab      = "Importance",
  col       = "#1B7837",
  cex.names = 0.8
)
abline(v = 0.2, lty = 2, col = "red")

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2))

x      <- dplyr::select(sampling_df, dplyr::any_of(best_metrics))
y      <- sampling_df$rh98
data_i <- data.frame(y, x, check.names = FALSE)

idx_train <- caret::createDataPartition(y = data_i$y, p = 0.70, list = FALSE)
trainData <- data_i[ idx_train, ]
testData  <- data_i[-idx_train, ]

fit_rf   <- fit_model(x = dplyr::select(trainData, dplyr::any_of(best_metrics)),
                      y = trainData$y,
                      rf_args = list(ntree = 100))
rf_model <- fit_rf$model

pred_train <- predict(rf_model, newdata = dplyr::select(trainData, dplyr::any_of(best_metrics)))
pred_test  <- predict(rf_model, newdata = dplyr::select(testData,  dplyr::any_of(best_metrics)))

results_train <- fit_metrics(trainData$y, as.numeric(pred_train))
results_test  <- fit_metrics(testData$y,  as.numeric(pred_test))

head(results_train)
head(results_test)

mosaic <- ee_build_AlphaEarth_embedding_terrain_stack(region, yr, yr)

ch_map <- map_create(
model    = fit_rf$model,
stack    = mosaic,
aoi      = region,
reducer  = "mosaic",
mode     = "auto",
to_float = TRUE
)

min_hcanopy <- 0
max_hcanopy <- 40

forest_palette <- colorRampPalette(
  c("#fcd88f", "#00441B", "#1B7837", "#A6DBA0", "#E7E1EF", "#762A83")
)(128)

pal_fun <- leaflet::colorNumeric(forest_palette, domain = c(min_hcanopy, max_hcanopy))

leaflet::leaflet() |>
  addEEImage(
    ch_map,
    bands   = "prediction_layer",
    group   = "Canopy Height",
    min     = min_hcanopy,
    max     = max_hcanopy,
    palette = forest_palette
  ) |>
  leaflet::addLegend(
    pal      = pal_fun,
    values   = c(min_hcanopy, max_hcanopy),
    opacity  = 1,
    title    = "Canopy Height (m)",
    position = "bottomleft"
  ) |>
  leaflet::setView(lng = centroid_lon, lat = centroid_lat, zoom = 15) |>
  leaflet::addLayersControl(
    overlayGroups = "Canopy Height",
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )

googledrive::drive_auth(cache = FALSE)

out <- map_download(
ee_image         = ch_map,
method           = "drive",
region           = region,
scale            = 10,
file_name_prefix = "prediction_2025",
dsn              = file.path(outdir, "prediction_2025.tif"),
drive_folder     = "EE_Exports",
monitor          = TRUE
)
terra::plot(out)

lapply(atl03_h5, close)
lapply(atl08_h5, close)