vignettes/firebehavioR.R

In EcoFire/firebehavioR: Prediction of Wildland Fire Behavior and Hazard

## ----setup, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ---- eval=FALSE---------------------------------------------------------
#  devtools::install_github("EcoFire/firebehavioR")

## ---- eval=FALSE---------------------------------------------------------
#  library("firebehavioR")

## ---- echo=FALSE---------------------------------------------------------
library("firebehavioR")

## ------------------------------------------------------------------------
data("coForest")
print(coForest)

## ------------------------------------------------------------------------
ex1 = cfis(fsg = coForest$cbh_m, u10 = 20, effm = 6, sfc = coForest$sfl_kgm2 * 10, 
    cbd = coForest$cbd_kgm3, id = 1)
print(ex1)

## ------------------------------------------------------------------------
ex2 = aggregate(x = ex1$cROS, by = list(treatmentStatus = coForest$status), FUN = mean)
print(ex2)

## ------------------------------------------------------------------------
ex3 = coForest[rep(seq_len(nrow(coForest)), 11), ]
ex3$u10 = sort(rep(10:20, 14))
head(ex3)

ex3$type = cfis(ex3$cbh_m, ex3$u10, 6, ex3$sfl_kgm2 * 10, ex3$cbd_kgm3, 1)$type
table(u10 = ex3$u10, ex3$type, ex3$status)

## ------------------------------------------------------------------------
ex4 = cfis(coForest[1, ]$cbh_m, u10 = 14, effm = 6, sfc = coForest[1, ]$sfl_kgm2, 
    cbd = coForest[1, ]$cbd_kgm3, id = 0:10)
ex4$id = 0:10

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Requisite conditions for successful spot fire ignition given the projected crown fire rate of spread."----
library("ggplot2")

ggplot(ex4, aes(sepDist, id)) + geom_ribbon(aes(ymin = pmin(sepDist, id), ymax = 0), 
    fill = "red", col = "red", alpha = 0.5) + geom_line(size = 1.5) + labs(x = "Critical seperation distance (m)", 
    y = "Ignition delay (min)") + annotate("text", x = 450, y = 3, label = "Spot fire ignition zone") + 
    theme_classic()

## ------------------------------------------------------------------------
data(fuelModels, fuelMoisture)

fuelModels['A10',]
exampSurfFuel = fuelModels['A10',]

fuelMoisture['D1L1',]
exampFuelMoisture = fuelMoisture['D1L1',]

exampCrownFuel = data.frame(
  CBD = coForest$cbd_kgm3,
  FMC = 100,
  CBH = coForest$cbh_m,
  CFL = coForest$cfl_kgm2
)

exampEnviro = data.frame(
  slope = 10,
  windspeed = 40,
  direction = 0,
  waf = 0.2
)

ex5 = rothermel(sapply(exampSurfFuel, rep, 14), sapply(exampFuelMoisture, rep, 14), exampCrownFuel, sapply(exampEnviro,rep,14))

## ------------------------------------------------------------------------
head(ex5$fireBehavior)

## ------------------------------------------------------------------------
head(ex5$detailSurface)
head(ex5$detailCrown)

## ------------------------------------------------------------------------
head(ex5$critInit)
head(ex5$critActive)
head(ex5$critCess)

## ------------------------------------------------------------------------
exampEnviro = data.frame(
  slope = 30,
  windspeed = 60,
  direction = 0,
  waf = 0.2
)

exampCrownFuel = data.frame(
  CBD = truncnorm::rtruncnorm(1000, a=0, coForest$cbd_kgm3[1], sd = coForest$cbd_kgm3[1]*.2),
  FMC = rep(100,1000),
  CBH = truncnorm::rtruncnorm(1000, a=0, coForest$cbh_m[1], sd = coForest$cbh_m[1]*.2),
  CFL = truncnorm::rtruncnorm(1000, a=0, coForest$cfl_kgm2[1], sd = coForest$cfl_kgm2[1]*.2)
)

ex6 = rep(NA,200)
for (i in 1:200){
  ex6[i] = rothermel(exampSurfFuel, exampFuelMoisture, exampCrownFuel[i,], exampEnviro, cfbForm = "sr")$fireBehavior[,3]  
}

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Distribution of hourly fuel moistures predicted by three methods using the ffm() function."----
ggplot(data.frame(ex6), aes(ex6)) + geom_density(fill='blue',alpha = 0.5,adjust = 1/5) + 
    xlab("Rate of spread (m/s)") + theme_classic() + coord_cartesian(expand=FALSE)

## ------------------------------------------------------------------------
fuelModels['A5',]
exampSurfFuel = fuelModels['A5',]

