CHM-Based Fuel Downscaling

Author

Mickey Campbell

Purpose

I just want to make sure I’m doing this right :). So, Nuria, I’d appreciate if you could take a look and let me know before I adapt this approach into the voxelFuel package.

Step 1: Setup

I’ll note that the two inputs from this example analysis came originally from Nuria, both representing Sycan marsh, one being a 10m foliage biomass map, and the other being a 1m CHM.

library(terra)

Warning: package 'terra' was built under R version 4.5.3

terra 1.9.1

# directories
main_dir <- "S:/ursa/campbell/fb_downscaling/data/from_nuria/FB_Downscaling"
in_dir <- file.path(main_dir, "input")
out_dir <- file.path(main_dir, "output")

# read in the data
fb_10m <- file.path(in_dir, "FB_Sycan.asc") |> rast()
chm_1m <- file.path(in_dir, "CHM_1m_Sycan.tif") |> rast()

# fix projection
crs(fb_10m) <- crs(chm_1m)

# clean names
names(fb_10m) <- "fb"
names(chm_1m) <- "chm"

# plot them out just to make sure they look ok
plot(fb_10m)

plot(chm_1m)

Step 2: Rescaling

Not that it matters hugely for the sake of example, but I believe the foliage biomass map is in kg/m². So, this step is designed to get total foliage biomass (in kg) per pixel.

# convert area-normalized biomass to total biomass per pixel and plot
fb_10m <- fb_10m * 10^2
plot(fb_10m)

Step 3: The fun part

Here’s where things get interesting. The logic here is:

Get “total height” (i.e., CHM pixel sum) within 10m superpixels. Result: chm_10m_sum.
Use nearest neighbor downsampling of this 10m total height raster to get back down to 1m. Result: chm_1m_sum.
Get the proportion of superpixel total height each original CHM pixel height represents. Result: chm_1m_prop.
Use nearest neighbor downsampling of the 10m foliage biomass raster to get down to 1m. Result: fb_1m_sub.
Multiply fine-scale proportional height raster with fine-scale foliage biomass raster to get apportioned foliage biomass. Result: fb_1m.

# aggregate chm to get total height within 10m chm superpixels
chm_10m_sum <- aggregate(chm_1m, 10, "sum")


|---------|---------|---------|---------|
=========================================

# disaggregate 10m chm superpixels to get subsampled fine-scale height totals
chm_1m_sum <- disagg(chm_10m_sum, 10)

# divide out to get fine-scale proportional height within superpixel
chm_1m_prop <- chm_1m / chm_1m_sum
chm_1m_prop <- ifel(chm_1m_sum == 0, 0, chm_1m_prop)

# disaggregate 10m biomass superpixels to get subsampled fine-scale biomass
fb_1m_sub <- disagg(fb_10m, 10)

# multiply out to apportion fine-scale biomass according to height proportion
fb_1m <- fb_1m_sub * chm_1m_prop

Step 4: Check your work

Just to double-check, aggregating fb_1m back up to 10m with a sum function should yield the same map as the original input.

# reaggregate for checking purposes
fb_10m_test <- aggregate(fb_1m, 10, "sum")


|---------|---------|---------|---------|
=========================================

# stack and plot
fb_10m_stack <- c(fb_10m, fb_10m_test)
names(fb_10m_stack) <- c("fb_orig", "fb_test")
plot(fb_10m_stack)

# plot the difference (should be all zeros)
plot(fb_10m_test - fb_10m)

Step 5: Build the model

Here I’ll combine the 1m CHM and 1m foliage biomass, take a subsample of 10k random non-NA pixels, and plot the data out. In this case, a linear model looks most appropriate (I know you were using an exponential model, but for this dataset, linear looks best…?).

# stack rasters
fb_chm_1m <- c(fb_1m, chm_1m)
names(fb_chm_1m) <- c("fb", "chm")

# randomly sample
samp <- spatSample(
  x = fb_chm_1m,
  size = 10000,
  na.rm = T,
  xy = T
)

# plot the data out
par(mar = c(4.5,4.5,0.5,0.5), las = 1)
plot(
  fb ~ chm, 
  data = samp, 
  xlab = "Canopy Height (m)", 
  ylab = "Foliage Biomass (kg)"
)

# build model
mod <- lm(fb ~ chm + 0, data = samp)

# add regression line and stats
abline(mod, lwd = 2, col = 2)
r2 <- summary(mod)$adj.r.squared |> formatC(2, format = "f")
r2 <- bquote(R^2==.(r2))
legend(
  "topleft",
  legend = r2, 
  bty = "n", 
  text.col = 2, 
  text.font = 2,
  x.intersp = 0
)

Step 6: Predict with the model and evaluate

The next step is to apply this model to the 1m CHM to generate a predicted foliage biomass map at 1m. To evaluate it, I will re-aggregate to 10m, and see how it compares to the original 10m biomass map, first just visually…

# predict at 1m
fb_1m_pred <- predict(chm_1m, mod)

# aggregate up to 10m
fb_10m_pred <- aggregate(fb_1m_pred, 10, "sum")


|---------|---------|---------|---------|
=========================================

# stack back with original map
fb_10m_comp <- c(fb_10m, fb_10m_pred)
names(fb_10m_comp) <- c("fb_orig", "fb_pred")

# map them out
max_fb <- global(fb_10m_comp, "max", na.rm = T) |> max()
plot(fb_10m_comp, range = c(0, max_fb))

Then, more quantitatively, with a 10k sample of pixels…

# plot pixel-wise agreement
fb_10m_comp_df <- as.data.frame(fb_10m_comp, na.rm = T)
fb_10m_comp_df_samp <- fb_10m_comp_df[sample(1:nrow(fb_10m_comp_df), 10000),]
par(mar = c(4.5,4.5,0.5,0.5), las = 1)
plot(
  fb_pred ~ fb_orig, 
  data = fb_10m_comp_df_samp,
  xlim = c(0, max_fb), 
  ylim = c(0, max_fb),
  xlab = "Original 10m Biomass",
  ylab = "Predicted 10m Biomass"
)
abline(0,1)
pred_obs_mod <- lm(fb_pred ~ fb_orig, data = fb_10m_comp_df_samp)
abline(pred_obs_mod, lwd = 2, col = 2)
r2 <- summary(mod)$adj.r.squared |> formatC(2, format = "f")
r2 <- bquote(R^2==.(r2))
legend(
  "topleft",
  legend = r2, 
  bty = "n", 
  text.col = 2, 
  text.font = 2,
  x.intersp = 0
)

And lastly, with histograms…

# plot out their histograms
par(mfrow = c(1,2), mar = c(4.5,4.5,1.5,0.5), las = 1)
hist(fb_10m_comp[[1]], xlim = c(0, max_fb))

Warning: [hist] a sample of 74% of the cells was used (of which 0% was NA)

hist(fb_10m_comp[[2]], xlim = c(0, max_fb))

Warning: [hist] a sample of 74% of the cells was used (of which 0% was NA)

Next Steps

Assuming this looks about right, the next step would then be to take these 1m foliage biomass pixels, and apportion them out vertically to voxels according to the proportional vertical distribution of points in each voxel column.