-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathDeepBioFusion_Preprocessing.R
More file actions
105 lines (88 loc) · 3.98 KB
/
DeepBioFusion_Preprocessing.R
File metadata and controls
105 lines (88 loc) · 3.98 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
# Load required libraries
library(raster)
library(rgdal)
library(sp)
library(lidR)
library(viridis)
# Set working directory and load LiDAR data
setwd("your/working/directory")
las <- readLAS('AOI.laz')
if (is.empty(las)) stop("LAS file is empty or corrupted.")
# Normalize heights and generate DTM/CHM
dtm <- grid_terrain(las, res = 0.2, algorithm = tin())
las <- normalize_height(las, dtm)
chm <- grid_canopy(las, res = 0.25, p2r(0.3))
# Plot and export DTM and CHM side-by-side
lidar_colors <- colorRampPalette(c("blue", "green", "yellow", "brown", "white"))(100)
pdf("outputs/dtm_chm_lidar_colors.pdf", width = 12, height = 6)
par(mfrow = c(1, 2))
plot(dtm, col = lidar_colors, main = "Digital Terrain Model (DTM)")
plot(chm, col = lidar_colors, main = "Canopy Height Model (CHM)")
dev.off()
# Segment trees
las <- segment_trees(las, li2012(), attribute = "treeID")
ttops <- find_trees(chm, lmf(4, hmin = 2))
# Extract convex crown hulls
crowns <- delineate_crowns(las, type = "convex")
plot(crowns)
# Load ancillary raster datasets
raster_sar <- raster("AOI_SAR.tif")
raster_optical <- raster("AOI_Optical.tif")
raster_species <- crop(raster("AOI_TreeSpecies.tif"), extent(las))
# Extract species class at treetop locations
tree_coords <- coordinates(ttops)
tree_species <- extract(raster_species, tree_coords)
# Replace '0' species values by neighbors (simple forward-backward fill)
for (i in 2:(length(tree_species)-1)) {
if (tree_species[i] == 0) {
tree_species[i] <- ifelse(tree_species[i-1] != 0, tree_species[i-1], tree_species[i+1])
}
}
# Define species-specific parameters
params <- list(
spruce = list(b1=2.2131, b2=0.3046, lna=-0.5244, k=1.013, b=8.8563, m=19, c=0, d=0.3879, asym=36),
pine = list(b1=2.2845, b2=0.3318, lna=-1.448, k=1.009, b=8.7399, m=16, c=0, d=0.5624, asym=28),
deciduous = list(b1=1.649, b2=0.373, lna=-2.1284, k=1.004, b=9.3375, m=11, c=0.0221, d=0.2838, asym=20)
)
# Helper functions
calculate_DBH <- function(species, Z) {
if (species == 0) return(0)
p <- params[[species]]
Z <- min(Z, p$asym)
return(((Z - 1.3)^(1/3) * p$b1) / (1 - p$b2 * (Z - 1.3)^(1/3)))
}
calculate_Biomass <- function(species, Z, DBH) {
if (species == 0) return(0)
p <- params[[species]]
Z <- min(Z, p$asym)
return(p$k * exp(p$lna + p$b * (DBH/(DBH + p$m)) + p$c * Z + p$d * log(Z)))
}
# Compute DBH and biomass
DBH <- numeric(length(ttops))
Biomass <- numeric(length(ttops))
for (i in seq_along(DBH)) {
sp <- switch(tree_species[i] + 1, "none", "spruce", "pine", "deciduous")
Z <- ttops$Z[i]
DBH[i] <- calculate_DBH(sp, Z)
Biomass[i] <- calculate_Biomass(sp, Z, DBH[i])
}
# Export results to CSV
tree_data <- data.frame(treeID = ttops$treeID, Z = ttops$Z, Species = tree_species, X = tree_coords[,1], Y = tree_coords[,2], DBH, Biomass)
write.csv(tree_data, "outputs/tree_data_with_biomass.csv", row.names = FALSE)
# -----------------------------------------------------------
# OPTIONAL: Generate 16x16 patches around each tree for training
# -----------------------------------------------------------
clip_tree_patches <- function(raster_data, coords, output_dir, size = 16) {
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
for (i in 1:nrow(coords)) {
ext <- extent(coords[i, 1] - size/2, coords[i, 1] + size/2,
coords[i, 2] - size/2, coords[i, 2] + size/2)
clip <- crop(raster_data, ext)
writeRaster(clip, filename = file.path(output_dir, paste0("tree_", i, ".tif")), format = "GTiff", overwrite = TRUE)
}
}
clip_tree_patches(raster_optical, tree_coords, "outputs/patches_optical16x16")
clip_tree_patches(raster_sar, tree_coords, "outputs/patches_sar16x16")
clip_tree_patches(raster("AOI_SARL.tif"), tree_coords, "outputs/patches_sarL_32x32", size = 32)
clip_tree_patches(raster("SARSentinel.tif"), tree_coords, "outputs/patches_sarC_32x32", size = 32)
clip_tree_patches(raster("SARX10.tif"), tree_coords, "outputs/patches_sarX10_32x32", size = 32)