To test out the moving window approach to downscaling fuel load maps versus the regression-aided approach in a 2D, canopy height model-based context.
I am going to use two inputs:
OR_Sycan_2021_cfl_kg_m2_proj.tif. It has a spatial resolution of 30m and the filename suggests it represents CFL in kg/m2.USGS_LPC_OR_SouthwestCentralSycan_2021_B21). The tile I am using is ...10TFN4948.laz.# load libraries
library(lidR)
library(terra)
# mute progress bars
options(lidR.progress = F)
# read in the lidar data, with some filters
las <- readLAS(
"S:/ursa/campbell/fb_downscaling/data/als/als_2021/a02_hgt/USGS_LPC_OR_SouthwestCentralSycan_2021_B21_10TFN4948.laz",
filter = "-drop_class 7 18 -drop_withheld -keep_z -1 50"
)
# read in the fuel data
fuel <- rast("S:/ursa/campbell/fb_downscaling/data/from_nuria/CFL_downscaling/input/OR_Sycan_2021_cfl_kg_m2_proj.tif")# plot the fuel map at full scale
plot(
fuel,
plg = list(
title = expression("Fuel load (kg m"^{-2}*")")
)
)
# add the lidar tile center point
x_center <- mean(range(las$X, na.rm = TRUE))
y_center <- mean(range(las$Y, na.rm = TRUE))
pt_center <- vect(
data.frame(x = x_center, y = y_center),
geom = c("x", "y"),
crs = projection(las)
) |>
project(fuel)
points(
pt_center,
pch = 21,
cex = 1.5,
bg = "orange",
col = "white",
lwd = 1.5
)
add_legend(
"topright",
legend = "Lidar tile",
pch = 21,
pt.cex = 1.5,
pt.bg = "orange",
col = "white",
pt.lwd = 1.5
)
# get lidar tile extent
bb <- st_bbox(las)
e <- ext(bb[1], bb[3], bb[2], bb[4]) |>
as.polygons() |>
set.crs(projection(las))
e_buff <- e |>
project(fuel) |>
buffer(res(fuel)[1]) |>
ext()
# crop fuel raster to extent
fuel_coarse <- crop(fuel, e_buff, touches = T, snap = "out")
# plot it out
plot(e)
plot(
fuel_coarse,
plg = list(
title = expression("Fuel load (kg m"^{-2}*")")
),
add = T,
)
lines(as.polygons(e), lwd = 2, col = "orange")
add_legend(
"topright",
legend = "Lidar tile",
pch = 0,
pt.cex = 2,
col = "orange",
pt.lwd = 2
)
In this step, I will convert the original map from biomass density (kg/m2) to biomass (kg).
# convert biomass density to biomass
fuel_coarse <- fuel_coarse *
res(fuel_coarse)[1] *
res(fuel_coarse)[2]
# plot it out
plot(
fuel_coarse,
plg = list(
title = "Fuel load (kg)"
)
)
By default, I’ll use the pitfree() algorithm built into lidR as it tends to produce nice-looking CHMs.
# interpolate chm
chm <- rasterize_canopy(las, res = 1, algorithm = pitfree())
# set pixels below 2m to 0 to isolate "canopy"
chm <- ifel(chm < 2, 0, chm)
# get chm boundary
chm_boundary <- as.polygons(ext(chm))
# plot chm
plot(
chm,
plg = list(
title = "Canopy Height (m)"
)
)
Here is where we begin the downscaling process. The first step in this process is to get the proportional height of each CHM pixel relative to the total CHM height in a focal neighborhood (moving window) around it.
For the sake of simplicity, I will start with a square kernel. Note that there is no way to get a square kernel that will exactly match the ratio of spatial resolutions in this case (30:1), since kernels are required to have odd dimensions. So, my solution is to create a kernel that is one pixel larger than the original ratio.
Using this moving window, I will calculate the focal sum of canopy heights:
# get coarse/fine resolution ratio as integer multiplier
res_mult <- as.integer(res(fuel_coarse) / res(chm))[1]
# get window size
w <- ifelse(res_mult %% 2 == 0, res_mult + 1, res_mult)
# get chm focal sum
chm_focal_sum <- focal(
x = chm,
w = w,
fun = "sum",
na.rm = T,
na.policy = "omit"
)
# plot it out
plot(
chm_focal_sum,
plg = list(
title = "Canopy Height Focal Sum (m)"
)
)
I will then divide individual pixel-level canopy heights to the focal sum to get focal proportional height. That is, each pixel will capture the proportion of canopy height relative to the sum thereof within its 31x31m surrounding area.
# get chm focal prop
chm_focal_prop <- chm / chm_focal_sum
# fix division by zero for cases where there are no pixels greater than 2m
chm_focal_prop <- ifel(chm_focal_sum == 0, 0, chm_focal_prop)
# plot it out
plot(
chm_focal_prop,
plg = list(
title = "Canopy Height Focal Proportion"
)
)
The next step is to use that exact same moving window kernel to get focal mean fuel load.
To do this, I first have to resample the original 30m fuel map down to 1m to match the resolution of the CHM. To ensure that they are in the same extent, CRS, and pixel grid, I will use the project() function with a nearest neighbor resampling algorithm. In effect, this should look exactly the same as Figure 3, but under the hood, there are 900 1m pixels with replicated values for each original 30m pixel.
# reproject and crop fuel rast to match chm
fuel_fine_near <- project(fuel_coarse, chm, method = "near")
# plot it out
plot(
fuel_fine_near,
plg = list(
title = "Fuel Load (kg)"
)
)
Next, I will get focal mean fuel load. In effect, this will just look like a smoothed out version of the previous figure.
# get fuel focal mean
fuel_focal_mean <- focal(
x = fuel_fine_near,
w = w,
fun = "mean",
na.rm = T,
na.policy = "omit"
)
# plot it out
plot(
fuel_focal_mean,
plg = list(
title = "Fuel Load (kg)"
)
)
As a double-check, the global means of these two maps (the nearest-neighbor 1m-resampled fuel map in Figure 7 and the focal mean-smoothed map in Figure 8) should be identical:
Within a margin of rounding error. Let’s move on…
The final(-ish) step in this process is to multiply proportional canopy height by focal mean fuel load. The key assumption here is that canopy height explains fuel load. Therefore, within a neighborhood surrounding each cell, the proportion of canopy height that that pixel possesses should equate to the proportion of fuel in that same neighborhood.
# get height-apportioned fuel
fuel_fine <- fuel_focal_mean * chm_focal_prop
# plot it out
plot(
fuel_fine,
plg = list(
title = "Fuel Load (kg)"
)
)
If you look closely, you will notice some “spikes” in fuel load. These are found with spatially isolated trees. This makes sense – for an extreme example, if there is only one tree found within a 31x31m area surrounding area, then all of that original biomass has to be assigned to that one tree. Potentially problematic?
One of the theoretical benefits of this approach is the fact that, when you re-aggregate the fine-scale map, it should approximately match the original, coarse-scale pixel values. I’ll test that theory here:
# reaggregate for testing
fuel_reagg <- resample(fuel_fine, fuel_coarse, method = "sum") |> crop(e)
# crop original fuel map to lidar extent
fuel_orig <- crop(fuel_coarse, e)
# stack them up
fuel_stack <- c(fuel_orig, fuel_reagg)
names(fuel_stack) <- c("original", "reaggregated")
# get global maximum for consistent visualization
fuel_max <- global(fuel_stack, "max") |> max()
# plot them out
plot(fuel_stack, range = c(0, fuel_max))
Interestingly, the re-aggregated results are noisier than the original map. The patterns are very close, but the re-aggregated results tend to yield more extreme highs and lows than the original map. This is also evident in a pixel-to-pixel scatterplot comparison:
# convert pixels to data.frame
df <- as.data.frame(fuel_stack)
x <- df$original
y <- df$reaggregated
# plot it out
plot(
x = x,
y = y,
pch = 16,
col = rgb(0,0,0,0.1),
xlim = range(df),
ylim = range(df),
xlab = "original",
ylab = "reaggregated",
las = 1
)
grid()
abline(0,1)
mod <- lm(y ~ x)
abline(mod, lwd = 2, col = 2)
r2 <- summary(mod)$r.squared |> formatC(2, format = "f")
rmse <- sqrt(mean((x - y)^2)) |> formatC(2, format = "f")
bias <- mean(y - x) |> formatC(2, format = "f")
r2_txt <- bquote(bold(R^2==.(r2)))
rmse_txt <- bquote(bold(RMSE==.(rmse)))
bias_txt <- bquote(bold(Bias==.(bias)))
legend(
"topleft",
legend = c(r2_txt, rmse_txt, bias_txt),
text.col = 2,
bty = "n",
x.intersp = 0
)
Let’s see if the added regression step could help improve invertibility, by potentially mitigating some of the extreme values found inherent to the moving window downscaling approach. To do this I’ll first explore the relationship between the 1m downscaled fuel map and the 1m CHM. I’ll fit a negative exponential model to it.
# generate spatial sample
samp <- spatSample(
x = fuel_fine,
size = 10000,
as.points = T
)
# extract chm values
samp <- extract(
x = chm,
y = samp,
ID = F,
bind = T
)
# rename
names(samp) <- c("fuel", "chm")
# plot it out
x <- samp$chm
y <- samp$fuel
plot(
x = x,
y = y,
pch = 16,
col = rgb(0,0,0,0.1),
xlab = "canopy height (m)",
ylab = "fuel load (kg)",
las = 1
)
mod <- nls(
y ~ a * (1 - exp(-b * x)),
start = list(
a = 0.1,
b = 0.1
)
)
newx <- seq(0, max(x, na.rm = T), length.out = 1000)
newy <- predict(mod, newdata = list(x = newx))
lines(newx, newy, lwd = 2, col = 2)
Pretty good-looking fit. Now I’ll apply this to the CHM to predict fuel load at fine-scale based on canopy height.
# apply model to chm
names(chm) <- "x"
fuel_fine_regr <- predict(
chm,
mod,
fun = function(model, data, ...) {
predict(model, newdata = as.data.frame(data), ...)
}
)
# plot it out
plot(
fuel_fine_regr,
plg = list(
title = "Fuel Load (kg)"
)
)
Note how the “isolated tree” effect is entirely gone from this new output – values make much more logical sense. Now let’s aggregate this back up to the original resolution for comparative purposes.
# reaggregate for testing
fuel_reagg_regr <- resample(fuel_fine_regr, fuel_coarse, method = "sum") |> crop(e)
# stack them up
fuel_stack <- c(fuel_orig, fuel_reagg_regr)
names(fuel_stack) <- c("original", "reaggregated")
# get global maximum for consistent visualization
fuel_max <- global(fuel_stack, "max") |> max()
# plot them out
plot(fuel_stack, range = c(0, fuel_max))
The re-aggregated results based on the regression model appear noisier than the pure moving window-based approach. Here is the scatterplot, for comparison:
# convert pixels to data.frame
df <- as.data.frame(fuel_stack)
x <- df$original
y <- df$reaggregated
# plot it out
plot(
x = x,
y = y,
pch = 16,
col = rgb(0,0,0,0.1),
xlim = range(df),
ylim = range(df),
xlab = "original",
ylab = "reaggregated",
las = 1
)
grid()
abline(0,1)
mod <- lm(y ~ x)
abline(mod, lwd = 2, col = 2)
r2 <- summary(mod)$r.squared |> formatC(2, format = "f")
rmse <- sqrt(mean((x - y)^2)) |> formatC(2, format = "f")
bias <- mean(y - x) |> formatC(2, format = "f")
r2_txt <- bquote(bold(R^2==.(r2)))
rmse_txt <- bquote(bold(RMSE==.(rmse)))
bias_txt <- bquote(bold(Bias==.(bias)))
legend(
"topleft",
legend = c(r2_txt, rmse_txt, bias_txt),
text.col = 2,
bty = "n",
x.intersp = 0
)
The conclusions are about what we already knew, but now hopefully with a little more visual/quantitative support to inform a discussion:
Possible suggestion for moving forward: The regression approach is basically an “add-on” step to the moving window approach, as you can see if you dig through the code above. And, it really does not add much in the way of processing times at the individual tile level. So, might it be worth adding an argument to the eventual downscaling function in voxelFuel that allows one to run the additional regression step? I would just hate to go all the way down the road of processing massive swaths of the Western US, ultimately yielding some voxel-based fuel load estimates that produce unusable fire behavior predictions. One of the big questions was about how to define the spatial window within which the regression model should be based. If we just applied this at the individual tile level, then perhaps that is a baked-in solution to the problem. It could introduce some slight tile-based edge effects, but perhaps those effects are less harmful in a fire behavior context than individual trees with unrealistically high fuel loads?