fuelMoisture['D1L1',]
exampFuelMoisture = fuelMoisture['D1L1',]

exampCrownFuel = data.frame(
  CBD = Inf,  FMC = 100,
  CBH = Inf,  CFL = Inf
)

exampEnviro = data.frame(
  slope = 40,  windspeed = 60,
  direction = 0, waf = 0.3
)

ex7 = expand.grid(winddir=seq(-180,180,10),slope = seq(0,30,5))
ex7$ros = NA
for (i in 1:259){
  exampEnviro$direction = ex7$winddir[i]
  exampEnviro$slope     = ex7$slope[i]
  ex7$ros[i] = rothermel(exampSurfFuel, exampFuelMoisture, exampCrownFuel, exampEnviro)$fireBehavior[,3]  
}

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Distribution of hourly fuel moistures predicted by three methods using the ffm() function."----
ggplot(ex7,aes(winddir,ros,color=slope,group=slope)) +geom_path() + theme_classic() + labs(x='Wind direction (deg)',y='Rate of spread (m/s)')

## ------------------------------------------------------------------------
ba = c(10, 15)
ht = c(12, 20)
tph = c(100, 300)
type = c("df", "pp", "mc","lp")
cbind(type,canFuel(ba, ht, tph, type))


## ------------------------------------------------------------------------
data(rrRAWS)
print(head(rrRAWS))

ff = rbind(
data.frame(ffm = ffm("simard",rrRAWS$rh, rrRAWS$temp_c)$fm1hr,method="simard"),
data.frame(ffm = ffm("wagner",rrRAWS$rh, rrRAWS$temp_c)$fm1hr,method="wagner"),
data.frame(ffm = ffm("anderson",rrRAWS$rh, rrRAWS$temp_c)$fm1hr,method="anderson")
)
ff$dateTime = rep(rrRAWS$dateTime,3)

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Distribution of hourly fuel moistures predicted by three methods using the ffm() function."----
ggplot(ff, aes(ffm, color = method, fill = method)) + geom_density(alpha = 0.5) + 
    xlab("Fine fuel moisture (%)") + theme_classic()

## ------------------------------------------------------------------------
rrRAWS_ma = subset(rrRAWS, rh <= 70 & temp_c >= 10 & rh > 42.5 - 1.25 * temp_c & 
    rh < 94.5 - 1.35 * temp_c)
ff.ma = data.frame(ffm = ffm("mcarthur", rrRAWS_ma$rh, rrRAWS_ma$temp_c)$fm1hr, method = "mcarthur", 
    dateTime = rrRAWS_ma$dateTime)
ff = rbind(data.frame(ffm = ffm("simard", rrRAWS$rh, rrRAWS$temp_c)$fm1hr, method = "simard"), 
    data.frame(ffm = ffm("wagner", rrRAWS$rh, rrRAWS$temp_c)$fm1hr, method = "wagner"), 
    data.frame(ffm = ffm("anderson", rrRAWS$rh, rrRAWS$temp_c)$fm1hr, method = "anderson"))
ff$dateTime = rep(rrRAWS$dateTime, 3)
ff = rbind(ff, ff.ma)
ff$dateTime <- strptime(ff$dateTime, "%m/%d/%Y %H:%M")

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Hourly fuel moistures predicted by four methods using the ffm() function."----
ffm.plot = ggplot(ff, aes(x = dateTime, y = ffm, color = method)) + geom_smooth(span = 0.1, 
    method = "loess", se = F) + theme_classic() + labs(x = "Time", y = "Fuel moisture (%)")
print(ffm.plot)

