The full, hopefully paper-ready analysis
In this, the next (and maybe last?) of the Monroe Mountain Fuel Consumption Modeling R Markdown series, I am to provide a full summary of the entire workflow that hopefully forms the basis of a paper. The objectives of this document, and more broadly, this study, are to:
# import libraries
library(terra)
library(rmarkdown)
library(knitr)
library(sf)
library(lidR)
library(lidRmetrics)
library(VSURF)
library(ranger)
library(dplyr)
library(colorspace)
library(RColorBrewer)
library(stringr)
# set random seed
set.seed(5757)
# define number of cores to be used
ncores <- 30L
# set number of lidar threads
set_lidr_threads(1L)
# suppress progress bars
options(lidR.progress = F)
# define modeling directory
mod.dir <- "S:/ursa/campbell/fasmee/modeling"# read in the lidar data
ctg.pre <- "S:/ursa/campbell/fasmee/data/lidar_data/monroe_mtn_lidar/a04_height" |>
readLAScatalog(filter = "-drop_withheld -drop_class 7 18 -unique_xyz -keep_z -1 50",
progress = F)
ctg.pst.shrt.a <- "S:/ursa/campbell/fasmee/data/lidar_data/monroe_mtn_postfire_lidar/Monroe2020/LAZ_1.4/a03_hgt" |>
readLAScatalog(filter = "-drop_withheld -drop_class 7 18 -unique_xyz -keep_z -1 50",
progress = F)
ctg.pst.shrt.b <- "S:/ursa/campbell/fasmee/data/lidar_data/monroe_mtn_postfire_lidar/Annabella2021/LAZ/a03_hgt" |>
readLAScatalog(filter = "-drop_withheld -drop_class 7 18 -unique_xyz -keep_z -1 50",
progress = F)
ctg.pst.long <- "S:/ursa2/campbell/ems4d/data/als/als_2023/LAZ/Eastburn_processed" |>
readLAScatalog(filter = "-drop_withheld -drop_class 7 18 -unique_xyz -keep_z -1 50",
progress = F)It is worth noting here that I have omitted one outlier plot that features a 300 Mg/ha increase in 1000-hour fuels. This may represent either an error in the data, or if it’s not, certaintly represents both an unlikely scenario and one that can have negative impacts on fuel consumption model performance.
# check to see if fuels data have already been processed
fuel.plots.file <- "S:/ursa/campbell/fasmee/data/from_ryan/plots_wFuels_clean.shp"
if (!file.exists(fuel.plots.file)){
# read in the fuel plot data
fuel.plots <- vect("S:/ursa/campbell/fasmee/data/from_ryan/plots_wFuels.shp")
# remove empty plots
fuel.plots <- fuel.plots[!is.na(fuel.plots$FID),]
# separate pre-fire and post-fire plots
fuel.plots.pre <- fuel.plots[grep("pre", fuel.plots$FID),]
fuel.plots.pst <- fuel.plots[grep("post", fuel.plots$FID),]
# remove unnecessary columns
del.cols <- c("FID", "X", "Y", "cbd", "cbh")
keep.cols <- names(fuel.plots)[!names(fuel.plots) %in% del.cols]
fuel.plots.pre <- fuel.plots.pre[,keep.cols]
fuel.plots.pst <- fuel.plots.pst[,keep.cols]
# move plot ID column from last to first
n <- ncol(fuel.plots.pre)
fuel.plots.pre <- fuel.plots.pre[,c(n,1:(n-1))]
fuel.plots.pst <- fuel.plots.pst[,c(n,1:(n-1))]
# create column name lookup table
cnlut <- matrix(
c("Plot", "plot",
"yr", "year",
"hr1", "hr1",
"hr10", "hr10",
"hr100", "hr100",
"hr1000", "hr1000",
"nonwody", "herb",
"shrubs", "shrub",
"trees", "seed",
"llm", "llm",
"ground", "duff",
"shrb_sd", "shrbsd",
"woody", "woody",
"uf", "nontre",
"cf", "tree",
"acf", "acf",
"tf", "tf",
"af", "af"),
byrow = T, ncol = 2
) |>
as.data.frame()
colnames(cnlut) <- c("old", "new")
# rename columns
for (i in 1:nrow(cnlut)){
old <- cnlut$old[i]
new <- cnlut$new[i]
names(fuel.plots.pre)[names(fuel.plots.pre) == old] <- new
names(fuel.plots.pst)[names(fuel.plots.pst) == old] <- new
}
# add "pre" and "pst" to column names for all attributes
dif.names <- names(fuel.plots.pre)[2:ncol(fuel.plots.pre)]
names(fuel.plots.pre)[2:ncol(fuel.plots.pre)] <-
paste0(names(fuel.plots.pre)[2:ncol(fuel.plots.pre)], "_pre")
names(fuel.plots.pst)[2:ncol(fuel.plots.pst)] <-
paste0(names(fuel.plots.pst)[2:ncol(fuel.plots.pst)], "_pst")
# merge them
fuel.plots <- merge(fuel.plots.pre, as.data.frame(fuel.plots.pst))
# pad single-digit plot ids with 0's and add "p" to keep as character
fuel.plots$plot <- str_pad(fuel.plots$plot, 3, "left", "0")
fuel.plots$plot <- paste0("p", fuel.plots$plot)
# sort of plot id
fuel.plots <- fuel.plots[order(fuel.plots$plot),]
# get differences between pre- and post-fire fuel loads
# (i.e., calc consumption)
for (dif.name in dif.names){
pre.fld <- paste0(dif.name, "_pre")
pst.fld <- paste0(dif.name, "_pst")
dif.fld <- paste0(dif.name, "_dif")
fuel.plots[,dif.fld] <- as.data.frame(fuel.plots[,pre.fld]) -
as.data.frame(fuel.plots[,pst.fld])
}
# remove the plot that had an extreme increase in ground fuels
fuel.plots <- fuel.plots[fuel.plots$hr1000_dif > -300,]
# write to file
writeVector(fuel.plots, fuel.plots.file, overwrite = T)
}
# read the data in
fuel.plots <- vect(fuel.plots.file)The modeling workflow is the same for all models generated in this analysis, and proceeds as follows:
# define modeling function
tuned.model <- function(mod.df){
mod.df <- na.omit(mod.df)
plot.ids <- mod.df$plot
folds <- substr(plot.ids, 2, 2)
mod.df$plot <- NULL
x.cols <- 2:ncol(mod.df)
y.col <- 1
v <- VSURF(
x = mod.df[,x.cols],
y = mod.df[,y.col],
RFimplem = "ranger",
parallel = T,
ncores = ncores,
verbose = F
)
x.cols <- x.cols[v$varselect.pred]
mod.df <- mod.df[,c(y.col, x.cols)]
cv.iters <- 100
hyp.grid <- data.frame(
mtry = sample(1:length(x.cols), cv.iters, T),
min.node.size = sample(1:10, cv.iters, T),
sample.fraction = runif(cv.iters, 0.1, 0.9),
mse = rep(NA, cv.iters)
)
for (i in seq(1,cv.iters)){
mod <- ranger(
dependent.variable.name = colnames(mod.df)[1],
data = mod.df,
mtry = hyp.grid$mtry[i],
sample.fraction = hyp.grid$sample.fraction[i],
min.node.size = hyp.grid$min.node.size[i],
num.threads = ncores
)
hyp.grid$mse[i] <- mod$prediction.error
}
mtry <- hyp.grid$mtry[which.min(hyp.grid$mse)]
min.node.size <- hyp.grid$min.node.size[which.min(hyp.grid$mse)]
sample.fraction <- hyp.grid$sample.fraction[which.min(hyp.grid$mse)]
df.tune <- data.frame(
hyperparameter = c("mtry", "min.node.size", "sample.fraction"),
tuned.value = c(mtry, min.node.size, sample.fraction)
)
for (fold in unique(folds)){
mod.df.valid <- mod.df[folds == fold,]
mod.df.train <- mod.df[folds != fold,]
mod <- ranger(
dependent.variable.name = colnames(mod.df.train)[1],
data = mod.df.train,
mtry = mtry,
sample.fraction = sample.fraction,
min.node.size = min.node.size,
num.threads = ncores
)
pred <- predict(mod, mod.df.valid)$prediction
if (fold == unique(folds)[1]){
df.pred.obs <- data.frame(
plot = plot.ids[folds == fold],
obs = mod.df.valid[,y.col],
prd = pred
)
} else {
df.pred.obs <- rbind(
df.pred.obs,
data.frame(
plot = plot.ids[folds == fold],
obs = mod.df.valid[,y.col],
prd = pred
)
)
}
}
full.mod <- ranger(
dependent.variable.name = colnames(mod.df)[1],
data = mod.df,
mtry = mtry,
sample.fraction = sample.fraction,
min.node.size = min.node.size,
importance = "permutation",
num.threads = ncores
)
df.imp <- importance(full.mod) |>
as.data.frame()
colnames(df.imp) <- "inc.mse"
df.imp$predictor <- row.names(df.imp)
df.imp$perc.inc.mse <- df.imp$inc.mse / full.mod$prediction.error
row.names(df.imp) <- 1:nrow(df.imp)
df.imp <- df.imp[,c("predictor", "inc.mse", "perc.inc.mse")]
return(list(df.pred.obs = df.pred.obs,
df.tune = df.tune,
df.imp = df.imp,
full.mod = full.mod))
}
# define predicted vs. observed assessment and plotting function
prd.obs.plot <- function(df.pred.obs, name, plot = T){
df.pred.obs <- na.omit(df.pred.obs)
x <- df.pred.obs$obs
y <- df.pred.obs$prd
mod <- lm(y ~ x)
a <- coef(mod)[[2]]
b <- coef(mod)[[1]]
r2 <- summary(mod)$adj.r.squared
rmse <- sqrt(mean((x-y)^2))
rrmse <- 100*rmse/mean(x)
bias <- mean(y-x)
rbias <- 100*bias/mean(x)
if (plot == T){
rg <- range(c(x,y), na.rm = T)
par(mar = c(5,5,2,1), las = 1)
plot(x, y, type = "n", xlab = "Observed", ylab = "Predicted", main = name,
ylim = rg, xlim = rg)
grid()
lines(x = c(-1e10, 1e10), y = c(-1e10, 1e10))
points(x, y, pch = 16, col = rgb(0,0,0,0.25), cex = 1.25)
new.x <- range(x)
new.y <- predict(mod, newdata = list(x = new.x))
lines(new.x, new.y, lwd = 2, col = 2)
r2.lab <- bquote(R^2 == .(round(r2,2)))
rrmse.lab <- bquote(symbol("%")*RMSE == .(round(rrmse)))
legend(
"topleft",
legend = c(r2.lab, rrmse.lab),
text.col = 2,
bty = "n",
x.intersp = 0
)
}
return(data.frame(var = name,
a = a,
b = b,
r2 = r2,
rmse = rmse,
rrmse = rrmse,
bias = bias,
rbias = rbias))
}The first major step in the analysis is to determine the best bandwidth for profile density change fuel consumption prediction. I will test every bandwidth between 1 and 25 cm at a 1 cm interval.
# define kernel density function
dens.fun <- function(z, bw){
z[z < 0] <- 0
d <- density(z, bw, from = 0, to = 30, n = length(seq(0,30,0.01)))
dy <- as.list(d$y)
dx <- as.integer(d$x * 100) |>
str_pad(4, "left", "0")
names(dy) <- paste0("d.", dx)
return(dy)
}
# reproject plots to match crs of pre and pst ctgs
fuel.plots.pre <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pre))
fuel.plots.pst.shrt.a <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.shrt.a)) |>
subset(!substr(plot, 2, 2) %in% c("5", "7"))
fuel.plots.pst.shrt.b <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.shrt.b)) |>
subset(substr(plot, 2, 2) %in% c("5", "7"))
fuel.plots.pst.long <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.long))
# begin bandwidth loop
bws <- seq(0.01,0.25,0.01)
i <- 0
n <- length(bws)
for (bw in bws){
# define output csv file
out.csv.pre <- paste0(mod.dir, "/pdc_pre_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv")
out.csv.pst.shrt <- paste0(mod.dir, "/pdc_pst_shrt_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv")
out.csv.dif.shrt <- paste0(mod.dir, "/pdc_dif_shrt_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv")
out.csv.pst.long <- paste0(mod.dir, "/pdc_pst_long_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv")
out.csv.dif.long <- paste0(mod.dir, "/pdc_dif_long_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv")
# check to see if it's already been run
if (!all(file.exists(out.csv.pre),
file.exists(out.csv.pst.shrt),
file.exists(out.csv.dif.shrt),
file.exists(out.csv.pst.long),
file.exists(out.csv.dif.long))){
# run plot metrics
dens.pre <- plot_metrics(
las = ctg.pre,
func = ~dens.fun(Z, bw),
geometry = fuel.plots.pre,
radius = 8
) |>
st_drop_geometry()
dens.pst.shrt.a <- plot_metrics(
las = ctg.pst.shrt.a,
func = ~dens.fun(Z, bw),
geometry = fuel.plots.pst.shrt.a,
radius = 8
) |>
st_drop_geometry()
dens.pst.shrt.b <- plot_metrics(
las = ctg.pst.shrt.b,
func = ~dens.fun(Z, bw),
geometry = fuel.plots.pst.shrt.b,
radius = 8
) |>
st_drop_geometry()
dens.pst.long <- plot_metrics(
las = ctg.pst.long,
func = ~dens.fun(Z, bw),
geometry = fuel.plots.pst.long,
radius = 8
) |>
st_drop_geometry()
# merge the a and b datasets for pst.shrt
dens.pst.shrt <- rbind(dens.pst.shrt.a, dens.pst.shrt.b)
# sort the data
dens.pre <- dens.pre[order(dens.pre$plot),]
dens.pst.shrt <- dens.pst.shrt[order(dens.pst.shrt$plot),]
dens.pst.long <- dens.pst.long[order(dens.pst.long$plot),]
# remove response variables, but keep plot id
dens.pre <- cbind(
dens.pre["plot"],
dens.pre[,!colnames(dens.pre) %in% colnames(fuel.plots.pre)]
)
dens.pst.shrt <- cbind(
dens.pst.shrt["plot"],
dens.pst.shrt[,!colnames(dens.pst.shrt) %in% colnames(fuel.plots.pre)]
)
dens.pst.long <- cbind(
dens.pst.long["plot"],
dens.pst.long[,!colnames(dens.pst.long) %in% colnames(fuel.plots.pre)]
)
# calculate differences
dens.dif.shrt <- cbind(
dens.pre["plot"],
dens.pre[,2:ncol(dens.pre)] - dens.pst.shrt[,2:ncol(dens.pst.shrt)]
)
dens.dif.long <- cbind(
dens.pre["plot"],
dens.pre[,2:ncol(dens.pre)] - dens.pst.long[,2:ncol(dens.pst.long)]
)
# change column names
colnames(dens.pre)[2:ncol(dens.pre)] <-
paste0(colnames(dens.pre)[2:ncol(dens.pre)], ".pre")
colnames(dens.pst.shrt)[2:ncol(dens.pst.shrt)] <-
paste0(colnames(dens.pst.shrt)[2:ncol(dens.pst.shrt)], ".pst.shrt")
colnames(dens.dif.shrt)[2:ncol(dens.dif.shrt)] <-
paste0(colnames(dens.dif.shrt)[2:ncol(dens.dif.shrt)], ".dif.shrt")
colnames(dens.pst.long)[2:ncol(dens.pst.long)] <-
paste0(colnames(dens.pst.long)[2:ncol(dens.pst.long)], ".pst.long")
colnames(dens.dif.long)[2:ncol(dens.dif.long)] <-
paste0(colnames(dens.dif.long)[2:ncol(dens.dif.long)], ".dif.long")
# write to csv
write.csv(dens.pre, out.csv.pre, row.names = F)
write.csv(dens.pst.shrt, out.csv.pst.shrt, row.names = F)
write.csv(dens.dif.shrt, out.csv.dif.shrt, row.names = F)
write.csv(dens.pst.long, out.csv.pst.long, row.names = F)
write.csv(dens.dif.long, out.csv.dif.long, row.names = F)
}
}Below are a few quick visual examples of different bandwidths to highlight the effect that bandwidth has on density profiles.
# define example bws
ex.bws <- c(0.02,0.05,0.10,0.20)
# define example row
ex.row <- 1
# set up plot
par(mfrow = c(2,2), mar = c(5,5,2,1), las = 1)
# loop through bws
for (ex.bw in ex.bws){
# read in the data
dens.pre <- read.csv(
paste0(
mod.dir, "/pdc_pre_bw_",
str_pad(ex.bw * 100, 3, "left", "0"), "cm.csv"
)
)
dens.pst.shrt <- read.csv(
paste0(
mod.dir, "/pdc_pst_shrt_bw_",
str_pad(ex.bw * 100, 3, "left", "0"), "cm.csv"
)
)
# get example row data
dens.pre.ex <- dens.pre[ex.row,2:ncol(dens.pre)]
dens.pst.shrt.ex <- dens.pst.shrt[ex.row,2:ncol(dens.pst.shrt)]
# create empty plot
plot(x = c(0,30), y = c(0,0.1),
type = "n", xlab = "Height (m)", ylab = "Density",
main = paste0("Bandwidth: ", ex.bw * 100, "cm"))
grid()
# add the data
polygon(x = c(0, seq(0,30,0.01), 30),
y = c(0, dens.pre.ex, 0),
col = adjust_transparency(2, 0.5),
border = NA)
lines(x = seq(0,30,0.01), y = dens.pre.ex, col = 2, lwd = 2)
polygon(x = c(0, seq(0,30,0.01), 30),
y = c(0, dens.pst.shrt.ex, 0),
col = adjust_transparency(4, 0.5),
border = NA)
lines(x = seq(0,30,0.01), y = dens.pst.shrt.ex, col = 4, lwd = 2)
# add legend
legend("topright", legend = c("Pre-fire", "Post-fire"), pch = 22,
col = c(2, 4), pt.lwd = 2, pt.cex = 2,
pt.bg = c(adjust_transparency(2, 0.5), adjust_transparency(4, 0.5)))
}
Below I will build models with each of the bandwidths aimed at predicting loads in each of 16 different measured fuel beds. I will do so using two different sets of predictor variables:
The hypothesis is that Set 2 may be superior as it informs the model what pre-fire fuel structure was, instead of relying purely on changes thereto as the basis of prediction.
# define output performance metrics csv and see if it's already been run
out.perf.mets.csv <- paste0(mod.dir, "/fuel_cons_vs_pdc_bw_opt_perf_mets.csv")
if (!file.exists(out.perf.mets.csv)){
# define response variables of interest
vars <- c("hr1", "hr10", "hr100", "hr1000", "herb", "shrub", "seed", "llm",
"duff", "shrbsd", "woody", "nontre", "tree", "acf", "tf", "af")
vars <- paste0(vars, "_dif")
# define bandwidths
bws <- seq(0.01,0.25,0.01)
# loop through response variables
i <- 0
n <- length(vars)
for (var in vars){
# print status
i <- i + 1
print(paste0(date(), " ", var, " (", i, "/", n, ")"))
# loop through bandwidths
j <- 0
m <- length(bws)
for (bw in bws){
# print status
j <- j + 1
print(paste0(date(), " ", bw, " (", j, "/", m, ")"))
# read in the data
dens.pre <- read.csv(
paste0(
mod.dir, "/pdc_pre_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv"
)
)
dens.dif.shrt <- read.csv(
paste0(
mod.dir, "/pdc_dif_shrt_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv"
)
)
dens.dif.long <- read.csv(
paste0(
mod.dir, "/pdc_dif_long_bw_",
str_pad(bw * 100, 3, "left", "0"), "cm.csv"
)
)
#-------------------------shrt
#-----------just dif
# create modeling data frame
x <- dens.dif.shrt
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output csvs and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_pred_obs_dif_shrt.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_tune_dif_shrt.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_imp_dif_shrt.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has not been run yet, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# generate performance metrics from the pred vs. obs data
perf.mets.bw <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.bw$type <- "dif.shrt"
perf.mets.bw$bw <- bw
if (var == vars[1] & bw == bws[1]){
perf.mets.bws <- perf.mets.bw
} else {
perf.mets.bws <- rbind(perf.mets.bws, perf.mets.bw)
}
#-----------pre & dif
# create modeling data frame
x <- merge(dens.pre, dens.dif.shrt)
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output csvs and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_pred_obs_pre_and_dif_shrt.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_tune_pre_and_dif_shrt.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_imp_pre_and_dif_shrt.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has not been run yet, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# generate performance metrics from the pred vs. obs data
perf.mets.bw <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.bw$type <- "pre.and.dif.shrt"
perf.mets.bw$bw <- bw
perf.mets.bws <- rbind(perf.mets.bws, perf.mets.bw)
#-------------------------long
#-----------just dif
# create modeling data frame
x <- dens.dif.long
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output csvs and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_pred_obs_dif_long.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_tune_dif_long.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_imp_dif_long.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has not been run yet, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# generate performance metrics from the pred vs. obs data
perf.mets.bw <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.bw$type <- "dif.long"
perf.mets.bw$bw <- bw
perf.mets.bws <- rbind(perf.mets.bws, perf.mets.bw)
#-----------pre & dif
# create modeling data frame
x <- merge(dens.pre, dens.dif.long)
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output csvs and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_pred_obs_pre_and_dif_long.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_tune_pre_and_dif_long.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_pdc_bw_", str_pad(bw * 100, 3, "left", "0"),
"cm_imp_pre_and_dif_long.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has not been run yet, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# generate performance metrics from the pred vs. obs data
perf.mets.bw <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.bw$type <- "pre.and.dif.long"
perf.mets.bw$bw <- bw
perf.mets.bws <- rbind(perf.mets.bws, perf.mets.bw)
}
}
# write to csv
write.csv(perf.mets.bws, out.perf.mets.csv, row.names = F)
}Below you will see the results of this bandwidth optimization test, comparing both the bandwidth and the difference in model performance between the two sets of predictors described above. Though the optimal bandwidth varies between different fuelbed models, there is a clear general trend that smaller bandwidths outperform larger bandwidths. Looking at the median among all of the models, the trend becomes even more clear. Furthermore, including pre-fire density estimates tends to substantially improve model predictions. The best bandwidth for the models that included the pre-fire density estimates was 1 cm.
# set up plot
layout(matrix(1:34, byrow = T, nrow = 17), widths = c(4,1))
par(oma = c(5,5,3,5), las = 1)
# define response variables of interest
vars <- c("hr1", "hr10", "hr100", "hr1000", "herb", "shrub", "seed", "llm",
"duff", "shrbsd", "woody", "nontre", "tree", "acf", "tf", "af",
"median")
vars <- paste0(vars, "_dif")
# read in the performance data
perf.mets.bws <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_pdc_bw_opt_perf_mets.csv")
)
# filter to just the shrt data
perf.mets.bws <- perf.mets.bws[
perf.mets.bws$type %in% c("dif.shrt", "pre.and.dif.shrt"),
]
# calculate median performance across all variables
perf.mets.bws.med <- perf.mets.bws |>
group_by(type, bw) |>
summarize(r2 = median(r2),
rrmse = median(rrmse),
.groups = "keep") |>
mutate(var = "median_dif") |>
as.data.frame()
# merge
perf.mets.bws.med <- bind_rows(perf.mets.bws, perf.mets.bws.med)
# loop through variables
i <- 0
for (var in vars){
# add one to counter
i <- i + 1
# get rrmse vals
rrmse.dif.shrt <-
perf.mets.bws.med$rrmse[perf.mets.bws.med$var == var &
perf.mets.bws.med$type == "dif.shrt"]
rrmse.pre.dif.shrt <-
perf.mets.bws.med$rrmse[perf.mets.bws.med$var == var &
perf.mets.bws.med$type == "pre.and.dif.shrt"]
# get range
rrmse.rng <- range(c(rrmse.dif.shrt, rrmse.pre.dif.shrt))
# set up plot
par(mar = c(0,0,0,4))
# create empty plot
plot(x = c(1,25), y = rrmse.rng, type = "n", xlab = NA, ylab = NA, xaxt = "n",
yaxt = "n")
# add vertical grid
grid()
# add lines
lines(rrmse.dif.shrt ~ seq(1,25), lwd = 2, col = 2)
lines(rrmse.pre.dif.shrt ~ seq(1,25), lwd = 2, col = 4)
# add vertical lines
abline(v = seq(1,25)[which.min(rrmse.dif.shrt)],
lwd = 7, col = adjust_transparency(2, 0.33))
abline(v = seq(1,25)[which.min(rrmse.pre.dif.shrt)],
lwd = 7, col = adjust_transparency(4, 0.33))
# add axes
if (i %% 2 == 0){
axis(2)
} else {
axis(4)
}
if (var == vars[length(vars)]){
axis(1, at = seq(1,25), labels = F, col.ticks = "lightgray")
axis(1)
}
# add box
if (var == vars[length(vars)]){
box(lwd = 2.5)
} else {
box()
}
# add legends
if (var == vars[1]){
legend(
x = par("usr")[1],
y = par("usr")[4] + 0.05 * diff(par("usr")[3:4]),
legend = c("Density Difference Only", "Lowest Error Bandwidth"),
col = c(2,adjust_transparency(2, 0.33)), lwd = c(2,7), xpd = NA,
xjust = 0, yjust = 0, bty = "n"
)
legend(
x = par("usr")[2],
y = par("usr")[4] + 0.05 * diff(par("usr")[3:4]),
legend = c("Density Difference & Pre-Fire Density",
"Lowest Error Bandwidth"),
col = c(4,adjust_transparency(4, 0.33)), lwd = c(2,7), xpd = NA,
xjust = 1, yjust = 0, bty = "n"
)
}
# add axis labels
if (var == vars[length(vars)]){
mtext("Bandwidth (cm)", 1, 2.5, outer = F)
mtext(expression(symbol("%")*RMSE), 2, 3, outer = T, las = 0)
}
#--------------------boxplot
# set up plot
par(rep(0,2.5))
# get sub df
perf.mets.bws.med.sub <- perf.mets.bws.med[perf.mets.bws.med$var == var,]
# boxplot
boxplot(rrmse ~ type, data = perf.mets.bws.med.sub, yaxt = "n", ylab = NA,
xaxt = "n", yaxt = "n",
col = c(adjust_transparency(2, 0.33), adjust_transparency(4, 0.33)),
border = c(2,4), pch = 16)
# add x axis
if (var == vars[length(vars)]){
axis(1, at = c(1,2), labels = F)
par(xpd = NA)
text(x = c(1,2), y = par("usr")[3] - 0.1 * diff(par("usr")[3:4]),
labels = c("Dens\nDiff\nOnly", "Dens\nDiff\n& Pre\nDens"),
adj = c(0.5,1))
par(xpd = F)
}
# add box
if (var == vars[length(vars)]){
box(lwd = 3)
} else {
box()
}
# add variable label
mtext(
text = sub("_dif", "", var),
side = 4,
line = 0.5,
font = 2
)
}
In 2019, Hu et al. introduced a novel method for mapping fire severity using the difference in the area under the curve of lidar-based height percentile profiles before and after a fire. The Profile Area Change (PAC) method showed promise. To use as a base of comparison against PDC, I will derive the height percentiles needed to calculate PAC below.
# define percentile metrics function
perc.fun <- function(z){
q <- quantile(z, seq(0,1,0.01)) |>
as.list()
names(q) <- paste0("p.", stringr::str_pad(seq(0,100,1), 3, "left", "0"))
return(q)
}
# reproject plots to match crs of pre and pst ctgs
fuel.plots.pre <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pre))
fuel.plots.pst.shrt.a <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.shrt.a)) |>
subset(!substr(plot, 2, 2) %in% c("5", "7"))
fuel.plots.pst.shrt.b <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.shrt.b)) |>
subset(substr(plot, 2, 2) %in% c("5", "7"))
fuel.plots.pst.long <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.long))
# define output csvs
out.csv.pre.shrt <- paste0(mod.dir, "/perc_pre_shrt.csv")
out.csv.pre.long <- paste0(mod.dir, "/perc_pre_long.csv")
out.csv.pst.shrt <- paste0(mod.dir, "/perc_pst_shrt.csv")
out.csv.pst.long <- paste0(mod.dir, "/perc_pst_long.csv")
# check to see if they've already been run
if (!all(file.exists(out.csv.pre.shrt),
file.exists(out.csv.pre.long),
file.exists(out.csv.pst.shrt),
file.exists(out.csv.pst.long))){
# run plot metrics
perc.pre <- plot_metrics(
las = ctg.pre,
func = ~perc.fun(Z),
geometry = fuel.plots.pre,
radius = 8
) |>
st_drop_geometry()
perc.pst.shrt.a <- plot_metrics(
las = ctg.pst.shrt.a,
func = ~perc.fun(Z),
geometry = fuel.plots.pst.shrt.a,
radius = 8
) |>
st_drop_geometry()
perc.pst.shrt.b <- plot_metrics(
las = ctg.pst.shrt.b,
func = ~perc.fun(Z),
geometry = fuel.plots.pst.shrt.b,
radius = 8
) |>
st_drop_geometry()
perc.pst.long <- plot_metrics(
las = ctg.pst.long,
func = ~perc.fun(Z),
geometry = fuel.plots.pst.long,
radius = 8
) |>
st_drop_geometry()
# merge the a and b datasets for pst.shrt
perc.pst.shrt <- rbind(perc.pst.shrt.a, perc.pst.shrt.b)
# sort the data
perc.pre <- perc.pre[order(perc.pre$plot),]
perc.pst.shrt <- perc.pst.shrt[order(perc.pst.shrt$plot),]
perc.pst.long <- perc.pst.long[order(perc.pst.long$plot),]
# remove response variables, but keep plot id
perc.pre <- cbind(
perc.pre["plot"],
perc.pre[,!colnames(perc.pre) %in% colnames(fuel.plots.pre)]
)
perc.pst.shrt <- cbind(
perc.pst.shrt["plot"],
perc.pst.shrt[,!colnames(perc.pst.shrt) %in% colnames(fuel.plots.pre)]
)
perc.pst.long <- cbind(
perc.pst.long["plot"],
perc.pst.long[,!colnames(perc.pst.long) %in% colnames(fuel.plots.pre)]
)
# set negative values to zero
perc.pre[perc.pre < 0] <- 0
perc.pst.shrt[perc.pst.shrt < 0] <- 0
perc.pst.long[perc.pst.long < 0] <- 0
# get plot-level maximum heights
max.hgt.pre <- apply(perc.pre[,2:ncol(perc.pre)], 1, max)
max.hgt.pst.shrt <- apply(perc.pst.shrt[,2:ncol(perc.pst.shrt)], 1, max)
max.hgt.pst.long <- apply(perc.pst.long[,2:ncol(perc.pst.long)], 1, max)
max.hgt.shrt <- apply(cbind(max.hgt.pre, max.hgt.pst.shrt), 1, max)
max.hgt.long <- apply(cbind(max.hgt.pre, max.hgt.pst.long), 1, max)
# make copies of perc.pre for normalizing to either shrt or long
perc.pre.shrt <- perc.pre
perc.pre.long <- perc.pre
# normalize heights to max
for (i in seq(1,nrow(perc.pre.shrt))){
perc.pre.shrt[i,2:ncol(perc.pre.shrt)] <-
perc.pre.shrt[i,2:ncol(perc.pre.shrt)] / max.hgt.shrt[i]
}
for (i in seq(1,nrow(perc.pre.long))){
perc.pre.long[i,2:ncol(perc.pre.long)] <-
perc.pre.long[i,2:ncol(perc.pre.long)] / max.hgt.long[i]
}
for (i in seq(1,nrow(perc.pst.shrt))){
perc.pst.shrt[i,2:ncol(perc.pst.shrt)] <-
perc.pst.shrt[i,2:ncol(perc.pst.shrt)] / max.hgt.shrt[i]
}
for (i in seq(1,nrow(perc.pst.long))){
perc.pst.long[i,2:ncol(perc.pst.long)] <-
perc.pst.long[i,2:ncol(perc.pst.long)] / max.hgt.long[i]
}
# write to files
write.csv(perc.pre.shrt, out.csv.pre.shrt, row.names = F)
write.csv(perc.pre.long, out.csv.pre.long, row.names = F)
write.csv(perc.pst.shrt, out.csv.pst.shrt, row.names = F)
write.csv(perc.pst.long, out.csv.pst.long, row.names = F)
}Below you will see an example of PAC for the first plot in the fuels database, demonstrating the concept of how PAC captures fuel consumption.
# read the data in
perc.pre.shrt <- read.csv(paste0(mod.dir, "/perc_pre_shrt.csv"))
perc.pst.shrt <- read.csv(paste0(mod.dir, "/perc_pst_shrt.csv"))
# plot an example out
ex.row <- 1
perc.pre.shrt.ex <- perc.pre.shrt[ex.row,2:ncol(perc.pre.shrt)]
perc.pst.shrt.ex <- perc.pst.shrt[ex.row,2:ncol(perc.pst.shrt)]
par(mar = c(5,5,1,1), las = 1)
plot(x = c(0,100), y = range(c(perc.pre.shrt.ex, perc.pst.shrt.ex)),
type = "n", xlab = "Percentile", ylab = "Normalized Height")
grid()
polygon(x = c(0, seq(0,100), 100),
y = c(0, perc.pre.shrt.ex, 0),
col = adjust_transparency(2, 0.5),
border = NA)
lines(x = seq(0,100), y = perc.pre.shrt.ex, col = 2, lwd = 2)
polygon(x = c(0, seq(0,100), 100),
y = c(0, perc.pst.shrt.ex, 0),
col = adjust_transparency(4, 0.5),
border = NA)
lines(x = seq(0,100), y = perc.pst.shrt.ex, col = 4, lwd = 2)
legend("topleft", legend = c("Pre-fire", "Post-fire"), col = c(2, 4),
fill = c(adjust_transparency(2, 0.5), adjust_transparency(4, 0.5)))
While PAC is a nice concept, as a singular measure, it cannot account for differences in consumption within different fuelbeds. To address this limitation, I am introduced Sliced Profile Area Change (SPAC), which begins exactly the same way PAC does, except instead of taking the entire area under the curve (AUC), the AUC is calculated only for a subset of the full vertical profile. The subsets are defined by height percentiles. Below, I will calculate AUC for every possible range of percentiles between the 0th and the 100th, at a 1 percentile interval (e.g., 0th - 1st, 0th - 2nd, 1st - 2nd, etc.). These will be used as candidate predictors for modeling consumption in the different measured fuelbed loads. The figure below illustrates the SPAC concept, including a depiction of all of the slices used in this analysis.
# read the data in
perc.pre.shrt <- read.csv(paste0(mod.dir, "/perc_pre_shrt.csv"))
perc.pst.shrt <- read.csv(paste0(mod.dir, "/perc_pst_shrt.csv"))
# plot an example out
ex.row <- 1
perc.pre.shrt.ex <- perc.pre.shrt[ex.row,2:ncol(perc.pre.shrt)]
perc.pst.shrt.ex <- perc.pst.shrt[ex.row,2:ncol(perc.pst.shrt)]
layout(matrix(c(1,2), nrow = 2), heights = c(1,3))
par(mar = c(0,5,1,1), las = 1)
n.slices <- 0
for (i in 0:100){
for (j in 0:100){
if (j > i){
n.slices <- n.slices + 1
}
}
}
plot(x = c(0,100), y = c(0,n.slices), type = "n", xlab = NA, ylab = NA,
xaxt = "n", yaxt = "n", bty = "n", yaxs = "i")
y <- 0
for (i in 0:100){
for (j in 0:100){
if (j > i){
y <- y + 1
lines(x = c(i,j), y = c(y,y), col = "lightgray")
}
}
}
ex.i <- 35
ex.j <- 65
y <- 0
for (i in 0:100){
for (j in 0:100){
if (j > i){
y <- y + 1
if (i == ex.i & j == ex.j){
lines(x = c(i,j), y = c(y,y), lwd = 2)
points(x = i, y = y, pch = 16)
points(x = j, y = y, pch = 16)
rect(xleft = i, ybottom = 0, xright = j, ytop = y,
col = rgb(0,0,0,0.25), border = NA)
lines(x = c(i,i), y = c(0,y), lty = 2)
lines(x = c(j,j), y = c(0,y), lty = 2)
text(x = mean(c(i,j)), y = y, pos = 1, font = 2, col = "white",
labels = paste0("(", i, "th - ", j, "th)"))
}
}
}
}
legend("topleft", legend = "Profile Slice", lty = 1, col = "lightgray",
bty = "n")
par(mar = c(5,5,0,1))
plot(x = c(0,100), y = range(c(perc.pre.shrt.ex, perc.pst.shrt.ex)),
type = "n", xlab = "Percentile", ylab = "Normalized Height")
grid()
polygon(x = c(0, seq(0,100), 100),
y = c(0, perc.pre.shrt.ex, 0),
col = adjust_transparency(2, 0.5),
border = NA)
lines(x = seq(0,100), y = perc.pre.shrt.ex, col = 2, lwd = 2)
polygon(x = c(0, seq(0,100), 100),
y = c(0, perc.pst.shrt.ex, 0),
col = adjust_transparency(4, 0.5),
border = NA)
lines(x = seq(0,100), y = perc.pst.shrt.ex, col = 4, lwd = 2)
rect(xleft = ex.i, ybottom = 0, xright = ex.j, ytop = 1e6,
col = rgb(0,0,0,0.25), border = NA)
lines(x = c(ex.i,ex.i), y = c(0,1e6), lty = 2)
lines(x = c(ex.j,ex.j), y = c(0,1e6), lty = 2)
legend("topleft", legend = c("Pre-fire", "Post-fire"), col = c(2, 4),
fill = c(adjust_transparency(2, 0.5), adjust_transparency(4, 0.5)))
The code below calculates SPAC for all of the target height percentile slices.
# create auc function
auc.slice.fun <- function(y, slice.min, slice.max){
x <- seq(0,100)
auc <- integrate(splinefun(x, y), slice.min, slice.max)$value
return(auc)
}
#---------------------shrt
# define output csv
out.csv.spac.shrt <- paste0(mod.dir, "/spac_shrt.csv")
# see if it's already been run
if (!file.exists(out.csv.spac.shrt)){
# begin slice loop (i = floor; j = ceiling)
successes <- 0
for (i in 0:100){
for (j in 0:100){
if (j > i){
test <- tryCatch({
auc.slice.pre.shrt <- apply(
perc.pre.shrt[,2:ncol(perc.pre.shrt)], 1,
function(x) auc.slice.fun(x, i, j)
)
auc.slice.pst.shrt <- apply(
perc.pst.shrt[,2:ncol(perc.pst.shrt)], 1,
function(x) auc.slice.fun(x, i, j)
)
spac.shrt <- auc.slice.pre.shrt - auc.slice.pst.shrt
}, error = function(e){
message(i, "-", j, ": ", e)
return(e)
})
if (inherits(test, "error")){
next
} else {
successes <- successes + 1
}
if (successes == 1){
spac.shrt.df <- data.frame(plot = perc.pre.shrt$plot, x = spac.shrt)
colnames(spac.shrt.df)[2] <- paste0(
"spac.shrt.", str_pad(i, 3, "left", "0"),
".", str_pad(j, 3, "left", "0")
)
} else {
spac.shrt.df <- cbind(spac.shrt.df, data.frame(x = spac.shrt))
colnames(spac.shrt.df)[ncol(spac.shrt.df)] <- paste0(
"spac.shrt.", str_pad(i, 3, "left", "0"),
".", str_pad(j, 3, "left", "0")
)
}
}
}
}
# write to file
write.csv(spac.shrt.df, out.csv.spac.shrt, row.names = F)
}
#---------------------long
# define output csv
out.csv.spac.long <- paste0(mod.dir, "/spac_long.csv")
# see if it's already been run
if (!file.exists(out.csv.spac.long)){
# begin slice loop (i = floor; j = ceiling)
successes <- 0
for (i in 0:100){
for (j in 0:100){
if (j > i){
test <- tryCatch({
auc.slice.pre.long <- apply(
perc.pre.long[,2:ncol(perc.pre.long)], 1,
function(x) auc.slice.fun(x, i, j)
)
auc.slice.pst.long <- apply(
perc.pst.long[,2:ncol(perc.pst.long)], 1,
function(x) auc.slice.fun(x, i, j)
)
spac.long <- auc.slice.pre.long - auc.slice.pst.long
}, error = function(e){
message(i, "-", j, ": ", e)
return(e)
})
if (inherits(test, "error")){
next
} else {
successes <- successes + 1
}
if (successes == 1){
spac.long.df <- data.frame(plot = perc.pre.long$plot, x = spac.long)
colnames(spac.long.df)[2] <- paste0(
"spac.long.", str_pad(i, 3, "left", "0"),
".", str_pad(j, 3, "left", "0")
)
} else {
spac.long.df <- cbind(spac.long.df, data.frame(x = spac.long))
colnames(spac.long.df)[ncol(spac.long.df)] <- paste0(
"spac.long.", str_pad(i, 3, "left", "0"),
".", str_pad(j, 3, "left", "0")
)
}
}
}
}
# write to file
write.csv(spac.long.df, out.csv.spac.long, row.names = F)
}In the code below, I take all of the SPAC slices and use them as candidate predictors for fuel consumption for all of the fuelbeds in the field data. As with PDC, I test two sets of predictors:
The hypothesis is that Set 2 may be superior as it informs the model what pre-fire fuel structure was, instead of relying purely on changes thereto as the basis of prediction.
# read in the data
spac.shrt.df <- read.csv(paste0(mod.dir, "/spac_shrt.csv"))
spac.long.df <- read.csv(paste0(mod.dir, "/spac_long.csv"))
perc.pre.shrt <- read.csv(paste0(mod.dir, "/perc_pre_shrt.csv"))
perc.pre.long <- read.csv(paste0(mod.dir, "/perc_pre_long.csv"))
# define output performance metrics csv and see if it's already been run
out.perf.mets.csv <- paste0(mod.dir, "/fuel_cons_vs_spac_perf_mets.csv")
if (!file.exists(out.perf.mets.csv)){
# define response variables of interest
vars <- c("hr1", "hr10", "hr100", "hr1000", "herb", "shrub", "seed", "llm",
"duff", "shrbsd", "woody", "nontre", "tree", "acf", "tf", "af")
vars <- paste0(vars, "_dif")
# loop through response variables
i <- 0
n <- length(vars)
for (var in vars){
# print status
i <- i + 1
print(paste0(date(), " ", var, " (", i, "/", n, ")"))
#-------------------------shrt
#-----------just spac
# create modeling data frame
x <- spac.shrt.df
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output files and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_spac_pred_obs_spac_shrt.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_spac_tune_spac_shrt.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_spac_imp_spac_shrt.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get the outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# save the outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has already been run, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from the pred vs. obs data
perf.mets.temp <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.temp$type <- "spac.shrt"
if (var == vars[1]){
perf.mets <- perf.mets.temp
} else {
perf.mets <- rbind(perf.mets, perf.mets.temp)
}
#-----------pre & dif
# create modeling data frame
x <- merge(perc.pre.shrt, spac.shrt.df)
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output files and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_spac_pred_obs_pre_and_spac_shrt.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_spac_tune_pre_and_spac_shrt.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_spac_imp_pre_and_spac_shrt.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# build model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has already been run, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from the pred vs. obs data
perf.mets.temp <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.temp$type <- "pre.and.spac.shrt"
perf.mets <- rbind(perf.mets, perf.mets.temp)
#-------------------------long
#-----------just spac
# create modeling data frame
x <- spac.long.df
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output files and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_spac_pred_obs_spac_long.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_spac_tune_spac_long.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_spac_imp_spac_long.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get the outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# save the outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has already been run, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from the pred vs. obs data
perf.mets.temp <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.temp$type <- "spac.long"
perf.mets <- rbind(perf.mets, perf.mets.temp)
#-----------pre & dif
# create modeling data frame
x <- merge(perc.pre.long, spac.long.df)
y <- fuel.plots[,c("plot", var)] |> as.data.frame()
mod.df <- merge(y, x)
# define output files and check if it has already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_spac_pred_obs_pre_and_spac_long.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_spac_tune_pre_and_spac_long.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_spac_imp_pre_and_spac_long.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# build model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has already been run, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from the pred vs. obs data
perf.mets.temp <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.temp$type <- "pre.and.spac.long"
perf.mets <- rbind(perf.mets, perf.mets.temp)
}
# write to csv
write.csv(perf.mets, out.perf.mets.csv, row.names = F)
}Below you can see a plot comparing the effect of including (or excluding) pre-fire height percentiles along with the SPAC values in the set of candidate predictors. Including pre-fire height percentiles appears to improve model performance slightly, but not by a large margin.
# read in the data
perf.mets <- read.csv(paste0(mod.dir, "/fuel_cons_vs_spac_perf_mets.csv"))
# subset to just shrt
perf.mets <- perf.mets[grep("shrt", perf.mets$type),]
# for plotting purposes, remove outliers
perf.mets <- perf.mets[perf.mets$rrmse < 500,]
# convert type to factor
perf.mets$type <- factor(
perf.mets$typ, levels = c("spac.shrt", "pre.and.spac.shrt"),
labels = c("SPAC Only", "SPAC & Pre-Perc")
)
# plot it out
boxplot(rrmse ~ type, data = perf.mets,
col = c(adjust_transparency(2, 0.33), adjust_transparency(4, 0.33)),
border = c(2,4), pch = 16, ylab = "%RMSE", xlab = NA)
Whereas PDC and SPAC are focused on using change in lidar structure as the primary basis of fuel consumption modeling, previous efforts have focused more on comparing pre-fire and post-fire modeled fuel loads. The assumption here is that if you can generate an accurate map of pre-fire fuel loading and an accurate map of post-fire fuel loading, then a simple subtraction between the two should yield an accurate fuel consumption estimate. Indeed, it has been shown to be successful by McCarley et al. (2024), including using these same field and lidar data. There are two different approaches for doing so:
Both of these models will be built with lidar structural metrics from Tompalski and Goodbody’s lidRmetrics R library. The code below derives these lidar structural metrics within point clouds clipped to 8 m-buffered plot areas, to match the approximate dimensions of the fuel plots.
# reproject plots to match crs of pre and pst ctgs
fuel.plots.pre <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pre))
fuel.plots.pst.shrt.a <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.shrt.a)) |>
subset(!substr(plot, 2, 2) %in% c("5", "7"))
fuel.plots.pst.shrt.b <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.shrt.b)) |>
subset(substr(plot, 2, 2) %in% c("5", "7"))
fuel.plots.pst.long <- fuel.plots |>
st_as_sf() |>
st_transform(st_crs(ctg.pst.long))
# define output csvs
out.csv.pre <- paste0(mod.dir, "/lidRmetrics_pre.csv")
out.csv.pst.shrt <- paste0(mod.dir, "/lidRmetrics_pst_shrt.csv")
out.csv.pst.long <- paste0(mod.dir, "/lidRmetrics_pst_long.csv")
# check to see if they've already been run
if (!all(file.exists(out.csv.pre),
file.exists(out.csv.pst.shrt),
file.exists(out.csv.pst.long))){
# run plot metrics
mets.pre <- plot_metrics(
las = ctg.pre,
func = ~lidRmetrics::metrics_set3(
x = X,
y = Y,
z = Z,
i = Intensity,
ReturnNumber = ReturnNumber,
NumberOfReturns = NumberOfReturns
),
geometry = fuel.plots.pre,
radius = 8
) |>
st_drop_geometry()
mets.pst.shrt.a <- plot_metrics(
las = ctg.pst.shrt.a,
func = ~lidRmetrics::metrics_set3(
x = X,
y = Y,
z = Z,
i = Intensity,
ReturnNumber = ReturnNumber,
NumberOfReturns = NumberOfReturns
),
geometry = fuel.plots.pst.shrt.a,
radius = 8
) |>
st_drop_geometry()
mets.pst.shrt.b <- plot_metrics(
las = ctg.pst.shrt.b,
func = ~lidRmetrics::metrics_set3(
x = X,
y = Y,
z = Z,
i = Intensity,
ReturnNumber = ReturnNumber,
NumberOfReturns = NumberOfReturns
),
geometry = fuel.plots.pst.shrt.b,
radius = 8
) |>
st_drop_geometry()
mets.pst.long <- plot_metrics(
las = ctg.pst.long,
func = ~lidRmetrics::metrics_set3(
x = X,
y = Y,
z = Z,
i = Intensity,
ReturnNumber = ReturnNumber,
NumberOfReturns = NumberOfReturns
),
geometry = fuel.plots.pst.long,
radius = 8
) |>
st_drop_geometry()
# merge the a and b datasets for the shrt data
mets.pst.shrt <- rbind(mets.pst.shrt.a, mets.pst.shrt.b)
# sort the data
mets.pre <- mets.pre[order(mets.pre$plot),]
mets.pst.shrt <- mets.pst.shrt[order(mets.pst.shrt$plot),]
mets.pst.long <- mets.pst.long[order(mets.pst.long$plot),]
# remove response variables, but keep plot id
mets.pre <- cbind(
mets.pre["plot"],
mets.pre[,!colnames(mets.pre) %in% colnames(fuel.plots.pre)]
)
mets.pst.shrt <- cbind(
mets.pst.shrt["plot"],
mets.pst.shrt[,!colnames(mets.pst.shrt) %in% colnames(fuel.plots.pre)]
)
mets.pst.long <- cbind(
mets.pst.long["plot"],
mets.pst.long[,!colnames(mets.pst.long) %in% colnames(fuel.plots.pre)]
)
# remove unscalable variables
del.cols <- c("n", "n_first", "n_intermediate", "n_last", "n_single",
"n_multiple", "n_return_1", "n_return_2", "n_return_3",
"n_return_4", "vn")
# remove predictors with no variability
sds.pre <- apply(mets.pre[,2:ncol(mets.pre)], 2, sd)
sds.pst.shrt <- apply(mets.pst.shrt[,2:ncol(mets.pst.shrt)], 2, sd)
sds.pst.long <- apply(mets.pst.long[,2:ncol(mets.pst.long)], 2, sd)
del.cols.pre <- c(
del.cols,
names(sds.pre)[sds.pre == 0 & !is.na(sds.pre)]
)
del.cols.pst.shrt <- c(
del.cols,
names(sds.pst.shrt)[sds.pst.shrt == 0 & !is.na(sds.pst.shrt)]
)
del.cols.pst.long <- c(
del.cols,
names(sds.pst.long)[sds.pst.long == 0 & !is.na(sds.pst.long)]
)
# remove predictors with na values
nas.pre <- apply(
mets.pre[,2:ncol(mets.pre)], 2, function(x) sum(is.na(x))
)
nas.pst.shrt <- apply(
mets.pst.shrt[,2:ncol(mets.pst.shrt)], 2, function(x) sum(is.na(x))
)
nas.pst.long <- apply(
mets.pst.long[,2:ncol(mets.pst.long)], 2, function(x) sum(is.na(x))
)
del.cols.pre <- c(
del.cols.pre, names(nas.pre)[nas.pre > 0]
)
del.cols.pst.shrt <- c(
del.cols.pst.shrt, names(nas.pst.shrt)[nas.pst.shrt > 0]
)
del.cols.pst.long <- c(
del.cols.pst.long, names(nas.pst.long)[nas.pst.long > 0]
)
# keep columns
keep.cols.pre <-
colnames(mets.pre)[!colnames(mets.pre) %in% del.cols.pre]
keep.cols.pst.shrt <-
colnames(mets.pst.shrt)[!colnames(mets.pst.shrt) %in% del.cols.pst.shrt]
keep.cols.pst.long <-
colnames(mets.pst.long)[!colnames(mets.pst.long) %in% del.cols.pst.long]
mets.pre <- mets.pre[,keep.cols.pre]
mets.pst.shrt <- mets.pst.shrt[,keep.cols.pst.shrt]
mets.pst.long <- mets.pst.long[,keep.cols.pst.long]
# write to files
write.csv(mets.pre, out.csv.pre, row.names = F)
write.csv(mets.pst.shrt, out.csv.pst.shrt, row.names = F)
write.csv(mets.pst.long, out.csv.pst.long, row.names = F)
}The code below runs the pre-post models.
# read in the data
mets.pre <- read.csv(paste0(mod.dir, "/lidRmetrics_pre.csv"))
mets.pst.shrt <- read.csv(paste0(mod.dir, "/lidRmetrics_pst_shrt.csv"))
mets.pst.long <- read.csv(paste0(mod.dir, "/lidRmetrics_pst_long.csv"))
# define output performance metrics csv and see if they've already been run
out.perf.mets.pre.csv <- paste0(
mod.dir, "/fuel_cons_vs_prepost_shrt_perf_mets_pre.csv"
)
out.perf.mets.pst.shrt.csv <- paste0(
mod.dir, "/fuel_cons_vs_prepost_shrt_perf_mets_pst_shrt.csv"
)
out.perf.mets.dif.shrt.csv <- paste0(
mod.dir, "/fuel_cons_vs_prepost_shrt_perf_mets_dif_shrt.csv"
)
out.perf.mets.pst.long.csv <- paste0(
mod.dir, "/fuel_cons_vs_prepost_long_perf_mets_pst_long.csv"
)
out.perf.mets.dif.long.csv <- paste0(
mod.dir, "/fuel_cons_vs_prepost_long_perf_mets_dif_long.csv"
)
if (!all(file.exists(out.perf.mets.pre.csv),
file.exists(out.perf.mets.pst.shrt.csv),
file.exists(out.perf.mets.dif.shrt.csv),
file.exists(out.perf.mets.pst.long.csv),
file.exists(out.perf.mets.dif.long.csv))){
# define response variables of interest
vars <- c("hr1", "hr10", "hr100", "hr1000", "herb", "shrub", "seed", "llm",
"duff", "shrbsd", "woody", "nontre", "tree", "acf", "tf", "af")
# loop through response variables
i <- 0
n <- length(vars)
for (var in vars){
# print status
i <- i + 1
print(paste0(date(), " ", var, " (", i, "/", n, ")"))
#-----------pre
# define variable of interest
var.pre <- paste0(var, "_pre")
# create modeling data frame
x <- mets.pre
y <- fuel.plots[,c("plot", var.pre)] |> as.data.frame()
mod.df <- merge(y, x)
# define output csv files and check if they've already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var.pre, "_vs_prepost_pred_obs.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var.pre, "_vs_prepost_tune.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var.pre, "_vs_prepost_imp.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has already been run, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from the pred vs. obs data
perf.mets.pre.temp <- prd.obs.plot(df.pred.obs, var.pre, plot = F)
# compile results
perf.mets.pre.temp$type <- "prepost.pre"
if (var == vars[1]){
perf.mets.pre <- perf.mets.pre.temp
} else {
perf.mets.pre <- rbind(perf.mets.pre, perf.mets.pre.temp)
}
#-----------pst.shrt
# define variable of interest
var.pst <- paste0(var, "_pst")
# create modeling data frame
x <- mets.pst.shrt
y <- fuel.plots[,c("plot", var.pst)] |> as.data.frame()
mod.df <- merge(y, x)
# define output files and see if they've already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var.pst, "_vs_prepost_pred_obs_shrt.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var.pst, "_vs_prepost_tune_shrt.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var.pst, "_vs_prepost_imp_shrt.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has already been run, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from the pred vs. obs data
perf.mets.pst.shrt.temp <- prd.obs.plot(df.pred.obs, var.pst, plot = F)
# compile results
perf.mets.pst.shrt.temp$type <- "prepost.pst.shrt"
if (var == vars[1]){
perf.mets.pst.shrt <- perf.mets.pst.shrt.temp
} else {
perf.mets.pst.shrt <- rbind(perf.mets.pst.shrt, perf.mets.pst.shrt.temp)
}
#-----------pst.long
# define variable of interest
var.pst <- paste0(var, "_pst")
# create modeling data frame
x <- mets.pst.long
y <- fuel.plots[,c("plot", var.pst)] |> as.data.frame()
mod.df <- merge(y, x)
# define output files and see if they've already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var.pst, "_vs_prepost_pred_obs_long.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var.pst, "_vs_prepost_tune_long.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var.pst, "_vs_prepost_imp_long.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if it has already been run, read in the pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from the pred vs. obs data
perf.mets.pst.long.temp <- prd.obs.plot(df.pred.obs, var.pst, plot = F)
# compile results
perf.mets.pst.long.temp$type <- "prepost.pst.long"
if (var == vars[1]){
perf.mets.pst.long <- perf.mets.pst.long.temp
} else {
perf.mets.pst.long <- rbind(perf.mets.pst.long, perf.mets.pst.long.temp)
}
#-----------dif.shrt
# define variable of interest
var.dif <- paste0(var, "_dif")
# read in the data
df.pred.obs.pre <- read.csv(
paste0(mod.dir, "/", var.pre, "_vs_prepost_pred_obs.csv")
)
df.pred.obs.pst.shrt <- read.csv(
paste0(mod.dir, "/", var.pst, "_vs_prepost_pred_obs_shrt.csv")
)
# get observed differences
df.obs <- fuel.plots[,c("plot", var.dif)] |> as.data.frame()
colnames(df.obs)[2] <- "obs"
# get predicted differences
df.pred.pre <- df.pred.obs.pre[,c("plot", "prd")]
colnames(df.pred.pre)[2] <- "prd.pre"
df.pred.pst.shrt <- df.pred.obs.pst.shrt[,c("plot", "prd")]
colnames(df.pred.pst.shrt)[2] <- "prd.pst.shrt"
df.pred <- merge(df.pred.pre, df.pred.pst.shrt, all = T)
df.pred$prd <- df.pred$prd.pre - df.pred$prd.pst.shrt
# merge the two
df.pred.obs <- merge(df.obs, df.pred)
df.pred.obs <- df.pred.obs[,c("plot", "obs", "prd")]
# write to file
write.csv(
df.pred.obs,
paste0(mod.dir, "/", var.dif, "_vs_prepost_pred_obs_shrt.csv"),
row.names = F
)
# get performance metrics
perf.mets.dif.shrt.temp <- prd.obs.plot(df.pred.obs, var.dif, plot = F)
# compile results
perf.mets.dif.shrt.temp$type <- "prepost.dif.shrt"
if (var == vars[1]){
perf.mets.dif.shrt <- perf.mets.dif.shrt.temp
} else {
perf.mets.dif.shrt <- rbind(perf.mets.dif.shrt, perf.mets.dif.shrt.temp)
}
#-----------dif.long
# define variable of interest
var.dif <- paste0(var, "_dif")
# read in the data
df.pred.obs.pre <- read.csv(
paste0(mod.dir, "/", var.pre, "_vs_prepost_pred_obs.csv")
)
df.pred.obs.pst.long <- read.csv(
paste0(mod.dir, "/", var.pst, "_vs_prepost_pred_obs_long.csv")
)
# get observed differences
df.obs <- fuel.plots[,c("plot", var.dif)] |> as.data.frame()
colnames(df.obs)[2] <- "obs"
# get predicted differences
df.pred.pre <- df.pred.obs.pre[,c("plot", "prd")]
colnames(df.pred.pre)[2] <- "prd.pre"
df.pred.pst.long <- df.pred.obs.pst.long[,c("plot", "prd")]
colnames(df.pred.pst.long)[2] <- "prd.pst.long"
df.pred <- merge(df.pred.pre, df.pred.pst.long, all = T)
df.pred$prd <- df.pred$prd.pre - df.pred$prd.pst.long
# merge the two
df.pred.obs <- merge(df.obs, df.pred)
df.pred.obs <- df.pred.obs[,c("plot", "obs", "prd")]
# write to file
write.csv(
df.pred.obs,
paste0(mod.dir, "/", var.dif, "_vs_prepost_pred_obs_long.csv"),
row.names = F
)
# get performance metrics
perf.mets.dif.long.temp <- prd.obs.plot(df.pred.obs, var.dif, plot = F)
# compile results
perf.mets.dif.long.temp$type <- "prepost.dif.long"
if (var == vars[1]){
perf.mets.dif.long <- perf.mets.dif.long.temp
} else {
perf.mets.dif.long <- rbind(perf.mets.dif.long, perf.mets.dif.long.temp)
}
}
# write to csvs
write.csv(perf.mets.pre, out.perf.mets.pre.csv, row.names = F)
write.csv(perf.mets.pst.shrt, out.perf.mets.pst.shrt.csv, row.names = F)
write.csv(perf.mets.pst.long, out.perf.mets.pst.long.csv, row.names = F)
write.csv(perf.mets.dif.shrt, out.perf.mets.dif.shrt.csv, row.names = F)
write.csv(perf.mets.dif.long, out.perf.mets.dif.long.csv, row.names = F)
}The code below runs the pooled models.
# read in the data
mets.pre <- read.csv(paste0(mod.dir, "/lidRmetrics_pre.csv"))
mets.pst.shrt <- read.csv(paste0(mod.dir, "/lidRmetrics_pst_shrt.csv"))
mets.pst.long <- read.csv(paste0(mod.dir, "/lidRmetrics_pst_long.csv"))
# add "pre" and "pst" to plot ids in lidRmetrics
mets.pre$plot <- paste0(mets.pre$plot, "_pre")
mets.pst.shrt$plot <- paste0(mets.pst.shrt$plot, "_pst")
mets.pst.long$plot <- paste0(mets.pst.long$plot, "_pst")
# select only the common columns between pre and post
comm.cols <- colnames(mets.pre)[colnames(mets.pre) %in%
colnames(mets.pst.shrt)]
comm.cols <- comm.cols[comm.cols %in% colnames(mets.pst.long)]
mets.pre <- mets.pre[,comm.cols]
mets.pst.shrt <- mets.pst.shrt[,comm.cols]
mets.pst.long <- mets.pst.long[,comm.cols]
# define output performance metrics csv and see if they've already been run
out.perf.mets.shrt.csv <- paste0(
mod.dir, "/fuel_cons_vs_pooled_perf_mets_shrt.csv"
)
out.perf.mets.long.csv <- paste0(
mod.dir, "/fuel_cons_vs_pooled_perf_mets_long.csv"
)
out.perf.mets.dif.shrt.csv <- paste0(
mod.dir, "/fuel_cons_vs_pooled_perf_mets_dif_shrt.csv"
)
out.perf.mets.dif.long.csv <- paste0(
mod.dir, "/fuel_cons_vs_pooled_perf_mets_dif_long.csv"
)
if (!all(file.exists(out.perf.mets.shrt.csv),
file.exists(out.perf.mets.long.csv),
file.exists(out.perf.mets.dif.shrt.csv),
file.exists(out.perf.mets.dif.long.csv))){
# define response variables of interest
vars <- c("hr1", "hr10", "hr100", "hr1000", "herb", "shrub", "seed", "llm",
"duff", "shrbsd", "woody", "nontre", "tree", "acf", "tf", "af")
# loop through response variables
i <- 0
n <- length(vars)
for (var in vars){
# print status
i <- i + 1
print(paste0(date(), " ", var, " (", i, "/", n, ")"))
#-----------shrt
# define variables of interest
var.pre <- paste0(var, "_pre")
var.pst <- paste0(var, "_pst")
# remove "pre" and "pst" from fuels data colmns
y.pre <- fuel.plots[,c("plot", var.pre)] |> as.data.frame()
y.pst <- fuel.plots[,c("plot", var.pst)] |> as.data.frame()
colnames(y.pre)[2] <- var
colnames(y.pst)[2] <- var
# add "pre" and "pst" to plot ids in fuels data
y.pre$plot <- paste0(y.pre$plot, "_pre")
y.pst$plot <- paste0(y.pst$plot, "_pst")
# create modeling data frame
x <- bind_rows(mets.pre, mets.pst.shrt)
y <- bind_rows(y.pre, y.pst)
mod.df <- merge(y, x)
# define output files and check if they've already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_pooled_pred_obs_shrt.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_pooled_tune_shrt.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_pooled_imp_shrt.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if they've already been run, read in pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from pred vs. obs data
perf.mets.shrt.temp <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.shrt.temp$type <- "pooled.shrt"
if (var == vars[1]){
perf.mets.shrt <- perf.mets.shrt.temp
} else {
perf.mets.shrt <- rbind(perf.mets.shrt, perf.mets.shrt.temp)
}
#-----------long
# remove "pre" and "pst" from fuels data colmns
y.pre <- fuel.plots[,c("plot", var.pre)] |> as.data.frame()
y.pst <- fuel.plots[,c("plot", var.pst)] |> as.data.frame()
colnames(y.pre)[2] <- var
colnames(y.pst)[2] <- var
# add "pre" and "pst" to plot ids in fuels data
y.pre$plot <- paste0(y.pre$plot, "_pre")
y.pst$plot <- paste0(y.pst$plot, "_pst")
# create modeling data frame
x <- bind_rows(mets.pre, mets.pst.long)
y <- bind_rows(y.pre, y.pst)
mod.df <- merge(y, x)
# define output files and check if they've already been run
out.csv.pred.obs <- paste0(
mod.dir, "/", var, "_vs_pooled_pred_obs_long.csv"
)
out.csv.tune <- paste0(
mod.dir, "/", var, "_vs_pooled_tune_long.csv"
)
out.csv.imp <- paste0(
mod.dir, "/", var, "_vs_pooled_imp_long.csv"
)
if (!all(file.exists(out.csv.pred.obs),
file.exists(out.csv.tune),
file.exists(out.csv.imp))){
# run the model and get outputs
mod <- tuned.model(mod.df)
df.pred.obs <- mod$df.pred.obs
df.tune <- mod$df.tune
df.imp <- mod$df.imp
# write outputs to files
write.csv(df.pred.obs, out.csv.pred.obs, row.names = F)
write.csv(df.tune, out.csv.tune, row.names = F)
write.csv(df.imp, out.csv.imp, row.names = F)
} else {
# if they've already been run, read in pred vs. obs data
df.pred.obs <- read.csv(out.csv.pred.obs)
}
# get performance metrics from pred vs. obs data
perf.mets.long.temp <- prd.obs.plot(df.pred.obs, var, plot = F)
# compile results
perf.mets.long.temp$type <- "pooled.long"
if (var == vars[1]){
perf.mets.long <- perf.mets.long.temp
} else {
perf.mets.long <- rbind(perf.mets.long, perf.mets.long.temp)
}
#-----------dif.shrt
# read in the data
df.pred.obs <- read.csv(
paste0(mod.dir, "/", var, "_vs_pooled_pred_obs_shrt.csv")
)
# remove obs
df.pred.obs$obs <- NULL
# separate back out the pre- and post-fire plots
df.pred.obs.pre <- df.pred.obs[grep("_pre", df.pred.obs$plot),]
df.pred.obs.pst <- df.pred.obs[grep("_pst", df.pred.obs$plot),]
# add "_pre" and "_pst" to prediction column names
colnames(df.pred.obs.pre)[2] <- paste0(colnames(df.pred.obs.pre)[2], "_pre")
colnames(df.pred.obs.pst)[2] <- paste0(colnames(df.pred.obs.pst)[2], "_pst")
# remove "_pre" and "_pst" from plot names
df.pred.obs.pre$plot <- sub("_pre", "", df.pred.obs.pre$plot)
df.pred.obs.pst$plot <- sub("_pst", "", df.pred.obs.pst$plot)
# merge pre and post fire predictions
y.pred <- merge(df.pred.obs.pre, df.pred.obs.pst, by = "plot", all = T)
# calculate difference and remove pre and pst fields
y.pred$prd <- y.pred$prd_pre - y.pred$prd_pst
y.pred$prd_pre <- NULL
y.pred$prd_pst <- NULL
# define variable of interest
var.dif <- paste0(var, "_dif")
# get observed differences
y.obs <- fuel.plots[,c("plot", var.dif)] |> as.data.frame()
colnames(y.obs)[2] <- "obs"
# merge the two
df.pred.obs <- merge(y.obs, y.pred)
df.pred.obs <- df.pred.obs[,c("plot", "obs", "prd")]
# write to file
write.csv(
df.pred.obs,
paste0(mod.dir, "/", var.dif, "_vs_pooled_pred_obs_shrt.csv"),
row.names = F
)
# get performance metrics
perf.mets.dif.shrt.temp <- prd.obs.plot(df.pred.obs, var.dif, plot = F)
# compile results
perf.mets.dif.shrt.temp$type <- "pooled.dif.shrt"
if (var == vars[1]){
perf.mets.dif.shrt <- perf.mets.dif.shrt.temp
} else {
perf.mets.dif.shrt <- rbind(perf.mets.dif.shrt, perf.mets.dif.shrt.temp)
}
#-----------dif.long
# read in the data
df.pred.obs <- read.csv(
paste0(mod.dir, "/", var, "_vs_pooled_pred_obs_long.csv")
)
# remove obs
df.pred.obs$obs <- NULL
# separate back out the pre- and post-fire plots
df.pred.obs.pre <- df.pred.obs[grep("_pre", df.pred.obs$plot),]
df.pred.obs.pst <- df.pred.obs[grep("_pst", df.pred.obs$plot),]
# add "_pre" and "_pst" to prediction column names
colnames(df.pred.obs.pre)[2] <- paste0(colnames(df.pred.obs.pre)[2], "_pre")
colnames(df.pred.obs.pst)[2] <- paste0(colnames(df.pred.obs.pst)[2], "_pst")
# remove "_pre" and "_pst" from plot names
df.pred.obs.pre$plot <- sub("_pre", "", df.pred.obs.pre$plot)
df.pred.obs.pst$plot <- sub("_pst", "", df.pred.obs.pst$plot)
# merge pre and post fire predictions
y.pred <- merge(df.pred.obs.pre, df.pred.obs.pst, by = "plot", all = T)
# calculate difference and remove pre and pst fields
y.pred$prd <- y.pred$prd_pre - y.pred$prd_pst
y.pred$prd_pre <- NULL
y.pred$prd_pst <- NULL
# define variable of interest
var.dif <- paste0(var, "_dif")
# get observed differences
y.obs <- fuel.plots[,c("plot", var.dif)] |> as.data.frame()
colnames(y.obs)[2] <- "obs"
# merge the two
df.pred.obs <- merge(y.obs, y.pred)
df.pred.obs <- df.pred.obs[,c("plot", "obs", "prd")]
# write to file
write.csv(
df.pred.obs,
paste0(mod.dir, "/", var.dif, "_vs_pooled_pred_obs_long.csv"),
row.names = F
)
# get performance metrics
perf.mets.dif.long.temp <- prd.obs.plot(df.pred.obs, var.dif, plot = F)
# compile results
perf.mets.dif.long.temp$type <- "pooled.dif.long"
if (var == vars[1]){
perf.mets.dif.long <- perf.mets.dif.long.temp
} else {
perf.mets.dif.long <- rbind(perf.mets.dif.long, perf.mets.dif.long.temp)
}
}
# write to csvs
write.csv(perf.mets.shrt, out.perf.mets.shrt.csv, row.names = F)
write.csv(perf.mets.long, out.perf.mets.long.csv, row.names = F)
write.csv(perf.mets.dif.shrt, out.perf.mets.dif.shrt.csv, row.names = F)
write.csv(perf.mets.dif.long, out.perf.mets.dif.long.csv, row.names = F)
}The results below are focused on comparison between the different modeling approaches. They represent models built with the short-term (1 year) post-fire lidar data, as this should represent the “best-case” scenario for fuel consumption mapping, given the temporal proximity to the fire events. It is very clear that PDC is dominant among the different modeling approaches, with higher R2 and lower %RMSE values nearly for every fuel bed. There is only one exception to this (tree), where SPAC proved superior. Among the two different PDC approaches, it appears that including pre-fire densities yielded improved predictive power in several instances, though the difference-only models still performed best in several others. Generally, pre-post and pooled outperformed SPAC.
# read in the data
perf.mets.pdc <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_pdc_bw_opt_perf_mets.csv")
)
perf.mets.spac <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_spac_perf_mets.csv")
)
perf.mets.prepost <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_prepost_shrt_perf_mets_dif_shrt.csv")
)
perf.mets.pooled <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_pooled_perf_mets_dif_shrt.csv")
)
# subset dfs as needed
perf.mets.pdc <- perf.mets.pdc[
perf.mets.pdc$bw == 0.01 & grepl("shrt", perf.mets.pdc$type),
]
perf.mets.spac <- perf.mets.spac[
grep("shrt", perf.mets.spac$type),
]
# split out the pre
perf.mets.pdc.preanddif <-
perf.mets.pdc[perf.mets.pdc$type == "pre.and.dif.shrt",]
perf.mets.pdc.difonly <-
perf.mets.pdc[perf.mets.pdc$type == "dif.shrt",]
perf.mets.spac.preanddif <-
perf.mets.spac[perf.mets.spac$type == "pre.and.spac.shrt",]
perf.mets.spac.difonly <-
perf.mets.spac[perf.mets.spac$type == "spac.shrt",]
# combine them to create rbind table
perf.mets.rbind <- bind_rows(
perf.mets.pdc.preanddif,
perf.mets.pdc.difonly,
perf.mets.spac.preanddif,
perf.mets.spac.difonly,
perf.mets.prepost,
perf.mets.pooled
)
# remove "_dif" from variable names
perf.mets.rbind$var <- sub("_dif", "", perf.mets.rbind$var)
# define response variables of interest
vars <- c("hr1", "hr10", "hr100", "hr1000", "herb", "shrub", "seed", "llm",
"duff", "shrbsd", "woody", "nontre", "tree", "acf", "tf", "af")
# simplify tables
perf.mets.pdc.preanddif <- perf.mets.pdc.preanddif[,c("var", "r2", "rrmse")]
perf.mets.pdc.difonly <- perf.mets.pdc.difonly[,c("var", "r2", "rrmse")]
perf.mets.spac.preanddif <- perf.mets.spac.preanddif[,c("var", "r2", "rrmse")]
perf.mets.spac.difonly <- perf.mets.spac.difonly[,c("var", "r2", "rrmse")]
perf.mets.prepost <- perf.mets.prepost[,c("var", "r2", "rrmse")]
perf.mets.pooled <- perf.mets.pooled[,c("var", "r2", "rrmse")]
# add method to r2 and rrmse colnames
colnames(perf.mets.pdc.preanddif)[2:3] <-
paste0(colnames(perf.mets.pdc.preanddif)[2:3], ".pdc.preanddif")
colnames(perf.mets.pdc.difonly)[2:3] <-
paste0(colnames(perf.mets.pdc.difonly)[2:3], ".pdc.difonly")
colnames(perf.mets.spac.preanddif)[2:3] <-
paste0(colnames(perf.mets.spac.preanddif)[2:3], ".spac.preanddif")
colnames(perf.mets.spac.difonly)[2:3] <-
paste0(colnames(perf.mets.spac.difonly)[2:3], ".spac.difonly")
colnames(perf.mets.prepost)[2:3] <-
paste0(colnames(perf.mets.prepost)[2:3], ".prepost")
colnames(perf.mets.pooled)[2:3] <-
paste0(colnames(perf.mets.pooled)[2:3], ".pooled")
# merge them to create "cbind" table
perf.mets.cbind <- merge(perf.mets.pdc.preanddif, perf.mets.pdc.difonly) |>
merge(perf.mets.spac.preanddif) |>
merge(perf.mets.spac.difonly) |>
merge(perf.mets.prepost) |>
merge(perf.mets.pooled)
# remove "_dif" from variable names
perf.mets.cbind$var <- sub("_dif", "", perf.mets.cbind$var)
# sort
perf.mets.cbind <- perf.mets.cbind[match(vars, perf.mets.cbind$var),]
# define colors
cols <- brewer.pal(10, "Paired")[c(1,2,5,6,9,10)]
# set up the plot
par(mfrow = c(1,2), las = 1)
# create empty plot
x <- rep(100,length(vars))
names(x) <- vars
r2.cols <- perf.mets.cbind[,grep("r2", colnames(perf.mets.cbind))]
dotchart(x, pch = 16, xlab = expression(R^2),
xlim = range(r2.cols))
# loop though variables
for (i in seq(1,length(vars))){
var <- vars[i]
points(x = perf.mets.cbind$r2.prepost[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[1])
points(x = perf.mets.cbind$r2.pooled[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[2])
points(x = perf.mets.cbind$r2.spac.difonly[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[3])
points(x = perf.mets.cbind$r2.spac.preanddif[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[4])
points(x = perf.mets.cbind$r2.pdc.difonly[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[5])
points(x = perf.mets.cbind$r2.pdc.preanddif[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[6])
}
# create empty plot
y <- rep(10000, length(vars))
names(y) <- vars
rrmse.cols <- perf.mets.cbind[
!perf.mets.cbind$var %in% c("herb", "shrub"),
grep("rrmse", colnames(perf.mets.cbind))
]
dotchart(y, pch = 16, xlab = expression(symbol("%")*RMSE),
xlim = range(rrmse.cols))
# loop though variables
for (i in seq(1,length(vars))){
var <- vars[i]
points(x = perf.mets.cbind$rrmse.prepost[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[1])
points(x = perf.mets.cbind$rrmse.pooled[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[2])
points(x = perf.mets.cbind$rrmse.spac.difonly[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[3])
points(x = perf.mets.cbind$rrmse.spac.preanddif[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[4])
points(x = perf.mets.cbind$rrmse.pdc.difonly[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[5])
points(x = perf.mets.cbind$rrmse.pdc.preanddif[perf.mets.cbind$var == var],
y = i, pch = 16, col = cols[6])
}
# add legend
legend(
"topright",
legend = c("PrePost", "Pooled", "SPAC Diff Only", "SPAC Pre & Diff",
"PDC Diff Only", "PDC Pre & Diff"),
pch = 16, col = cols
)
# figure out which approach performed best for each fuelbed
fn <- function(x){
m <- which.max(x)
cn <- colnames(r2.cols)[m]
var <- sub("r2\\.", "", cn)
return(var)
}
r2.best <- apply(r2.cols, 1, fn)
fn <- function(x){
m <- which.max(x)
cn <- colnames(rrmse.cols)[m]
var <- sub("rrmse\\.", "", cn)
return(var)
}
rrmse.best <- apply(r2.cols, 1, fn)
df.best <- data.frame(var = perf.mets.cbind$var,
r2.best = r2.best,
rrmse.best = rrmse.best)
lut <- data.frame(
abbv = c(
"pdc.preanddif",
"pdc.difonly",
"spac.preanddif",
"spac.difonly",
"prepost",
"pooled"
),
full = c(
"PDC Pre & Diff",
"PDC Diff Only",
"SPAC Pre & Diff",
"SPAC Diff Only",
"Pre-Post",
"Pooled"
)
)
df.best <- merge(df.best, lut, by.x = "r2.best", by.y = "abbv")
df.best$r2.best <- NULL
colnames(df.best)[colnames(df.best) == "full"] <- "r2.best"
df.best <- merge(df.best, lut, by.x = "rrmse.best", by.y = "abbv")
df.best$rrmse.best <- NULL
colnames(df.best)[colnames(df.best) == "full"] <- "rrmse.best"
df.best <- df.best[match(vars, df.best$var),]
kable(df.best,
col.names = c("Fuel Metric", "Best Model by R2", "Best Model by %RMSE"))| Fuel Metric | Best Model by R2 | Best Model by %RMSE | |
|---|---|---|---|
| 1 | hr1 | PDC Diff Only | PDC Diff Only |
| 2 | hr10 | PDC Diff Only | PDC Diff Only |
| 7 | hr100 | PDC Pre & Diff | PDC Pre & Diff |
| 10 | hr1000 | PDC Pre & Diff | PDC Pre & Diff |
| 8 | herb | PDC Pre & Diff | PDC Pre & Diff |
| 4 | shrub | PDC Diff Only | PDC Diff Only |
| 9 | seed | PDC Pre & Diff | PDC Pre & Diff |
| 13 | llm | PDC Pre & Diff | PDC Pre & Diff |
| 11 | duff | PDC Pre & Diff | PDC Pre & Diff |
| 5 | shrbsd | PDC Diff Only | PDC Diff Only |
| 12 | woody | PDC Pre & Diff | PDC Pre & Diff |
| 15 | nontre | PDC Pre & Diff | PDC Pre & Diff |
| 16 | tree | SPAC Diff Only | SPAC Diff Only |
| 6 | acf | PDC Diff Only | PDC Diff Only |
| 14 | tf | PDC Pre & Diff | PDC Pre & Diff |
| 3 | af | PDC Diff Only | PDC Diff Only |
# plot out their performance as boxplots
par(mfrow = c(1,2), mar = c(7,3,1,1), las = 1)
perf.mets.rbind$type <- factor(
perf.mets.rbind$type,
levels = c(
"prepost.dif.shrt", "pooled.dif.shrt",
"spac.shrt", "pre.and.spac.shrt",
"dif.shrt", "pre.and.dif.shrt"
)
)
bp <- boxplot(r2 ~ type, data = perf.mets.rbind, col = cols,
xaxt = "n", ylab = expression(R^2),
xlab = NA, outline = F)
tick <- seq_along(bp$names)
axis(1, at = tick, labels = F)
text(x = tick, y = par("usr")[3] - 0.025 * diff(par("usr")[3:4]),
labels = c(
"Pre-Post", "Pooled", "SPAC Diff Only", "SPAC Pre & Diff",
"PDC Diff Only", "PDC Pre & Diff"
),
srt = 45, xpd = T, adj = c(1,1))
bp <- boxplot(rrmse ~ type, data = perf.mets.rbind, col = cols,
xaxt = "n",
ylab = expression(symbol("%")*RMSE), xlab = NA, outline = F)
tick <- seq_along(bp$names)
axis(1, at = tick, labels = F)
text(x = tick, y = par("usr")[3] - 0.025 * diff(par("usr")[3:4]),
labels = c(
"Pre-Post", "Pooled", "SPAC Diff Only", "SPAC Pre & Diff",
"PDC Diff Only", "PDC Pre & Diff"
),
srt = 45, xpd = T, adj = c(1,1))
I ran all of the models above for both the short-term (1 year) post-fire lidar data as well as the long-term (2-3 year) post-fire data to assess how time-since-fire (regeneration, snag fall, etc.) may have affected model performance. Interestingly, the models performed very similarly. PDC and SPAC both tended to perform slightly better on the short-term post-fire data, but pre-post and pooled actually did the opposite.
#---------------shrt
# read in the data
perf.mets.pdc <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_pdc_bw_opt_perf_mets.csv")
)
perf.mets.spac <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_spac_perf_mets.csv")
)
perf.mets.prepost <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_prepost_shrt_perf_mets_dif_shrt.csv")
)
perf.mets.pooled <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_pooled_perf_mets_dif_shrt.csv")
)
# subset dfs as needed
perf.mets.pdc <- perf.mets.pdc[
perf.mets.pdc$bw == 0.01 & grepl("shrt", perf.mets.pdc$type),
]
perf.mets.spac <- perf.mets.spac[
grep("shrt", perf.mets.spac$type),
]
# split out the pre
perf.mets.pdc.preanddif <-
perf.mets.pdc[perf.mets.pdc$type == "pre.and.dif.shrt",]
perf.mets.pdc.difonly <-
perf.mets.pdc[perf.mets.pdc$type == "dif.shrt",]
perf.mets.spac.preanddif <-
perf.mets.spac[perf.mets.spac$type == "pre.and.spac.shrt",]
perf.mets.spac.difonly <-
perf.mets.spac[perf.mets.spac$type == "spac.shrt",]
# combine them to create rbind table
perf.mets.rbind.shrt <- bind_rows(
perf.mets.pdc.preanddif,
perf.mets.pdc.difonly,
perf.mets.spac.preanddif,
perf.mets.spac.difonly,
perf.mets.prepost,
perf.mets.pooled
)
# remove "_dif" from variable names
perf.mets.rbind.shrt$var <- sub("_dif", "", perf.mets.rbind.shrt$var)
# simplify
perf.mets.rbind.shrt <- perf.mets.rbind.shrt[,c("type", "var", "r2", "rrmse")]
# add "shrt" to r2 and rrmse columns
colnames(perf.mets.rbind.shrt)[3:4] <- paste0(
colnames(perf.mets.rbind.shrt)[3:4], ".shrt"
)
# remove "shrt" from types
perf.mets.rbind.shrt$type <- sub("\\.shrt", "", perf.mets.rbind.shrt$type)
#-------------long
# read in the data
perf.mets.pdc <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_pdc_bw_opt_perf_mets.csv")
)
perf.mets.spac <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_spac_perf_mets.csv")
)
perf.mets.prepost <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_prepost_long_perf_mets_dif_long.csv")
)
perf.mets.pooled <- read.csv(
paste0(mod.dir, "/fuel_cons_vs_pooled_perf_mets_dif_long.csv")
)
# subset dfs as needed
perf.mets.pdc <- perf.mets.pdc[
perf.mets.pdc$bw == 0.01 & grepl("long", perf.mets.pdc$type),
]
perf.mets.spac <- perf.mets.spac[
grep("long", perf.mets.spac$type),
]
# split out the pre
perf.mets.pdc.preanddif <-
perf.mets.pdc[perf.mets.pdc$type == "pre.and.dif.long",]
perf.mets.pdc.difonly <-
perf.mets.pdc[perf.mets.pdc$type == "dif.long",]
perf.mets.spac.preanddif <-
perf.mets.spac[perf.mets.spac$type == "pre.and.spac.long",]
perf.mets.spac.difonly <-
perf.mets.spac[perf.mets.spac$type == "spac.long",]
# combine them to create rbind table
perf.mets.rbind.long <- bind_rows(
perf.mets.pdc.preanddif,
perf.mets.pdc.difonly,
perf.mets.spac.preanddif,
perf.mets.spac.difonly,
perf.mets.prepost,
perf.mets.pooled
)
# remove "_dif" from variable names
perf.mets.rbind.long$var <- sub("_dif", "", perf.mets.rbind.long$var)
# simplify
perf.mets.rbind.long <- perf.mets.rbind.long[,c("type", "var", "r2", "rrmse")]
# add "long" to r2 and rrmse columns
colnames(perf.mets.rbind.long)[3:4] <- paste0(
colnames(perf.mets.rbind.long)[3:4], ".long"
)
# remove "long" from types
perf.mets.rbind.long$type <- sub("\\.long", "", perf.mets.rbind.long$type)
#-------------merge
# merge data
perf.mets.rbind.merge <- merge(perf.mets.rbind.shrt, perf.mets.rbind.long)
# get differences
perf.mets.rbind.merge$r2.diff <- perf.mets.rbind.merge$r2.shrt -
perf.mets.rbind.merge$r2.long
perf.mets.rbind.merge$rrmse.diff <- perf.mets.rbind.merge$rrmse.shrt -
perf.mets.rbind.merge$rrmse.long
# define response variables of interest
vars <- c("hr1", "hr10", "hr100", "hr1000", "herb", "shrub", "seed", "llm",
"duff", "shrbsd", "woody", "nontre", "tree", "acf", "tf", "af")
# plot out their performance as boxplots
par(mfrow = c(1,2), mar = c(7,5,1,1), las = 1)
perf.mets.rbind.merge$type <- factor(
perf.mets.rbind.merge$type,
levels = c(
"prepost.dif", "pooled.dif",
"spac", "pre.and.spac",
"dif", "pre.and.dif"
)
)
bp <- boxplot(r2.diff ~ type, data = perf.mets.rbind.merge, col = cols,
xaxt = "n", ylab = expression({R^2}["short"]-{R^2}["long"]),
xlab = NA, outline = F)
abline(h = 0, lty = 3)
tick <- seq_along(bp$names)
axis(1, at = tick, labels = F)
text(x = tick, y = par("usr")[3] - 0.025 * diff(par("usr")[3:4]),
labels = c(
"Pre-Post", "Pooled", "SPAC Diff Only", "SPAC Pre & Diff",
"PDC Diff Only", "PDC Pre & Diff"
),
srt = 45, xpd = T, adj = c(1,1))
bp <- boxplot(rrmse.diff ~ type, data = perf.mets.rbind.merge, col = cols,
xaxt = "n",
ylab = expression(
{symbol("%")*RMSE}["short"]-{symbol("%")*RMSE}["long"]
),
xlab = NA, outline = F)
abline(h = 0, lty = 3)
tick <- seq_along(bp$names)
axis(1, at = tick, labels = F)
text(x = tick, y = par("usr")[3] - 0.025 * diff(par("usr")[3:4]),
labels = c(
"Pre-Post", "Pooled", "SPAC Diff Only", "SPAC Pre & Diff",
"PDC Diff Only", "PDC Pre & Diff"
),
srt = 45, xpd = T, adj = c(1,1))
To my mind, there is maybe one or two more steps:
Anything else? If not, then we should probably start talking about journals for this type of work. I could see it landing in a remote sensing-focused journal, given the amount of technical lidar stuff going on. But, I could also see it being appropriate for a more fire-focused journal. In any case, it would be nice to have an idea upfront so I can emphasize the paper one way or the other as I start to write it. Another question would be who to include as coauthors, but we’ve got plenty of time for that.