## ------------------------------------------------------------------------
rh = rrRAWS$rh
temp = rrRAWS$temp_c
month = as.numeric(format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%m"))
hour = as.numeric(format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%H"))
ff.fbo = data.frame(ffm = ffm(method = "fbo", rh, temp, month, hour, asp = "N", slp = 10, 
    bla = "l", shade = "n")$fm1hr, method = "FBO", dateTime = strptime(rrRAWS$dateTime, 
    "%m/%d/%Y %H:%M"))
ff = rbind(ff, ff.fbo)

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Hourly fuel moistures predicted by method = 'fbo', in blue, relative to other methods, in grayscale."----
ffm.plot + scale_colour_grey() + geom_smooth(data = ff.fbo, aes(x = dateTime, y = ffm), 
    span = 0.1, method = "loess", se = F, color = "blue") + guides(color = FALSE) + 
    labs(x = "Time", y = "Fuel moisture (%)")

## ------------------------------------------------------------------------
waf(forestHt = 10, cr = 40, cc = 40)

## ------------------------------------------------------------------------
waf(forestHt = 10, sheltered = "y")

## ------------------------------------------------------------------------
waf(fuelDepth = 1)

## ------------------------------------------------------------------------
waf(fuelDepth = 1, forestHt = 10, cr = 10, cc = 10)

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Fire characteristics chart.",warning = FALSE----
fireChart("",100,10) +xlim(0,12000)


## ------------------------------------------------------------------------
rrRAWS.daily =   rrRAWS[format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%H:%M")=="14:35",]
indices = fireIndex(temp=rrRAWS.daily$temp_c, u= rrRAWS.daily$windSpeed_kmh, rh = rrRAWS.daily$rh)
normalize=function(x,low=0,high=1){low+(x-min(x,na.rm=T))*(high-low)/(max(x,na.rm=T)-min(x,na.rm=T))}
indices = data.frame(sapply(indices,normalize),  Date = strptime(rrRAWS.daily$dateTime, "%m/%d/%Y %H:%M"))
indices = setNames(reshape2::melt(indices,id="Date"),c("Date","Index","Value"))


## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Continuous fire index values applied to the rrRAWS dataset.",warning = FALSE----
ggplot(indices,aes(Date,Value,group=Index,color=Index))+geom_smooth(span = .1,  method = "loess", se = F) +theme_classic()+coord_cartesian(expand=F)

## ------------------------------------------------------------------------
data(rrRAWS)
daily.precip = rrRAWS
daily.precip$Date = strptime(daily.precip$dateTime, "%m/%d/%Y")
daily.precip = setNames(aggregate(daily.precip$precip_mm, by = list(as.character(daily.precip$Date)), 
    FUN = sum), c("Date", "DailyPrecip"))
rrRAWS.daily = rrRAWS[format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%H:%M") == 
    "14:35", ]
rrRAWS.daily$DailyPrecip = daily.precip$DailyPrecip
indices = fireIndexKBDI(temp = rrRAWS.daily$temp_c, precip = rrRAWS.daily$DailyPrecip, 
    map = 610, rh = rrRAWS.daily$rh, u = rrRAWS.daily$windSpeed_kmh)
indices = data.frame(sapply(indices,normalize), Date = strptime(rrRAWS.daily$dateTime, 
    "%m/%d/%Y %H:%M"))
indices = setNames(reshape2::melt(indices, id = "Date"), c("Date", "Index", "Value"))

## ----fig.height = 3, fig.width = 7, fig.align = "center", fig.cap = "Static fire index values applied to the rrRAWS dataset.",warning = FALSE----
ggplot(indices, aes(Date, Value, group = Index, color = Index)) + geom_smooth(span = 0.1, 
    method = "loess", se = F) + theme_classic() + coord_cartesian(expand = F)

## ---- echo=FALSE, results="asis"-----------------------------------------
table = xtable::xtable(data.frame(Method = c("KBDI", "Drought factor", "Forest Mark 5", 
    "Fosberg-KBDI", "Fuel moisture-KBDI", "Nesterov", "Nesterov-Modified", "Zdenko"), 
    Inputs = c("temp, precip, map", "temp, precip, map", "temp, precip, map, u, rh", 
        "temp, precip, map, u, rh", "temp, precip, map, u, rh", "temp, precip, rh", 
        "temp, precip, rh", "temp, precip, rh")))
print(table, type = "html", sanitize.text.function = identity)

## ---- echo=FALSE, results="asis"-----------------------------------------
table = xtable::xtable(data.frame(Variable = c("Fuel stratum gap", "10-m Open wind speed", 
    "Fuel moisture content", "Canopy bulk density", "Seperation distance", "Basal area", 
    "Average stand tree heights", "Trees per hectare", "Relative humidity", "Temperature", 
    "Month", "Hour", "Topographic aspect", "Topographic slope", "Heat per unit area", 
    "Fire rate of spread", "surface fuel load", "Surface area to volume", "Fuelbed depth", 
    "Moisture of extinction", "Heat content", "Canopy fuel load", "Wind direction", 
    "Wind adjustment factor", "Crown ratio", "Canopy cover"), Units = c("m", "km/hr", 
    "%", "kg/m^3^", "m", "m^2^/ha", "m", "trees/ha", "%", "degrees C", "Month of the year (1-12)", 
    "Hour of day (1-24)", "N,S,W, or E", "%", "kJ/m^2^", "m/min", "Mg/ha", "m^2/^m^3^", 
    "cm", "%", "J/g", "kg/m^2^", "0-360", "Ratio of 20-ft open wind speed to midflame wind speed", 
    "%", "%")))
print(table, type = "html", sanitize.text.function = identity)