| SOCR ≫ | DSPA ≫ | DSPA2 Topics ≫ |
DSPA2: Data Science and Predictive Analytics (UMich HS650)
Appendix 19A: Spatiotemporal Wildfire SDE Modeling
SOCR (Ivo Dinov and SOCR)
April 23, 2025
This DSPA Appendix shows examples of dynamic stochastic differential equation (SDE) models of wildfire. The main objective is to develop an advanced generative stochastic differential equation (SDE) model proposing a new approach of prediction of fire propagation.
Review the earlier Volterra-Lotka preditor-pray model and Hodgkin-Huxley (HH) model of neuronal neuron action potential propagation, which provide an oversimplified motivation for building a more advanced SDE cancer-model, which includes many controllable parameters accounting for the univariate distributions of dozens of observable clinical variables (biomarkers) and multivariate relations. The model (hyper)parameters can be tuned by Bayesian optimization using very sparsely sampled longitudinal data.
I’ll help you develop a stochastic differential equation (SDE) model for wildfire propagation based on the cancer treatment SDE model in the document. I’ll create an Rmd template with all the necessary components including model development, simulations, and visualizations.
1 Wildfire Propagation Model Using Stochastic Differential Equations
This DSPA Appendix develops an advanced generative stochastic differential equation (SDE) model for predicting wildfire propagation across landscapes. The model incorporates spatial and temporal dynamics, environmental factors, and stochastic elements that influence fire behavior. Similar to the DSPA cancer treatment SDE model, this wildfire propagation model will track multiple interacting variables
Fire Intensity (FI) - analogous to tumor Local Control (LC)
Fuel Consumption (FC) - similar to side effects (RP2)
Environmental Factors (E) - similar to biomarkers vector
The model will account for complex interactions between these variables while incorporating stochastic elements to represent unpredictable factors like wind gusts, terrain variations, and fuel heterogeneity.
2 SDE Model Design
Let’s define our state variables as
\[\mathbf{X}(t) = \begin{pmatrix} FI(t) \\ FC(t) \\ \mathbf{E}(t) \end{pmatrix},\]
where, \(FI(t)\) measures fire intensity at time \(t\), \(FC(t)\) measures the rate of fuel consumption, and \(\mathbf{E}(t)\) is a vector of environmental factors, e.g., wind speed, humidity, temperature, slope. We model the time evolution of hte wildfire as an SDE system
\[d \mathbf{X}(t) = \mathbf{f}(\mathbf{X}(t), \mathbf{c}(t), \boldsymbol{\theta})dt + \mathbf{g}(\mathbf{X}(t), \mathbf{c}(t), \boldsymbol{\theta})d\mathbf{W}(t),\]
where \(\mathbf{f}(\cdot)\) is the drift (deterministic) part, \(\mathbf{g}(\cdot)\) encodes the diffusion (noise) part, \(\boldsymbol{\theta}\) is a parameter set controlling fire spread rates, fuel consumption rates, etc., \(\mathbf{c}(t)\) represents control inputs (e.g., firefighting efforts, fuel breaks), and \(\mathbf{W}(t)\) is a vector Wiener process (Brownian motion).
Inspired by the Lotka-Volterra dynamics in the cancer model, we can define the system as
\[\begin{aligned} dFI &= \left[\alpha FI \left(1 - \frac{FI}{K_{FI}}\right) + \beta FI \cdot FC - \gamma_{\text{wind}} w(t) FI - \phi_{\text{control}} c(t) \right] dt + \sigma_{FI} FI \, dW_1\\ dFC &= \left[\delta FI \cdot FC - \epsilon FC + \zeta_{\text{dryness}} d(t) \right] dt + \sigma_{FC} \sqrt{FC} \, dW_2\\ dE_i &= \left[\mu_i(E_{i,\text{eq}} - E_i) + \eta_i FI - \nu_i FC \right] dt + \sigma_{E_i} \, dW_{i+2} \end{aligned},\]
where
\(\alpha\): Fire growth rate
\(K_{FI}\): Carrying capacity (maximum intensity)
\(\beta\): Rate at which fuel consumption increases fire intensity
\(\gamma_{\text{wind}}\): Effect of wind on fire spread
\(\phi_{\text{control}}\): Effectiveness of firefighting control efforts
\(\delta\): Rate at which fire intensity increases fuel consumption
\(\epsilon\): Natural fuel depletion rate
\(\zeta_{\text{dryness}}\): Effect of dryness on fuel consumption
\(\mu_i, \eta_i, \nu_i\): Environmental factor parameters
\(\sigma\): Noise intensities for each variable.
Our implementation is based on R.
3 Basic 2D SDE Model Implementation
First, let’s implement a simplified 2D model focusing just on Fire Intensity (FI) and Fuel Consumption (FC).
library(Sim.DiffProc)
library(plotly)
# Parameters
alphaP <- 0.4 # Fire growth rate
K_FI <- 2.0 # Maximum fire intensity
betaP <- 0.03 # Fuel-to-intensity conversion rate
deltaP <- 0.04 # Intensity-to-fuel-consumption rate
epsilonP <- 0.3 # Fuel depletion rate
gamma_wind <- 0.05 # Wind effect on fire spread
zeta_dry <- 0.02 # Dryness effect on fuel consumption
sigma_FI <- 0.08 # Noise for fire intensity
sigma_FC <- 0.06 # Noise for fuel consumption
# Wind and dryness functions (time-dependent)
w_t <- function(t) {
# Diurnal wind pattern with random gusts
base <- 0.5 + 0.3 * sin(2*pi*t/24) # 24-hour cycle
gusts <- 0.2 * (runif(1) > 0.8) # Occasional gusts
return(base + gusts)
}
d_t <- function(t) {
# Slowly decreasing humidity (increasing dryness)
return(0.3 + 0.1 * (1 - exp(-t/72)))
}
# Control efforts (firefighting)
c_t <- function(t) {
# Firefighting starts after detection (t > 8) and increases gradually
ifelse(t < 8, 0, 0.3 * (1 - exp(-(t-8)/24)))
}
# Define drift and diffusion for snssde2d()
fx <- expression(
alphaP * x * (1 - x/K_FI) + betaP * x * y - gamma_wind * 0.5 * x -
phi_control * ifelse(t < 8, 0, 0.3 * (1 - exp(-(t-8)/24))),
deltaP * x * y - epsilonP * y + zeta_dry * (0.3 + 0.1 * (1 - exp(-t/72)))
)
gx <- expression(
sigma_FI * x,
sigma_FC * sqrt(y)
)
# Set parameters not defined in expressions above
phi_control <- 0.2
# Simulate
set.seed(123)
mod <- snssde2d(
drift = fx,
diffusion = gx,
x0 = c(x = 0.1, y = 0.8), # Initial small fire, high fuel availability
t0 = 0,
T = 120, # 120 hours (5 days) simulation
N = 500,
M = 50 # 50 stochastic trajectories
)
# Visualize
fig <- plot_ly()
for (i in 1:5) { # Plot first 5 paths
fig <- fig %>%
add_trace(x = mod$time, y = mod$X[, i], name = paste("FI", i),
type = "scatter", mode = "lines", line = list(color = "orange")) %>%
add_trace(x = mod$time, y = mod$Y[, i], name = paste("FC", i),
line = list(color = "brown"), yaxis = "y2")
}
fig <- fig %>% layout(
title = "Fire Intensity (FI) vs. Fuel Consumption (FC) - 5 SDE Model Solution Paths across Time",
yaxis = list(title = "Fire Intensity",
tickfont = list(color = "orange"),
titlefont = list(color = "orange")),
yaxis2 = list(overlaying = "y", side = "right", title = "Fuel Consumption",
tickfont = list(color = "brown"),
titlefont = list(color = "brown")),
xaxis = list(title = "Time (Hours)"),
legend = list(x = 1.1, y = 0.9)
)
fig3.1 Parameter Estimation Using Bayesian Optimization
Now, let’s implement parameter estimation using Bayesian Optimization, assuming we have observed data.
library(rBayesianOptimization)
# Create synthetic observed data (as ground truth)
# In a real-world scenario, this would come from actual fire measurements
observed_data <- data.frame(
time = seq(from=mod$t0, to=mod$T, by=mod$Dt),
FI = rowMeans(mod$X),
FC = rowMeans(mod$Y)
)
# Objective function to minimize MSE
objective_func <- function(alphaP, betaP, deltaP) {
sim_mod <- snssde2d(
drift = expression(
alphaP * x * (1 - x/2.0) + betaP * x * y - 0.05 * 0.5 * x - 0.2 *
ifelse(t < 8, 0, 0.3 * (1 - exp(-(t-8)/24))),
deltaP * x * y - 0.3 * y + 0.02 * (0.3 + 0.1 * (1 - exp(-t/72)))
),
diffusion = expression(0.08*x, 0.06*sqrt(y)),
x0 = c(0.1, 0.8),
t0 = 0,
T = 120,
N = 500
)
# Calculate mean squared error between simulation and observed data
mse <- mean((sim_mod$X - observed_data$FI)^2 + (sim_mod$Y - observed_data$FC)^2)
return(list(Score = -mse)) # Negative because we want to maximize Score
}
# Bayesian optimization
opt_result <- BayesianOptimization(
objective_func,
bounds = list(
alphaP = c(0.2, 0.6), # Search range for fire growth rate
betaP = c(0.01, 0.05), # Search range for fuel effect
deltaP = c(0.02, 0.06) # Search range for intensity effect
),
init_points = 5,
n_iter = 20
)## elapsed = 0.00 Round = 1 alphaP = 0.440964 betaP = 0.04482155 deltaP = 0.05951986 Value = -0.01953817
## elapsed = 0.01 Round = 2 alphaP = 0.2091407 betaP = 0.01837854 deltaP = 0.05162415 Value = -0.03995364
## elapsed = 0.02 Round = 3 alphaP = 0.528225 betaP = 0.03315608 deltaP = 0.02034668 Value = -0.02953452
## elapsed = 0.00 Round = 4 alphaP = 0.2146278 betaP = 0.02775184 deltaP = 0.04330916 Value = -0.02437054
## elapsed = 0.02 Round = 5 alphaP = 0.2940198 betaP = 0.0246884 deltaP = 0.03299292 Value = -0.02191498
## elapsed = 0.00 Round = 6 alphaP = 0.2170595 betaP = 0.04943426 deltaP = 0.04007898 Value = -0.02723391
## elapsed = 0.00 Round = 7 alphaP = 0.3334215 betaP = 0.041091 deltaP = 0.05868681 Value = -0.02506321
## elapsed = 0.00 Round = 8 alphaP = 0.598943 betaP = 0.0100 deltaP = 0.0200 Value = -0.05390664
## elapsed = 0.02 Round = 9 alphaP = 0.6000 betaP = 0.04430274 deltaP = 0.0200 Value = -0.03204524
## elapsed = 0.01 Round = 10 alphaP = 0.5710333 betaP = 0.02681722 deltaP = 0.0200 Value = -0.02980931
## elapsed = 0.00 Round = 11 alphaP = 0.6000 betaP = 0.0500 deltaP = 0.04867891 Value = -0.02515796
## elapsed = 0.01 Round = 12 alphaP = 0.347903 betaP = 0.0429017 deltaP = 0.04318399 Value = -0.02221464
## elapsed = 0.01 Round = 13 alphaP = 0.404076 betaP = 0.0500 deltaP = 0.03226481 Value = -0.02832957
## elapsed = 0.00 Round = 14 alphaP = 0.2291243 betaP = 0.0500 deltaP = 0.0600 Value = -0.03983545
## elapsed = 0.00 Round = 15 alphaP = 0.6000 betaP = 0.03454424 deltaP = 0.05952342 Value = -0.03184339
## elapsed = 0.00 Round = 16 alphaP = 0.2000 betaP = 0.03503889 deltaP = 0.02491975 Value = -0.02586037
## elapsed = 0.00 Round = 17 alphaP = 0.4038183 betaP = 0.03074356 deltaP = 0.04173319 Value = -0.04008501
## elapsed = 0.00 Round = 18 alphaP = 0.6000 betaP = 0.04635544 deltaP = 0.04442064 Value = -0.02031306
## elapsed = 0.02 Round = 19 alphaP = 0.2532777 betaP = 0.02497065 deltaP = 0.0560211 Value = -0.02720903
## elapsed = 0.00 Round = 20 alphaP = 0.6000 betaP = 0.03843871 deltaP = 0.03253093 Value = -0.02823137
## elapsed = 0.01 Round = 21 alphaP = 0.5928061 betaP = 0.02093379 deltaP = 0.0200 Value = -0.02566581
## elapsed = 0.02 Round = 22 alphaP = 0.2000 betaP = 0.04504585 deltaP = 0.05006456 Value = -0.02574721
## elapsed = 0.00 Round = 23 alphaP = 0.6000 betaP = 0.04237101 deltaP = 0.0600 Value = -0.02503264
## elapsed = 0.00 Round = 24 alphaP = 0.2000 betaP = 0.0272826 deltaP = 0.0200 Value = -0.03644013
## elapsed = 0.02 Round = 25 alphaP = 0.6000 betaP = 0.02422915 deltaP = 0.04010826 Value = -0.03504895
##
## Best Parameters Found:
## Round = 1 alphaP = 0.440964 betaP = 0.04482155 deltaP = 0.05951986 Value = -0.01953817## alphaP betaP deltaP
## 0.44096395 0.04482155 0.05951986# Set the parameters to the estimates from Bayesian Optimization
alphaP <- opt_result$Best_Par[[1]] # Fire growth rate
betaP <- opt_result$Best_Par[[2]] # Fuel effect
deltaP <- opt_result$Best_Par[[3]] # Intensity effect
# Define drift and diffusion with optimized parameters
fx <- expression(
alphaP * x * (1 - x/K_FI) + betaP * x * y - gamma_wind * 0.5 * x -
phi_control * ifelse(t < 8, 0, 0.3 * (1 - exp(-(t-8)/24))),
deltaP * x * y - epsilonP * y + zeta_dry * (0.3 + 0.1 * (1 - exp(-t/72)))
)
gx <- expression(
sigma_FI * x,
sigma_FC * sqrt(y)
)
# Simulate with optimized parameters
set.seed(123)
mod_opt <- snssde2d(
drift = fx,
diffusion = gx,
x0 = c(x = 0.1, y = 0.8),
t0 = 0,
T = 120,
N = 500,
M = 50
)
# Visualize optimized model
fig_opt <- plot_ly()
for (i in 6:10) { # Plot a different set of paths
fig_opt <- fig_opt %>%
add_trace(x = mod_opt$time, y = mod_opt$X[, i], name = paste("FI", i),
type = "scatter", mode = "lines", line = list(color = "red")) %>%
add_trace(x = mod_opt$time, y = mod_opt$Y[, i], name = paste("FC", i),
line = list(color = "darkgreen"), yaxis = "y2")
}
# Add observed data
fig_opt <- fig_opt %>%
add_trace(x = observed_data$time, y = observed_data$FI, name = "Observed FI",
type = "scatter", mode = "lines", line = list(color = "black", dash = "dash", width = 3)) %>%
add_trace(x = observed_data$time, y = observed_data$FC, name = "Observed FC",
line = list(color = "black", dash = "dot", width = 3), yaxis = "y2")
fig_opt <- fig_opt %>% layout(
title = "Optimized Fire Intensity (FI) vs. Fuel Consumption (FC)",
yaxis = list(title = "Fire Intensity",
tickfont = list(color = "red"),
titlefont = list(color = "red")),
yaxis2 = list(overlaying = "y", side = "right", title = "Fuel Consumption",
tickfont = list(color = "darkgreen"),
titlefont = list(color = "darkgreen")),
xaxis = list(title = "Time (Hours)"),
legend = list(x = 1.1, y = 0.9)
)
fig_opt3.2 Spatial Extension for Wildfire Propagation
To better represent wildfire dynamics, we need to incorporate spatial propagation. Let’s create a lattice-based model where fire spreads across cells.
library(dplyr)
library(tidyr)
library(heatmaply)
# Create a lattice grid for spatial fire simulation
create_fire_grid <- function(rows = 20, cols = 20) {
grid <- expand.grid(x = 1:rows, y = 1:cols) %>%
mutate(
FI = 0, # Fire intensity
FC = runif(n(), 0.6, 1.0), # Fuel content (variable)
elevation = 0, # Flat terrain to start
wind_x = 0.2, # Constant rightward wind initially
wind_y = 0.1 # Slight upward wind component
)
# Add some topographic features (hills)
for(i in 1:3) {
# Random hill centers
center_x <- runif(1, 5, rows-5)
center_y <- runif(1, 5, cols-5)
height <- runif(1, 0.3, 0.7)
width <- runif(1, 3, 6)
# Add hill elevations
grid <- grid %>%
mutate(
dist_sq = (x - center_x)^2 + (y - center_y)^2,
elevation = elevation + 100 * height * exp(-dist_sq / (2 * width^2))
)
}
# Add a fuel break (road or river)
grid <- grid %>%
mutate(
FC = ifelse(abs(y - cols/2) < 2.5, 0.1, FC) # Horizontal break
)
# Normalize elevation for display
grid <- grid %>%
mutate(elevation = elevation / max(elevation))
# Start fire in bottom left corner
ignition_points <- grid %>%
filter(x <= 3, y <= 3)
grid$FI[grid$x <= 3 & grid$y <= 3] <- 0.5
return(grid)
}
# Initialize grid
rows <- 30
cols <- 30
grid <- create_fire_grid(rows, cols)
# Parameters for the spatial fire model
alpha_spread <- 0.15 # Base fire spread rate
beta_wind <- 0.2 # Wind influence factor
gamma_slope <- 0.25 # Slope influence factor
delta_fuel <- 0.3 # Fuel content influence
epsilon_burn <- 0.12 # Rate of fuel consumption
sigma_FI_spatial <- 0.04 # Stochastic component
dt <- 0.5 # Time step
# Function to simulate one time step
simulate_fire_step <- function(grid, t) {
# Create a copy for the updated values
new_grid <- grid
# Update each cell
for(i in 1:nrow(grid)) {
row <- grid[i,]
x <- row$x
y <- row$y
if(row$FI > 0 || any(grid$FI[grid$x >= x-1 & grid$x <= x+1 &
grid$y >= y-1 & grid$y <= y+1] > 0)) {
# Calculate neighboring fire influence
neighbors <- grid %>%
filter(x >= row$x-1, x <= row$x+1, y >= row$y-1, y <= row$y+1)
neighbor_influence <- sum(neighbors$FI) / 9 # Average including self
# Calculate slope effect (fire spreads faster uphill)
if(x < rows && y < cols) {
slope_x <- grid$elevation[grid$x == x+1 & grid$y == y] - row$elevation
slope_y <- grid$elevation[grid$x == x & grid$y == y+1] - row$elevation
} else {
slope_x <- 0
slope_y <- 0
}
slope_effect <- gamma_slope * (max(0, slope_x) + max(0, slope_y))
# Wind effect
wind_effect <- beta_wind * (row$wind_x + row$wind_y)
# Calculate deterministic component of fire intensity change
dFI_deterministic <- (
alpha_spread * neighbor_influence * (1 - row$FI) + # Spread from neighbors
delta_fuel * row$FC * row$FI + # Fuel effect
slope_effect + # Slope effect
wind_effect # Wind effect
) * dt
# Stochastic component
dFI_stochastic <- sigma_FI_spatial * row$FI * sqrt(dt) * rnorm(1)
# Update fire intensity
new_FI <- row$FI + dFI_deterministic + dFI_stochastic
new_FI <- max(0, min(1, new_FI)) # Bound between 0 and 1
# Update fuel consumption
dFC <- -epsilon_burn * row$FI * row$FC * dt
new_FC <- row$FC + dFC
new_FC <- max(0, new_FC) # Cannot be negative
# Store updated values
new_grid$FI[i] <- new_FI
new_grid$FC[i] <- new_FC
}
}
# Firefighting efforts activated after time t > 15
if(t > 15) {
# Concentrate efforts on the fire front (cells with medium intensity)
firefighting_cells <- which(new_grid$FI > 0.3 & new_grid$FI < 0.7)
if(length(firefighting_cells) > 0) {
# Apply suppression to the top 20% most critical cells
n_critical <- max(1, floor(0.2 * length(firefighting_cells)))
critical_cells <- sample(firefighting_cells, n_critical)
# Reduce fire intensity in these cells
new_grid$FI[critical_cells] <- new_grid$FI[critical_cells] * 0.8
}
}
return(new_grid)
}
# Run simulation for multiple time steps
time_steps <- 40
fire_history <- list()
fire_history[[1]] <- grid
for(t in 2:time_steps) {
cat("Simulating time step", t, "of", time_steps, "\n")
fire_history[[t]] <- simulate_fire_step(fire_history[[t-1]], t)
}## Simulating time step 2 of 40
## Simulating time step 3 of 40
## Simulating time step 4 of 40
## Simulating time step 5 of 40
## Simulating time step 6 of 40
## Simulating time step 7 of 40
## Simulating time step 8 of 40
## Simulating time step 9 of 40
## Simulating time step 10 of 40
## Simulating time step 11 of 40
## Simulating time step 12 of 40
## Simulating time step 13 of 40
## Simulating time step 14 of 40
## Simulating time step 15 of 40
## Simulating time step 16 of 40
## Simulating time step 17 of 40
## Simulating time step 18 of 40
## Simulating time step 19 of 40
## Simulating time step 20 of 40
## Simulating time step 21 of 40
## Simulating time step 22 of 40
## Simulating time step 23 of 40
## Simulating time step 24 of 40
## Simulating time step 25 of 40
## Simulating time step 26 of 40
## Simulating time step 27 of 40
## Simulating time step 28 of 40
## Simulating time step 29 of 40
## Simulating time step 30 of 40
## Simulating time step 31 of 40
## Simulating time step 32 of 40
## Simulating time step 33 of 40
## Simulating time step 34 of 40
## Simulating time step 35 of 40
## Simulating time step 36 of 40
## Simulating time step 37 of 40
## Simulating time step 38 of 40
## Simulating time step 39 of 40
## Simulating time step 40 of 40Now, let’s visualize the fire propagation over time:
# Function to plot a single time step as a heatmap
plot_fire_snapshot <- function(grid_data, title) {
# Reshape for heatmap
fire_matrix <- grid_data %>%
select(x, y, FI) %>%
pivot_wider(names_from = y, values_from = FI) %>%
select(-x) %>%
as.matrix()
# Create heatmap
heatmaply(fire_matrix,
colors = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
showticklabels = c(FALSE, FALSE),
main = title,
xlab = "",
ylab = "",
margins = c(0, 0, 0, 0))
}
# Create interactive time series of fire propagation
# Select a subset of times to display
selected_times <- seq(1, time_steps, by = 4)
# Plot individual snapshots
for(t in selected_times) {
plot_title <- paste("Fire Propagation at t =", t)
print(plot_fire_snapshot(fire_history[[t]], plot_title))
}
# Create animation frames for a comprehensive view
library(animation)
# Set up animation
ani.options(interval = 0.3)
# Function to plot a single frame
plot_fire_frame <- function(t) {
grid_data <- fire_history[[t]]
# Create separate matrices for FI, FC and elevation
FI_matrix <- grid_data %>%
select(x, y, FI) %>%
pivot_wider(names_from = y, values_from = FI) %>%
select(-x) %>%
as.matrix()
FC_matrix <- grid_data %>%
select(x, y, FC) %>%
pivot_wider(names_from = y, values_from = FC) %>%
select(-x) %>%
as.matrix()
elevation_matrix <- grid_data %>%
select(x, y, elevation) %>%
pivot_wider(names_from = y, values_from = elevation) %>%
select(-x) %>%
as.matrix()
# Set up plot layout (2x2 grid)
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))
# Plot fire intensity
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = paste("Fire Intensity at t =", t),
xlab = "", ylab = "")
# Plot fuel consumption
image(1:nrow(FC_matrix), 1:ncol(FC_matrix), FC_matrix,
col = colorRampPalette(c("white", "lightgreen", "darkgreen"))(100),
main = "Remaining Fuel",
xlab = "", ylab = "")
# Plot elevation
image(1:nrow(elevation_matrix), 1:ncol(elevation_matrix), elevation_matrix,
col = terrain.colors(100),
main = "Elevation",
xlab = "", ylab = "")
# Plot combined view (elevation with fire overlay)
image(1:nrow(elevation_matrix), 1:ncol(elevation_matrix), elevation_matrix,
col = terrain.colors(100),
main = "Combined View",
xlab = "", ylab = "")
# Overlay fire with transparency
# Add fire intensities where they exceed a threshold
fire_cells <- which(FI_matrix > 0.1, arr.ind = TRUE)
if(nrow(fire_cells) > 0) {
points(fire_cells[, 1], fire_cells[, 2],
pch = 15, cex = 0.5 + 1.5 * FI_matrix[fire_cells],
col = adjustcolor("red", alpha.f = 0.7))
}
}
# Save animation with corrected fps parameter
# Save animation with corrected parameters
saveGIF({
for(t in 1:time_steps) {
plot_fire_frame(t)
}
}, movie.name = "wildfire_propagation.gif", ani.width = 800, ani.height = 800,
interval = 0.25) # Use a value that results in fps being a factor of 100## [1] TRUE
# Advanced High-Dimensional
(nD) SDE Model
For a more realistic model, we’ll expand to a higher-dimensional SDE
system including environmental variables. This would typically be
implemented in Julia due to performance considerations,
similar to the cancer model.
# This is a conceptual outline, implementation would typically use Julia
# State variables:
# X(t) = [FI(t), FC(t), Wind_Speed(t), Wind_Direction(t), Temperature(t), Humidity(t), Precipitation(t)]
# SDE system:
# dFI = [...] dt + [...] dW_1
# dFC = [...] dt + [...] dW_2
# dWind_Speed = [...] dt + [...] dW_3
# dWind_Direction = [...] dt + [...] dW_4
# dTemperature = [...] dt + [...] dW_5
# dHumidity = [...] dt + [...] dW_6
# dPrecipitation = [...] dt + [...] dW_7
# The Julia implementation would:
# 1. Define the SDE problem
# 2. Use DifferentialEquations.jl for efficient numerical solutions
# 3. Parallelize over multiple parameter sets for optimization
# 4. Leverage GPU acceleration for spatial simulations
# 5. Use Bayesian optimization for parameter estimation3.3 Diagnostic Analytics and Model Validation
# Analysis of simulation results
# 1. Extract overall fire progression metrics
fire_progression <- data.frame(
Time = 1:time_steps,
TotalIntensity = sapply(fire_history, function(grid) sum(grid$FI)),
BurnedArea = sapply(fire_history, function(grid) sum(grid$FI > 0.1)),
AvgIntensity = sapply(fire_history, function(grid) mean(grid$FI[grid$FI > 0]))
)
# Visualize fire progression metrics
progression_plot <- plot_ly(fire_progression, x = ~Time) %>%
add_trace(y = ~TotalIntensity, name = "Total Intensity",
type = "scatter", mode = "lines", line = list(color = "red")) %>%
add_trace(y = ~BurnedArea, name = "Burned Area",
yaxis = "y2", line = list(color = "orange")) %>%
add_trace(y = ~AvgIntensity * 100, name = "Avg Intensity",
yaxis = "y3", line = list(color = "darkred")) %>%
layout(
title = "Fire Progression Metrics",
xaxis = list(title = "Time Steps"),
yaxis = list(title = "Total Intensity",
titlefont = list(color = "red"),
tickfont = list(color = "red")),
yaxis2 = list(title = "Burned Area (cells)",
overlaying = "y", side = "right",
titlefont = list(color = "orange"),
tickfont = list(color = "orange")),
yaxis3 = list(title = "Avg Intensity (%)",
overlaying = "y", side = "right", position = 0.85,
titlefont = list(color = "darkred"),
tickfont = list(color = "darkred")),
legend = list(x = 0.1, y = 0.9)
)
progression_plot# 2. Fire spread direction analysis
analyze_fire_direction <- function(grid_history) {
# Track fire centroid over time
centroids <- data.frame(
Time = 1:length(grid_history),
X_centroid = sapply(grid_history, function(grid) {
if(sum(grid$FI) > 0) {
return(weighted.mean(grid$x, grid$FI))
} else {
return(NA)
}
}),
Y_centroid = sapply(grid_history, function(grid) {
if(sum(grid$FI) > 0) {
return(weighted.mean(grid$y, grid$FI))
} else {
return(NA)
}
})
)
# Calculate direction of movement
centroids$Direction <- NA
for(i in 2:nrow(centroids)) {
if(!is.na(centroids$X_centroid[i-1]) && !is.na(centroids$X_centroid[i])) {
dx <- centroids$X_centroid[i] - centroids$X_centroid[i-1]
dy <- centroids$Y_centroid[i] - centroids$Y_centroid[i-1]
centroids$Direction[i] <- atan2(dy, dx) * 180 / pi # in degrees
centroids$Speed[i] <- sqrt(dx^2 + dy^2)
}
}
return(centroids)
}
fire_direction <- analyze_fire_direction(fire_history)
# Plot fire centroid movement
centroid_plot <- plot_ly(fire_direction) %>%
add_trace(x = ~X_centroid, y = ~Y_centroid, color = ~Time,
type = "scatter", mode = "markers+lines",
marker = list(size = 10),
line = list(width = 2),
text = ~paste("Time:", Time, "<br>Direction:", round(Direction, 1), "°"),
hoverinfo = "text") %>%
layout(
title = "Fire Centroid Movement",
xaxis = list(title = "X coordinate", range = c(1, rows)),
yaxis = list(title = "Y coordinate", range = c(1, cols)),
coloraxis = list(colorscale = "Inferno")
)
centroid_plot3.4 Multi-Core Parallel Processing and Model Optimization
For high-dimensional SDE models with spatial components, parallel processing becomes essential. Below is a conceptual outline of how this would be implemented.
# Outline for parallel processing (would be implemented in Julia)
# 1. Define parameter space for optimization
# parameters = [alpha_spread, beta_wind, gamma_slope, delta_fuel, epsilon_burn, ...]
# 2. Use parallel batch processing for Bayesian optimization
# - Each worker evaluates different parameter sets
# - Multiple stochastic trajectories per parameter set
# - Cost function based on fitting to observed fire data
# 3. Implement spatial parallelism
# - Divide spatial grid into sub-regions
# - Process each sub-region
# 4. Implement parallel ensemble simulations
# For each parameter set, run multiple stochastic realizations in parallel
parallel_ensemble_simulation <- function(params, n_ensembles = 20) {
# This would be implemented with future or foreach in R
# or with pmap in Julia
library(parallel)
library(doParallel)
# Set up parallel backend
cores <- detectCores() - 1
cl <- makeCluster(cores)
registerDoParallel(cl)
# Run ensemble of simulations in parallel
results <- foreach(i = 1:n_ensembles, .combine = 'rbind') %dopar% {
# Set random seed based on iteration
set.seed(123 + i)
# Run one full simulation
grid <- create_fire_grid(rows, cols)
history <- list(grid)
for(t in 2:time_steps) {
history[[t]] <- simulate_fire_step(history[[t-1]], t, params)
}
# Extract metrics
data.frame(
ensemble = i,
total_burned = sum(history[[time_steps]]$FI > 0.1),
max_intensity = max(sapply(history, function(g) max(g$FI))),
burn_duration = max(which(sapply(history, function(g) sum(g$FI)) > 0))
)
}
# Clean up
stopCluster(cl)
return(results)
}
# 5. Define cost function for parameter optimization
sde_cost_function <- function(params, observed_data) {
# Run ensemble
ensemble_results <- parallel_ensemble_simulation(params)
# Compare with observed data
# For example, compare burn patterns, fire intensity over time, etc.
# Return negative MSE or other cost metric
# This is a placeholder - actual implementation would depend on observed data format
mean_metrics <- colMeans(ensemble_results[, c("total_burned", "max_intensity", "burn_duration")])
# Calculate cost (e.g., weighted MSE against observed values)
weights <- c(0.4, 0.3, 0.3) # Relative importance of each metric
observed_values <- c(observed_data$total_burned, observed_data$max_intensity, observed_data$burn_duration)
cost <- sum(weights * ((mean_metrics - observed_values) / observed_values)^2)
return(cost)
}
# 6. Run Bayesian Optimization
bayesian_opt_fire_model <- function(observed_data) {
library(rBayesianOptimization)
# Wrapper function for BO
obj_func <- function(alpha_spread, beta_wind, gamma_slope, delta_fuel, epsilon_burn) {
params <- list(
alpha_spread = alpha_spread,
beta_wind = beta_wind,
gamma_slope = gamma_slope,
delta_fuel = delta_fuel,
epsilon_burn = epsilon_burn
)
cost <- sde_cost_function(params, observed_data)
return(list(Score = -cost)) # Negative because we want to maximize Score
}
# Run optimization
opt_result <- BayesianOptimization(
obj_func,
bounds = list(
alpha_spread = c(0.05, 0.25),
beta_wind = c(0.1, 0.4),
gamma_slope = c(0.1, 0.5),
delta_fuel = c(0.1, 0.5),
epsilon_burn = c(0.05, 0.2)
),
init_points = 10,
n_iter = 30
)
return(opt_result)
}3.5 Spatio-Temporal Visualization and Advanced Analytics
Let’s now implement more sophisticated visualizations and analytics for our wildfire model.
# 1. Create a 3D visualization of fire intensity over landscape
library(plotly)
# Function to create 3D surface plot at a specific time
plot_fire_3d <- function(grid_data, t=30) {
# Convert to matrices
elevation_matrix <- grid_data %>%
select(x, y, elevation) %>%
pivot_wider(names_from = y, values_from = elevation) %>%
select(-x) %>%
as.matrix()
fire_matrix <- grid_data %>%
select(x, y, FI) %>%
pivot_wider(names_from = y, values_from = FI) %>%
select(-x) %>%
as.matrix()
# Create scaling factor for fire visualization
z_scale <- 3 # Height multiplier for fire
# Create 3D surface plot
fig <- plot_ly() %>%
add_surface(
z = elevation_matrix,
colorscale = list(c(0, 1), c("tan", "darkgreen")),
opacity = 0.8,
showscale = FALSE, name ="Elevation"
)
# Add fire as heightened surface
fire_height <- elevation_matrix + fire_matrix * z_scale
# Only show fire where intensity > 0.1
for(i in 1:nrow(fire_matrix)) {
for(j in 1:ncol(fire_matrix)) {
if(fire_matrix[i,j] <= 0.1) {
fire_height[i,j] <- NA # Make non-burning areas transparent
}
}
}
fig <- fig %>%
add_surface(
z = fire_height,
colorscale = list(
c(0, 0.3, 0.6, 0.8, 1),
c("yellow", "orange", "red", "darkred", "black")
),
opacity = 0.9,
showscale = FALSE
)
# Add layout settings
fig <- fig %>% layout(
title = "3D Wildfire Visualization (Time=30 out of 40)",
scene = list(
aspectratio = list(x = 1, y = 1, z = 0.5),
xaxis = list(title = ""),
yaxis = list(title = ""),
zaxis = list(title = "")
)
)
return(fig)
}
# Create 3D visualization for a specific time point t=10
plot_fire_3d(fire_history[[30]], t=30)# 2. Create a risk assessment map based on simulation ensembles
create_risk_map <- function(ensemble_histories, threshold = 0.5) {
# Assuming ensemble_histories is a list of lists (outer: ensemble, inner: time steps)
# Create an empty risk map
risk_map <- create_fire_grid(rows, cols)
risk_map$risk <- 0
# Count how many ensembles burned each cell above threshold
n_ensembles <- length(ensemble_histories)
for(e in 1:n_ensembles) {
# Get final state of this ensemble
final_state <- ensemble_histories[[e]][[time_steps]]
# Mark cells that burned above threshold
for(i in 1:nrow(risk_map)) {
cell_x <- risk_map$x[i]
cell_y <- risk_map$y[i]
# Find corresponding cell in final state
final_FI <- final_state$FI[final_state$x == cell_x & final_state$y == cell_y]
if(length(final_FI) > 0 && final_FI > threshold) {
risk_map$risk[i] <- risk_map$risk[i] + 1
}
}
}
# Convert to probability
risk_map$risk <- risk_map$risk / n_ensembles
return(risk_map)
}
# For demonstration, create a synthetic ensemble history
# This would normally come from actual parallel runs
synthetic_ensemble_histories <- list()
for(e in 1:10) {
# Perturb the original history slightly for each ensemble
ensemble_history <- lapply(fire_history, function(grid) {
grid$FI <- pmax(0, pmin(1, grid$FI + rnorm(nrow(grid), 0, 0.05)))
return(grid)
})
synthetic_ensemble_histories[[e]] <- ensemble_history
}
# Create risk map
risk_map <- create_risk_map(synthetic_ensemble_histories, threshold = 0.3)
# Visualize risk map
risk_matrix <- risk_map %>%
select(x, y, risk) %>%
pivot_wider(names_from = y, values_from = risk) %>%
select(-x) %>%
as.matrix()
risk_plot <- plot_ly(
z = risk_matrix,
type = "heatmap",
colorscale = "Reds"
) %>%
layout(
title = "Wildfire Risk Map (Probability of Burning)",
xaxis = list(title = ""),
yaxis = list(title = "")
)
risk_plot# 3. Fire spread dynamics analysis
analyze_spread_rate <- function(grid_history) {
# Calculate burned area at each time step
burned_area <- sapply(grid_history, function(grid) sum(grid$FI > 0.1))
# Calculate spread rate (derivative of burned area)
spread_rate <- c(0, diff(burned_area))
# Create data frame
spread_data <- data.frame(
Time = 1:length(grid_history),
BurnedArea = burned_area,
SpreadRate = spread_rate
)
return(spread_data)
}
spread_data <- analyze_spread_rate(fire_history)
# Plot spread rate
spread_plot <- plot_ly(spread_data, x = ~Time) %>%
add_trace(y = ~BurnedArea, name = "Burned Area",
type = "scatter", mode = "lines", line = list(color = "orange")) %>%
add_trace(y = ~SpreadRate, name = "Spread Rate",
yaxis = "y2", line = list(color = "red")) %>%
layout(
title = "Fire Spread Dynamics",
xaxis = list(title = "Time Steps"),
yaxis = list(title = "Burned Area (cells)",
titlefont = list(color = "orange"),
tickfont = list(color = "orange")),
yaxis2 = list(title = "Spread Rate (cells/step)",
overlaying = "y", side = "right",
titlefont = list(color = "red"),
tickfont = list(color = "red")),
legend = list(x = 0.1, y = 0.9)
)
spread_plot3.6 Model Comparison and Validation with Real Data
we can empirically validate the complete SDE wildfire modeling system, provided we have real wildfire data available.
# This would typically use real wildfire data
# Here we outline the approach
# 1. Prepare validation data
prepare_validation_data <- function(file_path) {
# Load real fire perimeter data
# This could be from satellite imagery, fire agency reports, etc.
validation_data <- read.csv(file_path)
# Preprocess to match our model grid
# - Convert coordinates to grid cells
# - Extract burn progression over time
# - Calculate relevant metrics
return(validation_data)
}
# 2. Validation metrics
calculate_validation_metrics <- function(simulated_history, validation_data) {
# Calculate various metrics:
# - Total burned area accuracy
# - Fire perimeter match (e.g., Hausdorff distance)
# - Rate of spread correlation
# - Burn severity comparison
# Example metrics (placeholder)
metrics <- list(
area_accuracy = 0.85, # Ratio of correct/total area
perimeter_f1 = 0.78, # F1 score for perimeter prediction
spread_rate_corr = 0.72, # Correlation in spread rates
rmse = 0.15 # Root mean square error
)
return(metrics)
}
# 3. Model comparison function
compare_models <- function(model_params_list, validation_data) {
# Run each model parameterization
model_results <- list()
metrics_table <- data.frame()
for (i in 1:length(model_params_list)) {
params <- model_params_list[[i]]
# Run model with these parameters
history <- run_fire_simulation(params)
# Calculate validation metrics
metrics <- calculate_validation_metrics(history, validation_data)
# Store results
model_results[[i]] <- history
metrics_table <- rbind(metrics_table,
data.frame(model = i,
params = I(list(params)),
metrics = I(list(metrics))))
}
return(list(
results = model_results,
metrics = metrics_table
))
}Below is the complete end-to-end wildfire SDE model implementation.
# 1. Configure model parameters
model_config <- list(
# Grid parameters
grid_size = list(rows = 50, cols = 50),
cell_size = 100, # meters
# Fire dynamics parameters
alpha_spread = 0.15, # Base fire spread rate
beta_wind = 0.2, # Wind influence factor
gamma_slope = 0.25, # Slope influence factor
delta_fuel = 0.3, # Fuel content influence
epsilon_burn = 0.12, # Rate of fuel consumption
# Stochastic parameters
sigma_FI = 0.04, # Fire intensity noise
sigma_wind = 0.02, # Wind variability
# Simulation parameters
time_steps = 100,
dt = 0.5, # Time step size
# Firefighting parameters
firefighting_start = 20, # When firefighting begins
firefighting_strength = 0.2 # Effectiveness factor
)
# 2. Create the main simulation function
run_wildfire_simulation <- function(config, initial_ignition = NULL, ensemble_size = 1) {
# Initialize ensemble storage
ensemble_results <- list()
# Run each ensemble member
for(e in 1:ensemble_size) {
# Create initial grid
grid <- create_fire_grid(config$grid_size$rows, config$grid_size$cols)
# Set initial ignition if provided
if(!is.null(initial_ignition)) {
ignition_points <- grid %>%
filter(x >= initial_ignition$x_min & x <= initial_ignition$x_max &
y >= initial_ignition$y_min & y <= initial_ignition$y_max)
grid$FI[grid$x >= initial_ignition$x_min & grid$x <= initial_ignition$x_max &
grid$y >= initial_ignition$y_min & grid$y <= initial_ignition$y_max] <-
initial_ignition$intensity
}
# Initialize history
history <- list()
history[[1]] <- grid
# Run simulation
for(t in 2:config$time_steps) {
history[[t]] <- simulate_fire_step(history[[t-1]], t, config)
}
# Store ensemble member
ensemble_results[[e]] <- history
}
# Calculate ensemble statistics
if(ensemble_size > 1) {
ensemble_stats <- calculate_ensemble_statistics(ensemble_results)
return(list(
ensembles = ensemble_results,
stats = ensemble_stats
))
} else {
return(ensemble_results[[1]])
}
}
# 3. Integration with GIS and weather data
integrate_gis_data <- function(grid, gis_file) {
# Load GIS data (elevation, vegetation, etc.)
gis_data <- read_gis_data(gis_file)
# Interpolate to our grid
for(i in 1:nrow(grid)) {
x_coord <- grid$x[i] * grid$cell_size # Convert to real-world coordinates
y_coord <- grid$y[i] * grid$cell_size
# Update grid properties from GIS
grid$elevation[i] <- interpolate_elevation(gis_data, x_coord, y_coord)
grid$FC[i] <- interpolate_fuel(gis_data, x_coord, y_coord)
# etc...
}
return(grid)
}
# 4. Example complete workflow
wildfire_prediction_workflow <- function(area_name, gis_data_path, weather_data_path,
ignition_point, simulation_days,
ensemble_size = 50) {
# Load configuration
config <- get_default_config()
config$time_steps <- simulation_days * 24 # Hourly steps
# Create base grid
grid <- create_fire_grid(config$grid_size$rows, config$grid_size$cols)
# Integrate GIS data
grid <- integrate_gis_data(grid, gis_data_path)
# Set weather forecast
weather <- load_weather_forecast(weather_data_path)
config$weather <- weather
# Set ignition
ignition <- list(
x_min = ignition_point$x - 1,
x_max = ignition_point$x + 1,
y_min = ignition_point$y - 1,
y_max = ignition_point$y + 1,
intensity = 0.5
)
# Run ensemble simulation
results <- run_wildfire_simulation(config, ignition, ensemble_size)
# Generate reports
spread_report <- generate_spread_report(results)
risk_map <- create_risk_map(results$ensembles)
threat_assessment <- assess_infrastructure_threat(results, area_name)
# Return compiled results
return(list(
simulation = results,
spread_report = spread_report,
risk_map = risk_map,
threat_assessment = threat_assessment
))
}This wildfire propagation SDE model captures some of the key physical processes and stochastic elements of wildfire spread. The model includes
Dynamic Coupling: Fire intensity and fuel consumption are dynamically coupled, similar to predator-prey dynamics.
Environmental Factors: The model integrates environmental variables like wind, slope, and fuel characteristics.
Stochastic Elements: Random variations reflect the inherent unpredictability of fire behavior.
Spatial Representation: The model captures spatial heterogeneity in landscape and fire spread.
Parameter Estimation: Bayesian optimization provides a robust framework for tuning model parameters to match observed data.
4 Extended Wildfire SDE Model
The initial SDE model can be enhanced by exploring the following improvements
Finer Spatial Resolution: Implementing a higher-resolution grid for more detailed spatial modeling.
Advanced Weather Integration: Incorporating more sophisticated weather models and forecasts.
Fuel Type Classification: Differentiating between various vegetation and fuel types.
Spotting Behavior: Modeling long-distance ember transport and spot fires.
Fire-Atmosphere Interactions: Coupling with atmospheric models to capture fire-induced weather.
Machine Learning Enhancement: Using ML techniques to refine parameter estimation and improve predictions.
Real-Time Data Assimilation: Incorporating satellite and sensor data for real-time model updating.
Some of these enhancements are implemented below.
4.1 Advanced Weather Integration
Again, the code below assumes we have access to NCDF4 weather forecast data …
# Advanced Weather Integration
library(ncdf4) # For reading weather forecast data
# Function to load and process weather forecast data
process_weather_forecast <- function(forecast_file, grid_coords) {
# Read NetCDF weather forecast (typical format for meteorological data)
nc <- nc_open(forecast_file)
# Extract relevant variables
times <- ncvar_get(nc, "time")
lats <- ncvar_get(nc, "latitude")
lons <- ncvar_get(nc, "longitude")
# Extract weather variables
wind_speed <- ncvar_get(nc, "wind_speed")
wind_direction <- ncvar_get(nc, "wind_direction")
temperature <- ncvar_get(nc, "temperature")
relative_humidity <- ncvar_get(nc, "relative_humidity")
precipitation <- ncvar_get(nc, "precipitation")
# Close NetCDF connection
nc_close(nc)
# Convert to data frame for easier handling
weather_df <- expand.grid(
time_idx = 1:length(times),
lat_idx = 1:length(lats),
lon_idx = 1:length(lons)
)
# Add actual values
weather_df$time <- rep(times, each = length(lats) * length(lons))
weather_df$lat <- rep(rep(lats, each = length(lons)), length(times))
weather_df$lon <- rep(rep(lons, times = length(lats)), length(times))
# Add weather variables
weather_df$wind_speed <- as.vector(wind_speed)
weather_df$wind_direction <- as.vector(wind_direction)
weather_df$temperature <- as.vector(temperature)
weather_df$relative_humidity <- as.vector(relative_humidity)
weather_df$precipitation <- as.vector(precipitation)
# Calculate derived variables
weather_df$wind_x <- weather_df$wind_speed * cos(weather_df$wind_direction * pi/180)
weather_df$wind_y <- weather_df$wind_speed * sin(weather_df$wind_direction * pi/180)
# Calculate fire danger indices
weather_df$ffdi <- calculate_fire_danger_index(
weather_df$temperature,
weather_df$relative_humidity,
weather_df$wind_speed,
weather_df$precipitation
)
return(weather_df)
}
# Function to interpolate weather data to simulation grid at specific time
get_weather_at_time <- function(weather_df, grid, sim_time) {
# Find closest forecast time
time_idx <- which.min(abs(weather_df$time - sim_time))
# Filter weather data for this time
time_weather <- weather_df[weather_df$time_idx == time_idx,]
# For each grid cell, interpolate weather variables
for(i in 1:nrow(grid)) {
# Get cell coordinates
cell_lat <- grid$lat[i]
cell_lon <- grid$lon[i]
# Spatial interpolation (inverse distance weighting)
distances <- sqrt((time_weather$lat - cell_lat)^2 + (time_weather$lon - cell_lon)^2)
weights <- 1/pmax(distances, 0.001)^2
weights <- weights/sum(weights)
# Interpolate variables
grid$wind_x[i] <- sum(time_weather$wind_x * weights)
grid$wind_y[i] <- sum(time_weather$wind_y * weights)
grid$temperature[i] <- sum(time_weather$temperature * weights)
grid$relative_humidity[i] <- sum(time_weather$relative_humidity * weights)
grid$precipitation[i] <- sum(time_weather$precipitation * weights)
grid$ffdi[i] <- sum(time_weather$ffdi * weights)
}
return(grid)
}
# McArthur Fire Danger Index calculation
calculate_fire_danger_index <- function(temp, rh, wind, rain) {
# Drought factor (simplified)
drought <- pmin(10, pmax(0, 10 - rain * 10))
# Calculate FFDI using McArthur's equation
ffdi <- 2 * exp(0.05 * drought) * exp(0.9 * log(max(wind, 3)/3)) *
exp(0.14 * (temp - 20)) * exp(-0.04 * rh)
return(ffdi)
}
# Modify the SDE model to incorporate weather effects
update_sde_with_weather <- function(fx, gx, weather_effects) {
# Enhance drift term with weather effects
fx_weather <- expression(
# Fire intensity equation
alphaP * x * (1 - x/K_FI) + betaP * x * y -
gamma_wind * sqrt(wind_x^2 + wind_y^2) * x * (1 + delta_temp * (temperature - 20)/10) *
(1 - rh_effect * (relative_humidity/100)^2) -
phi_control * ifelse(t < 8, 0, 0.3 * (1 - exp(-(t-8)/24))) -
rain_suppression * precipitation,
# Fuel consumption equation
deltaP * x * y * (1 + temp_effect * (temperature - 20)/10) *
(1 - humidity_effect * (relative_humidity/100)) -
epsilonP * y + zeta_dry * (0.3 + 0.1 * (1 - exp(-t/72))) *
(1 - precipitation_effect * precipitation)
)
# Add stochastic effects for weather uncertainty
gx_weather <- expression(
sigma_FI * x * (1 + 0.2 * ffdi/50), # Higher variability with higher fire danger
sigma_FC * sqrt(y) * (1 + 0.1 * ffdi/50)
)
return(list(fx = fx_weather, gx = gx_weather))
}4.2 Fuel Type Classification
The SDE model below utilizes a more sophisticated fuel model with different vegetation types.
# Fuel Type Classification
# Define fuel types and their properties
define_fuel_types <- function() {
fuel_types <- list(
grass = list(
name = "Grassland",
load = 0.3, # kg/m²
moisture_drying = 0.3, # Drying rate coefficient
moisture_threshold = 15, # % moisture content threshold for burning
heat_content = 18000, # kJ/kg
spread_rate_coef = 0.8, # High spread rate
intensity_coef = 0.5, # Lower intensity
sustainability = 0.4, # Burns quickly
ignition_temp = 320 # Celsius
),
shrub = list(
name = "Shrubland",
load = 0.8,
moisture_drying = 0.2,
moisture_threshold = 25,
heat_content = 20000,
spread_rate_coef = 0.6,
intensity_coef = 0.7,
sustainability = 0.6,
ignition_temp = 340
),
forest = list(
name = "Forest",
load = 2.5,
moisture_drying = 0.1,
moisture_threshold = 30,
heat_content = 22000,
spread_rate_coef = 0.4,
intensity_coef = 0.9,
sustainability = 0.8,
ignition_temp = 360
),
slash = list(
name = "Logging Slash",
load = 3.5,
moisture_drying = 0.15,
moisture_threshold = 25,
heat_content = 21000,
spread_rate_coef = 0.5,
intensity_coef = 1.0,
sustainability = 0.9,
ignition_temp = 330
)
)
return(fuel_types)
}
# Function to assign fuel types based on land cover data
assign_fuel_types <- function(grid, landcover_file) {
# Read land cover data (e.g., from GeoTIFF)
landcover <- raster(landcover_file)
# Extract landcover values at grid cell locations
lc_values <- extract(landcover, cbind(grid$lon, grid$lat))
# Define mapping from land cover codes to fuel types
# This mapping would depend on the specific land cover classification used
lc_to_fuel <- c(
"1" = "grass", # e.g., for land cover code 1 (grassland)
"2" = "shrub", # e.g., for land cover code 2 (shrubland)
"3" = "forest", # e.g., for land cover code 3 (forest)
"4" = "slash" # e.g., for land cover code 4 (recently logged)
)
# Assign fuel type to each grid cell
grid$fuel_type <- lc_to_fuel[as.character(lc_values)]
# Replace NA values with default
grid$fuel_type[is.na(grid$fuel_type)] <- "grass"
# Get fuel type properties
fuel_types <- define_fuel_types()
# Add fuel properties to grid
for(i in 1:nrow(grid)) {
ft <- fuel_types[[grid$fuel_type[i]]]
grid$fuel_load[i] <- ft$load
grid$moisture_drying[i] <- ft$moisture_drying
grid$moisture_threshold[i] <- ft$moisture_threshold
grid$heat_content[i] <- ft$heat_content
grid$spread_rate_coef[i] <- ft$spread_rate_coef
grid$intensity_coef[i] <- ft$intensity_coef
grid$sustainability[i] <- ft$sustainability
grid$ignition_temp[i] <- ft$ignition_temp
}
return(grid)
}
# Update SDE model to incorporate fuel type effects
update_sde_with_fuel_types <- function(fx, gx) {
# Modify drift terms to incorporate fuel properties
fx_fuels <- expression(
# Fire intensity equation - now depends on fuel heat content and load
alphaP * x * (1 - x/K_FI) * spread_rate_coef +
betaP * x * y * intensity_coef * (heat_content/20000) * (fuel_load/1.0) -
gamma_wind * sqrt(wind_x^2 + wind_y^2) * x -
phi_control * ifelse(t < 8, 0, 0.3 * (1 - exp(-(t-8)/24))),
# Fuel consumption equation - now depends on sustainability
deltaP * x * y * sustainability -
epsilonP * y +
zeta_dry * (0.3 + 0.1 * (1 - exp(-t/72))) *
ifelse(relative_humidity < moisture_threshold,
(1 - relative_humidity/moisture_threshold)^2, 0)
)
# Modify diffusion terms
gx_fuels <- expression(
sigma_FI * x * spread_rate_coef,
sigma_FC * sqrt(y) * sustainability
)
return(list(fx = fx_fuels, gx = gx_fuels))
}4.3 Spotting Behavior
An implementation of spotting (long-distance ember transport) is demonstrated below.
# Spotting Behavior Model
# Function to simulate ember generation and transport
simulate_ember_spotting <- function(grid, wind_x, wind_y, fire_intensity) {
# Parameters
ember_prob_threshold <- 0.6 # Probability threshold for ember generation
max_spotting_distance <- 2000 # Maximum spotting distance in meters
cell_size <- 100 # Cell size in meters
max_cells <- max_spotting_distance / cell_size # Max distance in cells
# Initialize new ignition points
new_ignitions <- data.frame(x = numeric(0), y = numeric(0), intensity = numeric(0))
# Find cells with active fire
active_fire <- which(fire_intensity > 0.4)
# For each active fire cell, generate embers
for(idx in active_fire) {
# Cell coordinates
x <- grid$x[idx]
y <- grid$y[idx]
# Calculate ember generation probability based on fire intensity and fuel
ember_prob <- fire_intensity[idx] * grid$sustainability[idx] * (1 + grid$fuel_load[idx]/2)
# Determine number of embers
if(ember_prob > ember_prob_threshold) {
n_embers <- rpois(1, lambda = 3 * (ember_prob - ember_prob_threshold) * fire_intensity[idx])
# For each ember, determine landing location
for(e in 1:n_embers) {
# Sample spotting distance using exponential distribution
distance <- rexp(1, rate = 3/max_cells) # Mean distance = max_cells/3
distance <- min(distance, max_cells) # Cap at max distance
# Apply wind direction (with some randomness)
wind_dir <- atan2(wind_y[idx], wind_x[idx])
wind_speed <- sqrt(wind_x[idx]^2 + wind_y[idx]^2)
# Add random variation to direction
rand_dir <- wind_dir + rnorm(1, 0, pi/8) # Random variation in direction
# Calculate landing coordinates
dx <- distance * cos(rand_dir)
dy <- distance * sin(rand_dir)
# Apply additional drift based on wind speed
drift_factor <- wind_speed / 10 # Normalized drift
dx <- dx * (1 + drift_factor)
# New coordinates
new_x <- round(x + dx)
new_y <- round(y + dy)
# Check if within grid and not already burning
if(new_x >= 1 && new_x <= max(grid$x) &&
new_y >= 1 && new_y <= max(grid$y)) {
# Find grid index for this location
grid_idx <- which(grid$x == new_x & grid$y == new_y)
if(length(grid_idx) > 0 && fire_intensity[grid_idx] < 0.1) {
# Calculate initial intensity based on distance traveled and fuel at landing site
initial_intensity <- 0.2 * exp(-distance/20) * grid$sustainability[grid_idx]
# Add to new ignitions
new_ignitions <- rbind(new_ignitions, data.frame(
x = new_x,
y = new_y,
intensity = initial_intensity
))
}
}
}
}
}
# Remove duplicates (multiple embers landing in same cell)
if(nrow(new_ignitions) > 0) {
unique_ignitions <- aggregate(intensity ~ x + y, data = new_ignitions, max)
return(unique_ignitions)
} else {
return(NULL)
}
}
# Integrate spotting into the main simulation
modified_fire_step_with_spotting <- function(grid, t, params) {
# First, update fire using SDE model (same as before)
new_grid <- simulate_fire_step(grid, t)
# Then, simulate ember spotting
embers <- simulate_ember_spotting(
grid,
wind_x = new_grid$wind_x,
wind_y = new_grid$wind_y,
fire_intensity = new_grid$FI
)
# Process new ignitions
if(!is.null(embers) && nrow(embers) > 0) {
for(i in 1:nrow(embers)) {
# Find grid cell
idx <- which(new_grid$x == embers$x[i] & new_grid$y == embers$y[i])
if(length(idx) > 0) {
# Add new fire (maximum of existing and new ember fire)
new_grid$FI[idx] <- max(new_grid$FI[idx], embers$intensity[i])
}
}
}
return(new_grid)
}4.4 Fire-Atmosphere Interactions
Finally, let’s implement a basic fire-atmosphere coupling model.
# Fire-Atmosphere Coupling
# Function to simulate fire-induced weather modifications
update_atmosphere <- function(grid, fire_intensity) {
# Parameters for fire-atmosphere coupling
plume_height_coef <- 100 # Height coefficient (meters)
air_temp_rise_coef <- 5 # Temperature rise coefficient (°C)
convection_strength_coef <- 0.2 # Convection strength
# For each cell, calculate fire impacts on local atmosphere
for(i in 1:nrow(grid)) {
if(fire_intensity[i] > 0.1) {
# Calculate plume characteristics
plume_height <- plume_height_coef * fire_intensity[i]^1.5
temp_rise <- air_temp_rise_coef * fire_intensity[i]
# Identify cells affected by the plume (simplified as a circle)
plume_radius <- plume_height / 5 # Radius in meters
cell_radius <- plume_radius / grid$cell_size # Radius in cells
# Find cells within plume influence
x_center <- grid$x[i]
y_center <- grid$y[i]
for(j in 1:nrow(grid)) {
# Calculate distance to plume center
dist <- sqrt((grid$x[j] - x_center)^2 + (grid$y[j] - y_center)^2)
if(dist <= cell_radius) {
# Calculate influence based on distance
influence <- (1 - dist/cell_radius)^2
# Modify local weather conditions
# 1. Increase temperature
grid$temperature[j] <- grid$temperature[j] + temp_rise * influence
# 2. Decrease relative humidity
grid$relative_humidity[j] <- max(5, grid$relative_humidity[j] - temp_rise * influence)
# 3. Modify wind field (create convection)
if(dist > 0) { # Not the center cell
# Direction toward fire center
direction <- atan2(y_center - grid$y[j], x_center - grid$x[j])
# Convection speed depends on fire intensity and height
convection_speed <- convection_strength_coef * fire_intensity[i] * plume_height/1000
# Add convection component to wind
conv_x <- convection_speed * cos(direction) * influence
conv_y <- convection_speed * sin(direction) * influence
grid$wind_x[j] <- grid$wind_x[j] + conv_x
grid$wind_y[j] <- grid$wind_y[j] + conv_y
}
}
}
}
}
return(grid)
}
# Create the integrated fire-atmosphere step function
simulate_fire_atmosphere_step <- function(grid, t, params) {
# First, simulate basic fire spread
new_grid <- modified_fire_step_with_spotting(grid, t, params)
# Then, update atmosphere based on fire feedback
new_grid <- update_atmosphere(new_grid, new_grid$FI)
return(new_grid)
}4.5 Integrated Advanced Wildfire SDE Model
Next we build the entire advanced, enhanced, and integrated wildfire SDE model.
# Integrated Advanced Wildfire SDE Model
run_advanced_wildfire_simulation <- function(config, weather_file, landcover_file,
initial_ignition = NULL, ensemble_size = 1) {
# Initialize ensemble storage
ensemble_results <- list()
# Load and process external data
weather_forecast <- process_weather_forecast(weather_file, config$grid_coords)
# Run each ensemble member
for(e in 1:ensemble_size) {
# Create initial grid
grid <- create_fire_grid(config$grid_size$rows, config$grid_size$cols)
# Add geographical coordinates
grid$lat <- config$grid_coords$lat_min + (grid$y - 1) * config$grid_coords$lat_step
grid$lon <- config$grid_coords$lon_min + (grid$x - 1) * config$grid_coords$lon_step
# Assign fuel types based on land cover
grid <- assign_fuel_types(grid, landcover_file)
# Set initial ignition if provided
if(!is.null(initial_ignition)) {
ignition_points <- grid %>%
filter(x >= initial_ignition$x_min & x <= initial_ignition$x_max &
y >= initial_ignition$y_min & y <= initial_ignition$y_max)
grid$FI[grid$x >= initial_ignition$x_min & grid$x <= initial_ignition$x_max &
grid$y >= initial_ignition$y_min & grid$y <= initial_ignition$y_max] <-
initial_ignition$intensity
}
# Initialize history
history <- list()
history[[1]] <- grid
# Run simulation with advanced models
for(t in 2:config$time_steps) {
# Update weather conditions for current time
grid_with_weather <- get_weather_at_time(
weather_forecast,
history[[t-1]],
config$start_time + (t-1) * config$dt
)
# Simulate fire spread with all advanced components
history[[t]] <- simulate_fire_atmosphere_step(
grid_with_weather,
t,
config
)
# Apply any firefighting interventions
if(t >= config$firefighting_start) {
history[[t]] <- apply_firefighting(history[[t]], config$firefighting_params)
}
}
# Store ensemble member
ensemble_results[[e]] <- history
}
# Calculate ensemble statistics
if(ensemble_size > 1) {
ensemble_stats <- calculate_ensemble_statistics(ensemble_results)
return(list(
ensembles = ensemble_results,
stats = ensemble_stats
))
} else {
return(ensemble_results[[1]])
}
}
# Example configuration for advanced model
advanced_config <- list(
# Grid parameters
grid_size = list(rows = 50, cols = 50),
cell_size = 100, # meters
grid_coords = list(
lat_min = 34.5,
lon_min = -119.8,
lat_step = 0.001, # approximately 100m
lon_step = 0.001 # approximately 100m at this latitude
),
# Fire dynamics parameters
alphaP = 0.15, # Base fire spread rate
betaP = 0.2, # Fuel influence factor
gammaP = 0.25, # Slope influence factor
deltaP = 0.3, # Heat influence
epsilonP = 0.12, # Fuel consumption rate
# Weather effect parameters
gamma_wind = 0.2, # Wind influence
delta_temp = 0.05, # Temperature influence
rh_effect = 0.8, # Humidity influence
rain_suppression = 0.5, # Rain suppression effect
temp_effect = 0.03, # Temperature effect on fuel
humidity_effect = 0.7, # Humidity effect on fuel
precipitation_effect = 0.9, # Precipitation effect
# Stochastic parameters
sigma_FI = 0.04, # Fire intensity noise
sigma_FC = 0.06, # Fuel consumption noise
# Simulation parameters
time_steps = 72, # Run for 72 hours
dt = 1, # 1 hour time steps
start_time = as.POSIXct("2024-04-15 12:00:00", tz = "UTC"),
# Ember spotting parameters
spotting_params = list(
ember_prob_threshold = 0.6,
max_spotting_distance = 2000
),
# Fire-atmosphere parameters
atmosphere_params = list(
plume_height_coef = 100,
air_temp_rise_coef = 5,
convection_strength_coef = 0.2
),
# Firefighting parameters
firefighting_start = 12, # Start firefighting after 12 hours
firefighting_params = list(
strategy = "perimeter",
resources = 10,
effectiveness = 0.3
)
)
# Apply firefighting interventions
apply_firefighting <- function(grid, firefighting_params) {
# Implement firefighting logic based on strategy
if(firefighting_params$strategy == "perimeter") {
# Find fire perimeter cells
fire_cells <- which(grid$FI > 0.1)
perimeter_cells <- c()
for(idx in fire_cells) {
x <- grid$x[idx]
y <- grid$y[idx]
# Check if this is a perimeter cell (has at least one non-burning neighbor)
neighbors <- which(
grid$x >= x-1 & grid$x <= x+1 &
grid$y >= y-1 & grid$y <= y+1 &
!(grid$x == x & grid$y == y)
)
if(any(grid$FI[neighbors] < 0.1)) {
perimeter_cells <- c(perimeter_cells, idx)
}
}
# Prioritize perimeter cells by fire intensity
if(length(perimeter_cells) > 0) {
perimeter_df <- data.frame(
idx = perimeter_cells,
FI = grid$FI[perimeter_cells]
)
perimeter_df <- perimeter_df[order(-perimeter_df$FI),]
# Apply suppression to top priority cells based on available resources
n_target <- min(firefighting_params$resources, nrow(perimeter_df))
target_cells <- perimeter_df$idx[1:n_target]
# Reduce fire intensity in target cells
grid$FI[target_cells] <- grid$FI[target_cells] * (1 - firefighting_params$effectiveness)
}
}
return(grid)
}4.6 Advanced Visualization for the Integrated Model
Let’s showcase some of the solutions of the advanced integrated SDE model.
# Advanced visualization to show all components
visualize_advanced_fire_model <- function(simulation_results, time_step) {
# Extract data for the specific time step
grid_data <- simulation_results[[time_step]]
# Create layout for multiple panels
par(mfrow = c(2, 3), mar = c(3, 3, 2, 1), oma = c(0, 0, 2, 0))
# 1. Fire Intensity with Wind Vectors
# Create matrices for visualization
FI_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FI_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FI[i]
}
# Plot fire intensity
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = "Fire Intensity & Wind",
xlab = "", ylab = "")
# Add wind vectors (subsample for clarity)
subsample <- seq(1, nrow(grid_data), by = 5)
arrows(
grid_data$x[subsample],
grid_data$y[subsample],
grid_data$x[subsample] + grid_data$wind_x[subsample] * 2,
grid_data$y[subsample] + grid_data$wind_y[subsample] * 2,
length = 0.05, col = "white"
)
# 2. Fuel Consumption
FC_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FC_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FC[i]
}
image(1:nrow(FC_matrix), 1:ncol(FC_matrix), FC_matrix,
col = colorRampPalette(c("white", "lightgreen", "darkgreen"))(100),
main = "Fuel Remaining",
xlab = "", ylab = "")
# 3. Fuel Types
fuel_type_numeric <- as.numeric(factor(grid_data$fuel_type))
fuel_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
fuel_matrix[grid_data$x[i], grid_data$y[i]] <- fuel_type_numeric[i]
}
image(1:nrow(fuel_matrix), 1:ncol(fuel_matrix), fuel_matrix,
col = c("yellow", "orange", "darkgreen", "brown"),
main = "Fuel Types",
xlab = "", ylab = "")
legend("bottomright",
legend = c("Grass", "Shrub", "Forest", "Slash"),
fill = c("yellow", "orange", "darkgreen", "brown"),
cex = 0.6, bty = "n")
# 4. Temperature with Fire-Atmosphere Effects
temp_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
temp_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$temperature[i]
}
image(1:nrow(temp_matrix), 1:ncol(temp_matrix), temp_matrix,
col = colorRampPalette(c("blue", "green", "yellow", "orange", "red"))(100),
main = "Temperature (°C)",
xlab = "", ylab = "")
# 5. Relative Humidity
rh_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
rh_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$relative_humidity[i]
}
image(1:nrow(rh_matrix), 1:ncol(rh_matrix), rh_matrix,
col = colorRampPalette(c("red", "orange", "yellow", "lightblue", "blue"))(100),
main = "Relative Humidity (%)",
xlab = "", ylab = "")
# 6. Spotting Activity (New Ignitions)
# Identify spot fires - cells that weren't burning in previous time step
if(time_step > 1) {
prev_grid <- simulation_results[[time_step-1]]
spot_fires <- data.frame(
x = numeric(0),
y = numeric(0)
)
for(i in 1:nrow(grid_data)) {
if(grid_data$FI[i] > 0.1) {
prev_idx <- which(prev_grid$x == grid_data$x[i] & prev_grid$y == grid_data$y[i])
if(length(prev_idx) > 0 && prev_grid$FI[prev_idx] < 0.1) {
# This is a new ignition
spot_fires <- rbind(spot_fires, data.frame(
x = grid_data$x[i],
y = grid_data$y[i]
))
}
}
}
# Create visualization base (copy of fire intensity)
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = "Spotting Activity",
xlab = "", ylab = "")
# Add spot fire markers
if(nrow(spot_fires) > 0) {
points(spot_fires$x, spot_fires$y, pch = 16, col = "white", cex = 1.2)
points(spot_fires$x, spot_fires$y, pch = 8, col = "black", cex = 1)
}
} else {
# For first time step, just show fire intensity
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = "Initial Fire State",
xlab = "", ylab = "")
}
# Add overall title
mtext(paste("Advanced Wildfire Model - Time Step", time_step),
outer = TRUE, cex = 1.5)
}
# Function to create 3D visualization of wildfire with atmospheric effects
visualize_fire_atmosphere_3D <- function(simulation_results, time_step) {
# Extract data for the specific time step
grid_data <- simulation_results[[time_step]]
# Create matrices
FI_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
temp_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
elevation_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FI_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FI[i]
temp_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$temperature[i]
elevation_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$elevation[i]
}
# Normalize temperature for color mapping
temp_norm <- (temp_matrix - min(temp_matrix)) / (max(temp_matrix) - min(temp_matrix))
# Create 3D plot using plotly
library(plotly)
# Base terrain
fig <- plot_ly(
z = elevation_matrix * 10, # Scale elevation for visibility
type = "surface",
colorscale = list(c(0, 1), c("tan", "darkgreen")),
showscale = FALSE,
opacity = 0.8
)
# Add fire as a raised surface where intensity > 0.1
fire_height <- elevation_matrix * 10 # Same as base terrain
fire_color <- matrix(NA, nrow = nrow(FI_matrix), ncol = ncol(FI_matrix))
for(i in 1:nrow(FI_matrix)) {
for(j in 1:ncol(FI_matrix)) {
if(FI_matrix[i,j] > 0.1) {
# Height proportional to fire intensity
fire_height[i,j] <- fire_height[i,j] + FI_matrix[i,j] * 3
# Color based on intensity (red to black)
fire_color[i,j] <- FI_matrix[i,j]
}
}
}
fig <- fig %>% add_surface(
z = fire_height,
surfacecolor = fire_color,
colorscale = list(
c(0, 0.5, 1),
c("yellow", "red", "black")
),
cmin = 0.1,
cmax = 1,
showscale = TRUE,
colorbar = list(title = "Fire Intensity"),
opacity = 0.9
)
# Add smoke plumes using cone shapes
# Find cells with significant fire
hot_spots <- which(FI_matrix > 0.5, arr.ind = TRUE)
if(nrow(hot_spots) > 0) {
for(i in 1:min(nrow(hot_spots), 10)) { # Limit to 10 plumes for performance
x_pos <- hot_spots[i, 1]
y_pos <- hot_spots[i, 2]
# Intensity at this position
intensity <- FI_matrix[x_pos, y_pos]
# Plume height based on intensity
plume_height <- intensity * 15
# Create smoke particles as scattered points
n_particles <- 100
# Wind direction at this location
grid_idx <- which(grid_data$x == x_pos & grid_data$y == y_pos)
wind_x <- grid_data$wind_x[grid_idx]
wind_y <- grid_data$wind_y[grid_idx]
# Generate random particles in a cone shape
heights <- runif(n_particles, 0, plume_height)
angles <- runif(n_particles, 0, 2*pi)
radii <- runif(n_particles, 0, heights/5) # Cone shape
# Adjust for wind
wind_factor <- sqrt(wind_x^2 + wind_y^2) * 2
wind_angle <- atan2(wind_y, wind_x)
x_offset <- radii * cos(angles) + heights * wind_factor * cos(wind_angle)
y_offset <- radii * sin(angles) + heights * wind_factor * sin(wind_angle)
# Add particles
fig <- fig %>% add_trace(
x = x_pos + x_offset,
y = y_pos + y_offset,
z = elevation_matrix[x_pos, y_pos] * 10 + heights,
type = "scatter3d",
mode = "markers",
marker = list(
color = "gray",
opacity = 0.5,
size = 3
),
showlegend = FALSE
)
}
}
# Add temperature as a semi-transparent layer
# We'll create points in a grid and color them by temperature
x_seq <- seq(1, nrow(temp_matrix), by = 2)
y_seq <- seq(1, ncol(temp_matrix), by = 2)
temp_points <- expand.grid(x = x_seq, y = y_seq)
temp_points$z <- 12 # Fixed height for temperature layer
temp_points$temp <- NA
for(i in 1:nrow(temp_points)) {
x <- temp_points$x[i]
y <- temp_points$y[i]
temp_points$temp[i] <- temp_matrix[x, y]
}
fig <- fig %>% add_trace(
x = temp_points$x,
y = temp_points$y,
z = temp_points$z,
type = "scatter3d",
mode = "markers",
marker = list(
color = temp_points$temp,
colorscale = "Hot",
opacity = 0.3,
size = 5
),
name = "Air Temperature"
)
# Add layout settings
fig <- fig %>% layout(
title = paste("3D Wildfire with Atmospheric Effects - Time Step", time_step),
scene = list(
aspectratio = list(x = 1, y = 1, z = 0.5),
xaxis = list(title = ""),
yaxis = list(title = ""),
zaxis = list(title = "Elevation"),
camera = list(eye = list(x = 1.5, y = -1.5, z = 1))
)
)
return(fig)
}
# Create advanced animation with all components
create_advanced_fire_animation <- function(simulation_results, output_file = "advanced_wildfire.gif") {
library(animation)
# Set animation options
ani.options(interval = 0.25)
# Save animation
saveGIF({
for(t in 1:length(simulation_results)) {
visualize_advanced_fire_model(simulation_results, t)
}
}, movie.name = output_file, ani.width = 1200, ani.height = 800)
}
# Create a simpler but effective animation function
create_fire_animation <- function(simulation_results, output_file = "wildfire_simple.gif") {
library(animation)
# Set animation options - use a compatible interval
ani.options(interval = 0.2) # 5 frames per second (factor of 100)
# Save animation with the corrected approach
saveGIF({
for(t in 1:length(simulation_results)) {
# Create a simple single-panel plot for each frame
# Extract the fire intensity matrix for current time step
grid_data <- simulation_results[[t]]
# Create matrices for visualization
FI_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FI_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FI[i]
}
# Simple plot with just fire intensity
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = paste("Wildfire Simulation - Hour", t-1),
xlab = "X", ylab = "Y")
# Add elevation contours if available
if("elevation" %in% names(grid_data)) {
elevation_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
elevation_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$elevation[i]
}
contour(1:nrow(elevation_matrix), 1:ncol(elevation_matrix), elevation_matrix,
add = TRUE, col = "black", lwd = 0.5)
}
}
}, movie.name = output_file, ani.width = 800, ani.height = 600, clean = TRUE)
return(output_file)
}
# Generate synthetic data for the advanced visualization
generate_synthetic_fire_data <- function(time_steps = 40, grid_size = 30) {
# Set seed for reproducibility
set.seed(123)
# Create base grid
grid <- expand.grid(x = 1:grid_size, y = 1:grid_size)
# Add synthetic elevation with hills
elevation <- matrix(0, nrow = grid_size, ncol = grid_size)
# Add some hills
for(i in 1:3) {
center_x <- runif(1, 5, grid_size-5)
center_y <- runif(1, 5, grid_size-5)
height <- runif(1, 0.3, 0.7)
width <- runif(1, 3, 6)
for(x in 1:grid_size) {
for(y in 1:grid_size) {
dist_sq <- (x - center_x)^2 + (y - center_y)^2
elevation[x, y] <- elevation[x, y] + height * exp(-dist_sq / (2 * width^2))
}
}
}
# Normalize elevation
elevation <- elevation / max(elevation)
# Add elevation to grid
for(i in 1:nrow(grid)) {
grid$elevation[i] <- elevation[grid$x[i], grid$y[i]]
}
# Assign fuel types based on elevation
grid$fuel_type <- "grass" # Default
for(i in 1:nrow(grid)) {
if(grid$elevation[i] > 0.7) {
grid$fuel_type[i] <- "forest"
} else if(grid$elevation[i] > 0.3) {
grid$fuel_type[i] <- "shrub"
}
# Add some randomness
if(runif(1) < 0.2) {
grid$fuel_type[i] <- sample(c("grass", "shrub", "forest"), 1)
}
}
# Add initial properties
grid$FI <- 0 # Fire intensity
grid$FC <- 1 # Fuel consumption (full at start)
grid$temperature <- 25 # Ambient temperature
grid$relative_humidity <- 40 # RH percentage
grid$wind_x <- 0.2 # Wind components
grid$wind_y <- 0.1
grid$precipitation <- 0 # No rain initially
# Add fuel properties based on fuel types
for(i in 1:nrow(grid)) {
if(grid$fuel_type[i] == "grass") {
grid$fuel_load[i] <- 0.3
grid$moisture_threshold[i] <- 15
grid$spread_rate_coef[i] <- 0.8
grid$intensity_coef[i] <- 0.5
grid$sustainability[i] <- 0.4
} else if(grid$fuel_type[i] == "shrub") {
grid$fuel_load[i] <- 0.8
grid$moisture_threshold[i] <- 25
grid$spread_rate_coef[i] <- 0.6
grid$intensity_coef[i] <- 0.7
grid$sustainability[i] <- 0.6
} else { # forest
grid$fuel_load[i] <- 2.5
grid$moisture_threshold[i] <- 30
grid$spread_rate_coef[i] <- 0.4
grid$intensity_coef[i] <- 0.9
grid$sustainability[i] <- 0.8
}
}
# Start fire in bottom left corner
ignition_points <- which(grid$x <= 3 & grid$y <= 3)
grid$FI[ignition_points] <- 0.5
# Create simulation results container
simulation_results <- list()
simulation_results[[1]] <- grid
# Parameters for fire spread
alpha_spread <- 0.15 # Base fire spread rate
beta_wind <- 0.2 # Wind influence factor
gamma_slope <- 0.25 # Slope influence factor
delta_fuel <- 0.3 # Fuel content influence
epsilon_burn <- 0.12 # Rate of fuel consumption
# Simulate simple fire spread for remaining time steps
for(t in 2:time_steps) {
new_grid <- simulation_results[[t-1]]
# Update weather - add a weather front moving through
if(t > 15 && t < 25) {
# Weather front arrives
new_grid$wind_x <- 0.3 + (t - 15) * 0.02
new_grid$wind_y <- 0.2 + (t - 15) * 0.01
new_grid$temperature <- 25 - (t - 15) * 0.5
new_grid$relative_humidity <- 40 + (t - 15) * 2
# Light precipitation
if(t > 18 && t < 22) {
new_grid$precipitation <- 0.2
} else {
new_grid$precipitation <- 0
}
} else if(t >= 25) {
# After front passes
new_grid$wind_x <- 0.4
new_grid$wind_y <- 0.3
new_grid$temperature <- 20
new_grid$relative_humidity <- 60
new_grid$precipitation <- 0
}
# Spread fire
for(i in 1:nrow(new_grid)) {
# Skip if cell already has maximum fire
if(new_grid$FI[i] >= 1.0) next
row <- new_grid[i,]
x <- row$x
y <- row$y
# Check if this cell or neighbors are burning
if(row$FI > 0 || any(new_grid$FI[new_grid$x >= x-1 & new_grid$x <= x+1 &
new_grid$y >= y-1 & new_grid$y <= y+1] > 0)) {
# Calculate neighboring fire influence
neighbors <- new_grid[new_grid$x >= x-1 & new_grid$x <= x+1 &
new_grid$y >= y-1 & new_grid$y <= y+1,]
neighbor_influence <- sum(neighbors$FI) / 9 # Average including self
# Calculate slope effect (fire spreads faster uphill)
if(x < grid_size && y < grid_size) {
slope_x <- new_grid$elevation[new_grid$x == x+1 & new_grid$y == y] - row$elevation
slope_y <- new_grid$elevation[new_grid$x == x & new_grid$y == y+1] - row$elevation
slope_effect <- gamma_slope * (max(0, slope_x) + max(0, slope_y))
} else {
slope_effect <- 0
}
# Wind effect
wind_effect <- beta_wind * (row$wind_x + row$wind_y)
# Calculate fire intensity change based on all factors
dFI <- alpha_spread * neighbor_influence * (1 - row$FI) * row$spread_rate_coef +
delta_fuel * row$FC * row$intensity_coef +
slope_effect + wind_effect -
0.5 * row$precipitation # Rain suppression
# Update fire intensity
new_grid$FI[i] <- min(1.0, max(0, row$FI + dFI))
# Update fuel consumption
if(new_grid$FI[i] > 0.1) {
new_grid$FC[i] <- max(0, row$FC - epsilon_burn * row$FI * row$sustainability)
}
# Add fire-atmosphere coupling - increase temperature near fire
if(new_grid$FI[i] > 0.3) {
new_grid$temperature[i] <- new_grid$temperature[i] + 10 * new_grid$FI[i]
new_grid$relative_humidity[i] <- max(10, new_grid$relative_humidity[i] - 20 * new_grid$FI[i])
}
}
}
# Add spotting (occasionally create new fire spots)
if(runif(1) < 0.2) { # 20% chance each time step
# Find cells with high intensity fire
hot_spots <- which(new_grid$FI > 0.7)
if(length(hot_spots) > 0) {
# Pick a random hot spot
source_idx <- sample(hot_spots, 1)
source_x <- new_grid$x[source_idx]
source_y <- new_grid$y[source_idx]
# Calculate ember landing spot (biased by wind)
distance <- rexp(1, rate = 0.5) * 3 + 3 # 3-6 cells typically
angle <- atan2(new_grid$wind_y[source_idx], new_grid$wind_x[source_idx]) +
rnorm(1, 0, 0.5) # Some randomness
target_x <- round(source_x + distance * cos(angle))
target_y <- round(source_y + distance * sin(angle))
# Check if valid location and not already burning
if(target_x >= 1 && target_x <= grid_size &&
target_y >= 1 && target_y <= grid_size) {
target_idx <- which(new_grid$x == target_x & new_grid$y == target_y)
if(length(target_idx) > 0 && new_grid$FI[target_idx] < 0.1 &&
new_grid$FC[target_idx] > 0.3) {
# Create new fire spot
new_grid$FI[target_idx] <- 0.3 # Initial ember intensity
}
}
}
}
simulation_results[[t]] <- new_grid
}
return(simulation_results)
}
# Create proper visualization function for advanced model
visualize_advanced_fire_model <- function(simulation_results, time_step) {
# Extract data for the specific time step
grid_data <- simulation_results[[time_step]]
# Create layout for multiple panels
par(mfrow = c(2, 3), mar = c(3, 3, 2, 1), oma = c(0, 0, 2, 0))
# 1. Fire Intensity with Wind Vectors
# Create matrices for visualization
FI_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FI_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FI[i]
}
# Plot fire intensity
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = "Fire Intensity & Wind",
xlab = "", ylab = "")
# Add wind vectors (subsample for clarity)
subsample <- seq(1, nrow(grid_data), by = 5)
arrows(
grid_data$x[subsample],
grid_data$y[subsample],
grid_data$x[subsample] + grid_data$wind_x[subsample] * 2,
grid_data$y[subsample] + grid_data$wind_y[subsample] * 2,
length = 0.05, col = "white"
)
# 2. Fuel Consumption
FC_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FC_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FC[i]
}
image(1:nrow(FC_matrix), 1:ncol(FC_matrix), FC_matrix,
col = colorRampPalette(c("white", "lightgreen", "darkgreen"))(100),
main = "Fuel Remaining",
xlab = "", ylab = "")
# 3. Fuel Types
fuel_type_numeric <- as.numeric(factor(grid_data$fuel_type))
fuel_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
fuel_matrix[grid_data$x[i], grid_data$y[i]] <- fuel_type_numeric[i]
}
image(1:nrow(fuel_matrix), 1:ncol(fuel_matrix), fuel_matrix,
col = c("yellow", "orange", "darkgreen", "brown")[1:length(unique(fuel_type_numeric))],
main = "Fuel Types",
xlab = "", ylab = "")
legend("bottomright",
legend = unique(grid_data$fuel_type),
fill = c("yellow", "orange", "darkgreen", "brown")[1:length(unique(fuel_type_numeric))],
cex = 0.6, bty = "n")
# 4. Temperature with Fire-Atmosphere Effects
temp_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
temp_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$temperature[i]
}
image(1:nrow(temp_matrix), 1:ncol(temp_matrix), temp_matrix,
col = colorRampPalette(c("blue", "green", "yellow", "orange", "red"))(100),
main = "Temperature (°C)",
xlab = "", ylab = "")
# 5. Relative Humidity
rh_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
rh_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$relative_humidity[i]
}
image(1:nrow(rh_matrix), 1:ncol(rh_matrix), rh_matrix,
col = colorRampPalette(c("red", "orange", "yellow", "lightblue", "blue"))(100),
main = "Relative Humidity (%)",
xlab = "", ylab = "")
# 6. Spotting Activity (New Ignitions) or Elevation for first frame
if(time_step > 1) {
prev_grid <- simulation_results[[time_step-1]]
# Plot fire intensity again as base
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = "Spotting Activity",
xlab = "", ylab = "")
# Find cells that weren't burning in previous time step but are now
spot_fires <- data.frame(
x = numeric(0),
y = numeric(0)
)
for(i in 1:nrow(grid_data)) {
if(grid_data$FI[i] > 0.1) {
prev_idx <- which(prev_grid$x == grid_data$x[i] & prev_grid$y == grid_data$y[i])
if(length(prev_idx) > 0 && prev_grid$FI[prev_idx] < 0.1) {
# This is a new ignition
spot_fires <- rbind(spot_fires, data.frame(
x = grid_data$x[i],
y = grid_data$y[i]
))
}
}
}
# Add spot fire markers if any
if(nrow(spot_fires) > 0) {
points(spot_fires$x, spot_fires$y, pch = 16, col = "white", cex = 1.2)
points(spot_fires$x, spot_fires$y, pch = 8, col = "black", cex = 1)
}
} else {
# For first time step, show elevation
elevation_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
elevation_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$elevation[i]
}
image(1:nrow(elevation_matrix), 1:ncol(elevation_matrix), elevation_matrix,
col = terrain.colors(100),
main = "Elevation",
xlab = "", ylab = "")
contour(1:nrow(elevation_matrix), 1:ncol(elevation_matrix), elevation_matrix,
add = TRUE, col = "black", lwd = 0.5)
}
# Add overall title
mtext(paste("Advanced Wildfire Model - Time Step", time_step),
outer = TRUE, cex = 1.5)
}
# Now create the animation with synthetic data
generate_and_animate_advanced_model <- function() {
# Generate synthetic data
simulation_data <- generate_synthetic_fire_data(time_steps = 40, grid_size = 30)
# Create animation with the complete data
create_advanced_fire_animation(simulation_data, "advanced_wildfire_simulation.gif")
return("Animation created: advanced_wildfire_simulation.gif")
}
# Run this function to generate and visualize the complete advanced model
result <- generate_and_animate_advanced_model()
result## [1] "Animation created: advanced_wildfire_simulation.gif"
This simulation
shows generating a complete synthetic dataset that contains all
the required variables and sets appropriate animation parameters.
Specifically, the synthetically generated data incorporates all these
advanced features
Weather integration with a front passing through
Different fuel types with varying properties
Ember spotting behavior
Fire-atmosphere interactions that modify local temperature and humidity.
The animation depicts 6 graph panels with time dynamics of the SDE wildfire model.
4.7 Advanced Model Integration with Real-World Data
We can incorporate real-world data into the advanced SDE wildfire model.
Below, we show a more dynamic and non-linear wildfire simulation by modifying the earlier simulation parameters to create more complex behavior throughout the entire time periodm \(0\leq t\leq 40\).
#################################################################
# FIRST
# USING SIM/SYNTH DATA ....
# Enhanced synthetic data generation with more dynamic fire behavior
generate_dynamic_fire_data <- function(time_steps = 40, grid_size = 30) {
# Set seed for reproducibility
set.seed(123)
# Create base grid
grid <- expand.grid(x = 1:grid_size, y = 1:grid_size)
# Add synthetic elevation with more complex terrain
elevation <- matrix(0, nrow = grid_size, ncol = grid_size)
# Create a ridgeline running diagonally
for(x in 1:grid_size) {
for(y in 1:grid_size) {
# Ridge along diagonal
ridge_dist <- abs((x-y)/sqrt(2))
elevation[x, y] <- elevation[x, y] + 0.6 * exp(-ridge_dist^2 / 18)
}
}
# Add several hills with varying heights
for(i in 1:5) {
center_x <- runif(1, 5, grid_size-5)
center_y <- runif(1, 5, grid_size-5)
height <- runif(1, 0.2, 0.8)
width <- runif(1, 3, 8)
for(x in 1:grid_size) {
for(y in 1:grid_size) {
dist_sq <- (x - center_x)^2 + (y - center_y)^2
elevation[x, y] <- elevation[x, y] + height * exp(-dist_sq / (2 * width^2))
}
}
}
# Add a river valley
river_x <- round(grid_size * 0.7)
river_width <- 1.5
for(x in 1:grid_size) {
for(y in 1:grid_size) {
river_dist <- abs(x - river_x)
valley_effect <- -0.3 * exp(-river_dist^2 / (2 * river_width^2))
elevation[x, y] <- elevation[x, y] + valley_effect
}
}
# Normalize elevation
elevation <- (elevation - min(elevation)) / (max(elevation) - min(elevation))
# Add elevation to grid
for(i in 1:nrow(grid)) {
grid$elevation[i] <- elevation[grid$x[i], grid$y[i]]
}
# Assign fuel types based on elevation and patterns
grid$fuel_type <- "grass" # Default
# Create a moisture gradient (wetter near river)
grid$moisture <- 0.5 # Base moisture
for(i in 1:nrow(grid)) {
# Higher elevation = drier
elev_effect <- -0.3 * grid$elevation[i]
# Distance from river = drier
river_dist <- abs(grid$x[i] - river_x)
river_effect <- -0.2 * (1 - exp(-river_dist^2 / 25))
# Combine effects
grid$moisture[i] <- 0.5 + elev_effect + river_effect
grid$moisture[i] <- max(0.1, min(0.9, grid$moisture[i])) # Bound between 0.1-0.9
}
# Assign fuel types based on elevation and moisture
for(i in 1:nrow(grid)) {
# Forest on higher, wetter areas
if(grid$elevation[i] > 0.6 && grid$moisture[i] > 0.4) {
grid$fuel_type[i] <- "forest"
}
# Shrub on medium elevations or drier high areas
else if(grid$elevation[i] > 0.3 || (grid$elevation[i] > 0.6 && grid$moisture[i] <= 0.4)) {
grid$fuel_type[i] <- "shrub"
}
# Grass in low, dry areas
else if(grid$moisture[i] < 0.3) {
grid$fuel_type[i] <- "grass"
}
# Slash in disturbed areas (add some random patches)
if(runif(1) < 0.05 && grid$x[i] < grid_size*0.6) { # Slash on left side of map
grid$fuel_type[i] <- "slash"
}
}
# Create a firebreak (road or river)
firebreak_cells <- which(abs(grid$x - river_x) < 1)
grid$fuel_type[firebreak_cells] <- "none"
# Add initial properties
grid$FI <- 0 # Fire intensity
grid$FC <- 1 # Fuel consumption (full at start)
grid$temperature <- 25 # Ambient temperature
grid$relative_humidity <- 40 # RH percentage
grid$wind_x <- 0.2 # Wind components
grid$wind_y <- 0.1
grid$precipitation <- 0 # No rain initially
# Add fuel properties based on fuel types
for(i in 1:nrow(grid)) {
if(grid$fuel_type[i] == "grass") {
grid$fuel_load[i] <- 0.3 * (1 - 0.5 * grid$moisture[i]) # Less fuel when wet
grid$moisture_threshold[i] <- 15
grid$spread_rate_coef[i] <- 0.8
grid$intensity_coef[i] <- 0.5
grid$sustainability[i] <- 0.4
} else if(grid$fuel_type[i] == "shrub") {
grid$fuel_load[i] <- 0.8 * (1 - 0.3 * grid$moisture[i])
grid$moisture_threshold[i] <- 25
grid$spread_rate_coef[i] <- 0.6
grid$intensity_coef[i] <- 0.7
grid$sustainability[i] <- 0.6
} else if(grid$fuel_type[i] == "forest") {
grid$fuel_load[i] <- 2.5 * (1 - 0.2 * grid$moisture[i])
grid$moisture_threshold[i] <- 30
grid$spread_rate_coef[i] <- 0.4
grid$intensity_coef[i] <- 0.9
grid$sustainability[i] <- 0.8
} else if(grid$fuel_type[i] == "slash") {
grid$fuel_load[i] <- 3.5 * (1 - 0.2 * grid$moisture[i])
grid$moisture_threshold[i] <- 25
grid$spread_rate_coef[i] <- 0.5
grid$intensity_coef[i] <- 1.0
grid$sustainability[i] <- 0.9
} else { # none (firebreak)
grid$fuel_load[i] <- 0.01
grid$moisture_threshold[i] <- 99
grid$spread_rate_coef[i] <- 0.01
grid$intensity_coef[i] <- 0.01
grid$sustainability[i] <- 0.01
}
}
# Start multiple ignition points for more complex fire behavior
# Initial fire in bottom left
ignition1 <- which(grid$x <= 3 & grid$y <= 3)
grid$FI[ignition1] <- 0.5
# Create simulation results container
simulation_results <- list()
simulation_results[[1]] <- grid
# Parameters for fire spread with higher variability
alpha_spread <- 0.2 # Base fire spread rate
beta_wind <- 0.25 # Wind influence factor
gamma_slope <- 0.3 # Slope influence factor - increased for more terrain effect
delta_fuel <- 0.35 # Fuel content influence
epsilon_burn <- 0.15 # Rate of fuel consumption
sigma_stochastic <- 0.08 # Stochastic component - increased for more variability
# Create dynamic weather shifts for the simulation period
# Wind shifts will create more interesting fire patterns
wind_directions <- list()
for(t in 1:time_steps) {
# Create several weather regimes during the simulation
if(t < 10) {
# Initial winds from southwest
wind_directions[[t]] <- list(
x = 0.2 + 0.05 * sin(t/3),
y = 0.1 + 0.03 * sin(t/2)
)
} else if(t < 20) {
# Wind shifts to westerly and strengthens
wind_directions[[t]] <- list(
x = 0.4 + 0.1 * sin(t/4),
y = -0.05 + 0.08 * sin(t/3)
)
} else if(t < 30) {
# Wind shifts to northwesterly
wind_directions[[t]] <- list(
x = 0.3 + 0.07 * sin(t/3),
y = -0.25 + 0.05 * sin(t/5)
)
} else {
# Wind shifts to northerly
wind_directions[[t]] <- list(
x = 0.1 + 0.05 * sin(t/4),
y = -0.4 + 0.06 * sin(t/3)
)
}
}
# Track where we've added delayed ignitions
delayed_ignitions <- list()
# Simulate fire spread for remaining time steps
for(t in 2:time_steps) {
new_grid <- simulation_results[[t-1]]
# Update weather - with more shifts throughout simulation
# Apply wind from pre-calculated patterns
new_grid$wind_x <- wind_directions[[t-1]]$x
new_grid$wind_y <- wind_directions[[t-1]]$y
# Temperature and humidity changes
if(t >= 10 && t < 20) {
# Weather front arrives - cooler and more humid
new_grid$temperature <- 25 - (t - 10) * 0.5
new_grid$relative_humidity <- 40 + (t - 10) * 2
# Light precipitation if in middle of front
if(t >= 13 && t <= 16) {
new_grid$precipitation <- 0.3
} else {
new_grid$precipitation <- 0
}
} else if(t >= 20 && t < 30) {
# Front clears, hot and dry conditions return
new_grid$temperature <- 20 + (t - 20) * 0.7
new_grid$relative_humidity <- 60 - (t - 20) * 3
new_grid$precipitation <- 0
} else if(t >= 30) {
# Very hot and dry - extreme fire weather
new_grid$temperature <- 28 + (t - 30) * 0.2
new_grid$relative_humidity <- max(10, 30 - (t - 30))
new_grid$precipitation <- 0
}
# Add a delayed ignition point across the barrier at t=25
if(t == 25) {
# Start a new fire across the river/firebreak
new_ignition_point <- which(grid$x > river_x+2 &
grid$x < river_x+5 &
grid$y > 10 & grid$y < 15 &
grid$fuel_type != "none")
if(length(new_ignition_point) > 0) {
new_grid$FI[sample(new_ignition_point, 1)] <- 0.6
delayed_ignitions[[1]] <- t # Record when we added this
}
}
# Spread fire with enhanced dynamics
for(i in 1:nrow(new_grid)) {
# Skip if cell already has maximum fire or is a firebreak
if(new_grid$FI[i] >= 1.0 || new_grid$fuel_type[i] == "none") next
row <- new_grid[i,]
x <- row$x
y <- row$y
# Check if this cell or neighbors are burning
if(row$FI > 0 || any(new_grid$FI[new_grid$x >= x-1 & new_grid$x <= x+1 &
new_grid$y >= y-1 & new_grid$y <= y+1] > 0)) {
# Calculate neighboring fire influence - more weight to upwind neighbors
neighbors <- new_grid[new_grid$x >= x-1 & new_grid$x <= x+1 &
new_grid$y >= y-1 & new_grid$y <= y+1,]
# Calculate wind direction relative to each neighbor
neighbor_weights <- numeric(nrow(neighbors))
for(n in 1:nrow(neighbors)) {
dx <- neighbors$x[n] - x
dy <- neighbors$y[n] - y
# Skip self
if(dx == 0 && dy == 0) {
neighbor_weights[n] <- 0
next
}
# Calculate dot product with wind vector
wind_alignment <- -1 * (dx * row$wind_x + dy * row$wind_y) /
sqrt(dx^2 + dy^2) / sqrt(row$wind_x^2 + row$wind_y^2)
# Upwind neighbors have more influence
neighbor_weights[n] <- 1 + max(0, wind_alignment)
}
# Normalize weights and calculate influence
if(sum(neighbor_weights) > 0) {
neighbor_weights <- neighbor_weights / sum(neighbor_weights)
neighbor_influence <- sum(neighbors$FI * neighbor_weights)
} else {
neighbor_influence <- mean(neighbors$FI)
}
# Calculate slope effect (fire spreads faster uphill)
slope_effect <- 0
for(n in 1:nrow(neighbors)) {
if(neighbors$x[n] == x && neighbors$y[n] == y) next # Skip self
dx <- neighbors$x[n] - x
dy <- neighbors$y[n] - y
# Elevation difference
elev_diff <- neighbors$elevation[n] - row$elevation
# Only upslope matters, and effect is proportional to steepness
if(elev_diff > 0) {
# Steeper slopes have bigger effect
slope_effect <- slope_effect + gamma_slope * elev_diff * 5
}
}
# Wind effect - stronger effect with higher wind speed
wind_speed <- sqrt(row$wind_x^2 + row$wind_y^2)
wind_effect <- beta_wind * wind_speed * (1 + 0.5 * wind_speed)
# Fuel effect - varies by fuel type and moisture
fuel_effect <- delta_fuel * row$fuel_load * row$FC * row$intensity_coef *
max(0, 1 - (row$relative_humidity / row$moisture_threshold))
# Calculate fire intensity change based on all factors
dFI_deterministic <- alpha_spread * neighbor_influence * (1 - row$FI) * row$spread_rate_coef +
fuel_effect + slope_effect + wind_effect -
5 * row$precipitation # Stronger rain suppression
# Add stochastic component - more randomness around firebreaks
is_near_firebreak <- any(new_grid$fuel_type[new_grid$x >= x-1 & new_grid$x <= x+1 &
new_grid$y >= y-1 & new_grid$y <= y+1] == "none")
stochastic_factor <- sigma_stochastic
if(is_near_firebreak) stochastic_factor <- stochastic_factor * 2
dFI_stochastic <- stochastic_factor * rnorm(1, 0, max(0.01, row$FI))
# Update fire intensity
new_grid$FI[i] <- min(1.0, max(0, row$FI + dFI_deterministic + dFI_stochastic))
# Update fuel consumption - varies by intensity and fuel type
if(new_grid$FI[i] > 0.1) {
consumption_rate <- epsilon_burn * row$sustainability *
(1 + 0.5 * (new_grid$FI[i] - 0.1) / 0.9) # More consumption at higher intensity
new_grid$FC[i] <- max(0, row$FC - consumption_rate)
}
# Add fire-atmosphere coupling - increase temperature near fire
if(new_grid$FI[i] > 0.3) {
# Temperature increase proportional to fire intensity
temp_increase <- 15 * new_grid$FI[i]
humid_decrease <- 30 * new_grid$FI[i]
new_grid$temperature[i] <- new_grid$temperature[i] + temp_increase
new_grid$relative_humidity[i] <- max(5, new_grid$relative_humidity[i] - humid_decrease)
# Also affect nearby cells (simple diffusion)
for(j in 1:nrow(neighbors)) {
if(neighbors$x[j] == x && neighbors$y[j] == y) next # Skip self
neighbor_idx <- which(new_grid$x == neighbors$x[j] & new_grid$y == neighbors$y[j])
if(length(neighbor_idx) > 0) {
new_grid$temperature[neighbor_idx] <- new_grid$temperature[neighbor_idx] + 0.3 * temp_increase
new_grid$relative_humidity[neighbor_idx] <- max(5, new_grid$relative_humidity[neighbor_idx] - 0.3 * humid_decrease)
}
}
}
}
}
# Enhanced spotting - more likely with stronger winds and more intense fire
# Multiple potential spots per time step
n_spotting_attempts <- rpois(1, 2) # Average 2 spotting attempts per time step
for(spot in 1:n_spotting_attempts) {
# Probability increases with wind speed
wind_speed <- mean(sqrt(new_grid$wind_x^2 + new_grid$wind_y^2))
spot_prob <- 0.2 + 0.3 * min(1, wind_speed) # 20-50% chance based on wind
if(runif(1) < spot_prob) {
# Find cells with high intensity fire
hot_spots <- which(new_grid$FI > 0.6)
if(length(hot_spots) > 0) {
# Weight selection by fire intensity
spot_weights <- new_grid$FI[hot_spots]
spot_weights <- spot_weights / sum(spot_weights)
# Pick a source based on weights
source_idx <- sample(hot_spots, 1, prob = spot_weights)
source_x <- new_grid$x[source_idx]
source_y <- new_grid$y[source_idx]
# Calculate ember landing spot (strongly biased by wind)
# Distance depends on wind speed and fire intensity
base_distance <- 3 + rexp(1, rate = 0.2)
wind_multiplier <- 1 + 2 * wind_speed
intensity_multiplier <- 1 + new_grid$FI[source_idx]
distance <- base_distance * wind_multiplier * intensity_multiplier
# Direction is wind direction with some randomness
wind_angle <- atan2(new_grid$wind_y[source_idx], new_grid$wind_x[source_idx])
random_deviation <- rnorm(1, 0, pi/8) # More focused in wind direction
angle <- wind_angle + random_deviation
target_x <- round(source_x + distance * cos(angle))
target_y <- round(source_y + distance * sin(angle))
# Check if valid location and not already burning intensely
if(target_x >= 1 && target_x <= grid_size &&
target_y >= 1 && target_y <= grid_size) {
target_idx <- which(new_grid$x == target_x & new_grid$y == target_y)
if(length(target_idx) > 0 && new_grid$FI[target_idx] < 0.4 &&
new_grid$fuel_type[target_idx] != "none" &&
new_grid$FC[target_idx] > 0.2) {
# Create new fire spot with variable intensity
# Higher intensity when: dry fuels, strong winds, and from intense source
initial_intensity <- 0.2 + 0.3 * new_grid$FI[source_idx] *
(1 - new_grid$relative_humidity[target_idx]/100) *
min(1, wind_speed)
new_grid$FI[target_idx] <- min(0.8, max(new_grid$FI[target_idx], initial_intensity))
# Record spot fire creation for visualization
delayed_ignitions[[length(delayed_ignitions) + 1]] <- list(
time = t,
from_x = source_x,
from_y = source_y,
to_x = target_x,
to_y = target_y
)
}
}
}
}
}
# Firebreaks can occasionally fail when fire is intense and conditions are extreme
if(t >= 30) { # In the extreme conditions phase
# Find firebreak cells adjacent to intense fire
for(i in firebreak_cells) {
x <- new_grid$x[i]
y <- new_grid$y[i]
# Check neighbors for intense fire
neighbors <- which(new_grid$x >= x-1 & new_grid$x <= x+1 &
new_grid$y >= y-1 & new_grid$y <= y+1 &
new_grid$FI > 0.8)
if(length(neighbors) > 0) {
# Chance of firebreak failure increases with nearby fire intensity and wind
nearby_intensity <- max(new_grid$FI[neighbors])
wind_speed <- sqrt(new_grid$wind_x[i]^2 + new_grid$wind_y[i]^2)
failure_prob <- 0.01 * nearby_intensity * (1 + 2 * wind_speed)
if(runif(1) < failure_prob) {
# Firebreak fails - add some fuel and allow fire to spread
new_grid$fuel_type[i] <- "grass"
new_grid$fuel_load[i] <- 0.2
new_grid$FC[i] <- 1.0
new_grid$spread_rate_coef[i] <- 0.8
new_grid$intensity_coef[i] <- 0.5
new_grid$sustainability[i] <- 0.4
# Record failure for visualization
delayed_ignitions[[length(delayed_ignitions) + 1]] <- list(
time = t,
type = "firebreak_failure",
x = x,
y = y
)
}
}
}
}
simulation_results[[t]] <- new_grid
}
return(list(
results = simulation_results,
delayed_ignitions = delayed_ignitions
))
}
# Enhanced visualization function to show spotted fires and firebreak failures
visualize_enhanced_fire_model <- function(simulation_data, time_step) {
# Extract data for the specific time step
simulation_results <- simulation_data$results
delayed_ignitions <- simulation_data$delayed_ignitions
grid_data <- simulation_results[[time_step]]
# Create layout for multiple panels
par(mfrow = c(2, 3), mar = c(3, 3, 2, 1), oma = c(0, 0, 2, 0))
# 1. Fire Intensity with Wind Vectors
# Create matrices for visualization
FI_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FI_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FI[i]
}
# Plot fire intensity
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = "Fire Intensity & Wind",
xlab = "", ylab = "")
# Add wind vectors (subsample for clarity)
subsample <- seq(1, nrow(grid_data), by = 4)
arrows(
grid_data$x[subsample],
grid_data$y[subsample],
grid_data$x[subsample] + grid_data$wind_x[subsample] * 3,
grid_data$y[subsample] + grid_data$wind_y[subsample] * 3,
length = 0.05, col = "white"
)
# 2. Fuel Consumption
FC_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
FC_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$FC[i]
}
image(1:nrow(FC_matrix), 1:ncol(FC_matrix), FC_matrix,
col = colorRampPalette(c("black", "red", "orange", "yellow", "lightgreen", "darkgreen"))(100),
main = "Fuel Remaining",
xlab = "", ylab = "")
# 3. Fuel Types
fuel_type_numeric <- as.numeric(factor(grid_data$fuel_type))
fuel_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
fuel_matrix[grid_data$x[i], grid_data$y[i]] <- fuel_type_numeric[i]
}
fuel_colors <- c("blue", "yellow", "orange", "darkgreen", "brown")[1:length(unique(fuel_type_numeric))]
fuel_labels <- unique(grid_data$fuel_type)
image(1:nrow(fuel_matrix), 1:ncol(fuel_matrix), fuel_matrix,
col = fuel_colors,
main = "Fuel Types",
xlab = "", ylab = "")
legend("bottomright",
legend = fuel_labels,
fill = fuel_colors,
cex = 0.6, bty = "n")
# 4. Temperature with Fire-Atmosphere Effects
temp_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
temp_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$temperature[i]
}
image(1:nrow(temp_matrix), 1:ncol(temp_matrix), temp_matrix,
col = colorRampPalette(c("blue", "green", "yellow", "orange", "red"))(100),
main = "Temperature (°C)",
xlab = "", ylab = "")
# 5. Relative Humidity
rh_matrix <- matrix(0, nrow = max(grid_data$x), ncol = max(grid_data$y))
for(i in 1:nrow(grid_data)) {
rh_matrix[grid_data$x[i], grid_data$y[i]] <- grid_data$relative_humidity[i]
}
image(1:nrow(rh_matrix), 1:ncol(rh_matrix), rh_matrix,
col = colorRampPalette(c("red", "orange", "yellow", "lightblue", "blue"))(100),
main = "Relative Humidity (%)",
xlab = "", ylab = "")
# 6. Fire Progression and Spotting
# Show the base fire intensity again
image(1:nrow(FI_matrix), 1:ncol(FI_matrix), FI_matrix,
col = colorRampPalette(c("darkgreen", "yellow", "orange", "red", "darkred"))(100),
main = "Fire Progression & Spotting",
xlab = "", ylab = "")
# Add contour for previous time steps' fire perimeter
if(time_step > 1) {
# Get previous steps at intervals
interval <- max(1, floor(time_step / 5))
prev_steps <- seq(1, time_step-1, by = interval)
for(prev_t in prev_steps) {
prev_grid <- simulation_results[[prev_t]]
# Create perimeter matrix
prev_FI <- matrix(0, nrow = max(prev_grid$x), ncol = max(prev_grid$y))
for(i in 1:nrow(prev_grid)) {
prev_FI[prev_grid$x[i], prev_grid$y[i]] <- prev_grid$FI[i]
}
# Add contour
contour(1:nrow(prev_FI), 1:ncol(prev_FI), prev_FI,
levels = 0.2, add = TRUE,
drawlabels = FALSE,
col = adjustcolor("white", alpha.f = 0.5),
lwd = 0.5)
}
}
# Add spot fire indicators and firebreak failures
if(length(delayed_ignitions) > 0) {
for(i in 1:length(delayed_ignitions)) {
ignition <- delayed_ignitions[[i]]
# Only show ignitions up to current time step
if(is.list(ignition) && "time" %in% names(ignition)) {
if(ignition$time <= time_step) {
if("type" %in% names(ignition) && ignition$type == "firebreak_failure") {
# Mark firebreak failure
points(ignition$x, ignition$y, pch = 15, col = "white", cex = 1.2)
points(ignition$x, ignition$y, pch = 0, col = "black", cex = 1.2)
} else if(all(c("from_x", "from_y", "to_x", "to_y") %in% names(ignition))) {
# Mark spot fire with arrow
arrows(ignition$from_x, ignition$from_y,
ignition$to_x, ignition$to_y,
length = 0.1, col = "white", lwd = 1
)
# Circle the ember landing spot
points(ignition$to_x, ignition$to_y, pch = 16, col = "white", cex = 1.2)
points(ignition$to_x, ignition$to_y, pch = 8, col = "black", cex = 1)
} else if(is.numeric(ignition)) {
# Simple marker for older ignition format
text(10, 2, paste("New fire at t =", ignition), col = "white", cex = 0.8)
}
}
}
}
}
# Add overall title
mtext(paste("Advanced Wildfire Model - Time Step", time_step),
outer = TRUE, cex = 1.5)
}
# Create the enhanced animation function
create_enhanced_fire_animation <- function(simulation_data, output_file = "enhanced_wildfire.gif") {
library(animation)
# Set animation options with correct frame rate
ani.options(interval = 0.2) # 5 frames per second
# Save animation with the corrected approach
saveGIF({
for(t in 1:length(simulation_data$results)) {
visualize_enhanced_fire_model(simulation_data, t)
}
}, movie.name = output_file, ani.width = 1200, ani.height = 800, clean = TRUE)
return(output_file)
}
# Run the enhanced simulation and create animation
generate_and_animate_enhanced_model <- function() {
# Generate enhanced synthetic data
cat("Generating enhanced wildfire simulation data...\n")
simulation_data <- generate_dynamic_fire_data(time_steps = 40, grid_size = 30)
cat("Creating animation with the enhanced data...\n")
output_file <- create_enhanced_fire_animation(simulation_data, "enhanced_wildfire_simulation.gif")
# Bonus: Create diagnostic plots to show the dynamic nature
create_diagnostic_plots(simulation_data)
return(paste("Animation created:", output_file))
}
# Create diagnostic plots to highlight the dynamic behavior
create_diagnostic_plots <- function(simulation_data) {
# Extract data and calculate metrics over time
results <- simulation_data$results
time_steps <- length(results)
# Initialize data frames for metrics with correct length
metrics <- data.frame(
time = 1:time_steps,
burned_area = numeric(time_steps),
avg_intensity = numeric(time_steps),
max_intensity = numeric(time_steps),
spread_rate = numeric(time_steps),
active_perimeter = numeric(time_steps)
)
# Create weather data frame separately with the correct length
metrics$weather <- data.frame(
temperature = numeric(time_steps),
humidity = numeric(time_steps),
wind_speed = numeric(time_steps)
)
# Calculate metrics for each time step
for(t in 1:time_steps) {
grid <- results[[t]]
# Burned area
metrics$burned_area[t] <- sum(grid$FI > 0.1)
# Fire intensity
fire_cells <- which(grid$FI > 0.1)
if(length(fire_cells) > 0) {
metrics$avg_intensity[t] <- mean(grid$FI[fire_cells])
metrics$max_intensity[t] <- max(grid$FI)
}
# Active perimeter length
perimeter_cells <- 0
for(i in which(grid$FI > 0.1)) {
x <- grid$x[i]
y <- grid$y[i]
# Check if cell has at least one non-burning neighbor
neighbors <- which(
grid$x >= x-1 & grid$x <= x+1 &
grid$y >= y-1 & grid$y <= y+1 &
!(grid$x == x & grid$y == y)
)
if(length(neighbors) > 0 && any(grid$FI[neighbors] <= 0.1)) {
perimeter_cells <- perimeter_cells + 1
}
}
metrics$active_perimeter[t] <- perimeter_cells
# Spread rate
if(t > 1) {
metrics$spread_rate[t] <- metrics$burned_area[t] - metrics$burned_area[t-1]
}
# Average weather conditions
metrics$weather$temperature[t] <- mean(grid$temperature)
metrics$weather$humidity[t] <- mean(grid$relative_humidity)
metrics$weather$wind_speed[t] <- mean(sqrt(grid$wind_x^2 + grid$wind_y^2))
}
# Rest of the function remains the same...
# Create plots
png("wildfire_dynamic_metrics.png", width = 1200, height = 800)
par(mfrow = c(2, 2), mar = c(4, 4, 3, 4))
# Plot 1: Fire Growth and Spread Rate
plot(metrics$time, metrics$burned_area, type = "l", col = "orange", lwd = 2,
xlab = "Time Step", ylab = "Burned Area (cells)",
main = "Fire Growth and Spread Rate")
par(new = TRUE)
plot(metrics$time, metrics$spread_rate, type = "l", col = "red", lwd = 2,
xlab = "", ylab = "", axes = FALSE)
axis(4, col = "red", col.axis = "red")
mtext("Spread Rate (cells/step)", side = 4, line = 2, col = "red")
legend("topleft", legend = c("Burned Area", "Spread Rate"),
col = c("orange", "red"), lwd = 2, bty = "n")
# Plot 2: Fire Intensity
plot(metrics$time, metrics$avg_intensity, type = "l", col = "darkred", lwd = 2,
xlab = "Time Step", ylab = "Fire Intensity",
main = "Fire Intensity Dynamics", ylim = c(0, 1))
lines(metrics$time, metrics$max_intensity, col = "red", lwd = 2, lty = 2)
legend("topleft", legend = c("Average Intensity", "Maximum Intensity"),
col = c("darkred", "red"), lwd = 2, lty = c(1, 2), bty = "n")
# Plot 3: Weather Conditions
plot(metrics$time, metrics$weather$temperature, type = "l", col = "red", lwd = 2,
xlab = "Time Step", ylab = "Temperature (°C)",
main = "Weather Conditions", ylim = c(0, max(metrics$weather$temperature) * 1.1))
par(new = TRUE)
plot(metrics$time, metrics$weather$humidity, type = "l", col = "blue", lwd = 2,
xlab = "", ylab = "", axes = FALSE,
ylim = c(0, max(metrics$weather$humidity) * 1.1))
axis(4, col = "blue", col.axis = "blue")
mtext("Relative Humidity (%)", side = 4, line = 2, col = "blue")
legend("topright", legend = c("Temperature", "Humidity"),
col = c("red", "blue"), lwd = 2, bty = "n")
# Plot 4: Active Fire Perimeter and Wind Speed
plot(metrics$time, metrics$active_perimeter, type = "l", col = "purple", lwd = 2,
xlab = "Time Step", ylab = "Active Perimeter (cells)",
main = "Fire Perimeter and Wind Speed")
par(new = TRUE)
plot(metrics$time, metrics$weather$wind_speed, type = "l", col = "cyan4", lwd = 2,
xlab = "", ylab = "", axes = FALSE)
axis(4, col = "cyan4", col.axis = "cyan4")
mtext("Wind Speed", side = 4, line = 2, col = "cyan4")
legend("topleft", legend = c("Active Perimeter", "Wind Speed"),
col = c("purple", "cyan4"), lwd = 2, bty = "n")
dev.off()
# Create a second diagnostic plot showing fire behavior by fuel type
png("wildfire_fuel_type_analysis.png", width = 1200, height = 600)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
# Calculate burned area by fuel type
fuel_types <- unique(results[[1]]$fuel_type)
n_fuels <- length(fuel_types)
burned_by_type <- matrix(0, nrow = time_steps, ncol = n_fuels)
colnames(burned_by_type) <- fuel_types
for(t in 1:time_steps) {
grid <- results[[t]]
for(f in 1:n_fuels) {
fuel_type <- fuel_types[f]
fuel_cells <- which(grid$fuel_type == fuel_type)
if(length(fuel_cells) > 0) {
burned_by_type[t, f] <- sum(grid$FI[fuel_cells] > 0.1)
}
}
}
# Plot burned area by fuel type
matplot(1:time_steps, burned_by_type, type = "l", lty = 1:n_fuels,
col = rainbow(n_fuels), lwd = 2,
xlab = "Time Step", ylab = "Burned Area (cells)",
main = "Fire Spread by Fuel Type")
legend("topleft", legend = fuel_types, col = rainbow(n_fuels),
lty = 1:n_fuels, lwd = 2, bty = "n")
# Calculate average fire intensity by fuel type
intensity_by_type <- matrix(0, nrow = time_steps, ncol = n_fuels)
colnames(intensity_by_type) <- fuel_types
for(t in 1:time_steps) {
grid <- results[[t]]
for(f in 1:n_fuels) {
fuel_type <- fuel_types[f]
fuel_cells <- which(grid$fuel_type == fuel_type & grid$FI > 0.1)
if(length(fuel_cells) > 0) {
intensity_by_type[t, f] <- mean(grid$FI[fuel_cells])
}
}
}
# Plot intensity by fuel type
matplot(1:time_steps, intensity_by_type, type = "l", lty = 1:n_fuels,
col = rainbow(n_fuels), lwd = 2,
xlab = "Time Step", ylab = "Average Fire Intensity",
main = "Fire Intensity by Fuel Type")
legend("topleft", legend = fuel_types, col = rainbow(n_fuels),
lty = 1:n_fuels, lwd = 2, bty = "n")
dev.off()
return(c("wildfire_dynamic_metrics.png", "wildfire_fuel_type_analysis.png"))
}
# Run the enhanced simulation and create the animation
result <- generate_and_animate_enhanced_model()## Generating enhanced wildfire simulation data...
## Creating animation with the enhanced data...## [1] "Animation created: enhanced_wildfire_simulation.gif"#################################################################
# SECOND
# FOR USING REAL DATA ....Functions to integrate real-world data
# 1. Load actual weather forecast from NOAA
load_noaa_weather_forecast <- function(lat, lon, api_key = NULL) {
library(httr)
library(jsonlite)
# Construct API URL (using National Weather Service API)
# This doesn't require an API key
url <- sprintf("https://api.weather.gov/points/%.4f,%.4f", lat, lon)
# Get forecast endpoint
response <- GET(url)
if(status_code(response) == 200) {
data <- fromJSON(content(response, "text"))
forecast_url <- data$properties$forecast
# Get forecast data
forecast_response <- GET(forecast_url)
if(status_code(forecast_response) == 200) {
forecast_data <- fromJSON(content(forecast_response, "text"))
# Extract relevant information
periods <- forecast_data$properties$periods
# Process forecast periods
weather_df <- data.frame(
time = as.POSIXct(periods$startTime),
temperature = periods$temperature,
wind_speed = NA,
wind_direction = NA,
humidity = NA,
precipitation = NA
)
# Extract wind data
for(i in 1:nrow(weather_df)) {
# Parse wind speed (e.g., "10 mph")
wind_parts <- strsplit(periods$windSpeed[i], " ")[[1]]
weather_df$wind_speed[i] <- as.numeric(wind_parts[1])
# Parse wind direction (e.g., "NW", "SSE")
weather_df$wind_direction[i] <- periods$windDirection[i]
# Convert direction to degrees
dir_map <- c("N" = 0, "NNE" = 22.5, "NE" = 45, "ENE" = 67.5,
"E" = 90, "ESE" = 112.5, "SE" = 135, "SSE" = 157.5,
"S" = 180, "SSW" = 202.5, "SW" = 225, "WSW" = 247.5,
"W" = 270, "WNW" = 292.5, "NW" = 315, "NNW" = 337.5)
weather_df$wind_direction_deg[i] <- dir_map[weather_df$wind_direction[i]]
# Extract humidity and precipitation if available
if("relativeHumidity" %in% names(periods)) {
weather_df$humidity[i] <- periods$relativeHumidity[i]
}
if("probabilityOfPrecipitation" %in% names(periods)) {
weather_df$precipitation[i] <- periods$probabilityOfPrecipitation$value[i] / 100
}
}
return(weather_df)
}
}
warning("Could not load NOAA forecast data")
return(NULL)
}
# 2. Load actual fuel data from LANDFIRE
load_landfire_data <- function(bbox, resolution = 30) {
library(raster)
library(rgdal)
# LANDFIRE REST service URL for fuel model
url <- paste0("https://lfps.usgs.gov/arcgis/services/landfire/US_200EVT/MapServer/WMSServer",
"?request=GetMap",
"&service=WMS",
"&version=1.1.1",
"&layers=US_200EVT",
"&styles=",
"&format=image/tiff",
"&transparent=true",
"&srs=EPSG:4326",
"&bbox=", bbox[1], ",", bbox[2], ",", bbox[3], ",", bbox[4],
"&width=", ceiling((bbox[3]-bbox[1])/resolution*111000),
"&height=", ceiling((bbox[4]-bbox[2])/resolution*111000))
# Download temporary file
temp_file <- tempfile(fileext = ".tif")
download.file(url, temp_file, mode = "wb", quiet = TRUE)
# Read raster
landfire_raster <- raster(temp_file)
# Load LANDFIRE fuel model lookup table
# This would typically be a local file with the LANDFIRE classifications
# For this example, we'll create a simplified version
fuel_lookup <- data.frame(
value = c(0, 101:109, 111:119, 121:129, 131:139, 141:149, 151:159, 161:169, 171:179, 181:189, 191:199),
fuel_type = c("NoData",
rep("grass", 9),
rep("shrub", 9),
rep("forest", 9),
rep("forest", 9),
rep("forest", 9),
rep("forest", 9),
rep("shrub", 9),
rep("slash", 9),
rep("slash", 9),
rep("grass", 9))
)
# Return both raster and lookup table
return(list(
raster = landfire_raster,
lookup = fuel_lookup
))
}
# 3. Integrate satellite hotspot detection
integrate_satellite_hotspots <- function(grid, date) {
library(rjson)
# Construct bounding box from grid
min_lat <- min(grid$lat)
max_lat <- max(grid$lat)
min_lon <- min(grid$lon)
max_lon <- max(grid$lon)
# Format date as required (YYYY-MM-DD)
date_str <- format(date, "%Y-%m-%d")
# Query NASA FIRMS API for fire hotspots
# In a real implementation, you would register for an API key
# This is a placeholder URL
url <- sprintf("https://firms.modaps.eosdis.nasa.gov/api/area/csv/[API_KEY]/VIIRS_SNPP_NRT/%f,%f,%f,%f/%s",
min_lon, min_lat, max_lon, max_lat, date_str)
# For demonstration, we'll create synthetic hotspot data
# In a real implementation, you would parse the actual API response
hotspots <- data.frame(
lat = runif(5, min_lat, max_lat),
lon = runif(5, min_lon, max_lon),
confidence = sample(70:100, 5, replace = TRUE)
)
# Convert hotspots to grid coordinates
hotspots$x <- NA
hotspots$y <- NA
for(i in 1:nrow(hotspots)) {
# Find closest grid cell
dists <- sqrt((grid$lat - hotspots$lat[i])^2 + (grid$lon - hotspots$lon[i])^2)
closest_idx <- which.min(dists)
hotspots$x[i] <- grid$x[closest_idx]
hotspots$y[i] <- grid$y[closest_idx]
}
# Return hotspot data
return(hotspots)
}
# 4. Create an integrated model run with real-world data
run_integrated_wildfire_model <- function(bbox, start_date, duration_hours, resolution = 30) {
# 1. Set up grid based on bounding box
grid_size <- ceiling(c(
lon_cells = (bbox[3] - bbox[1]) * 111000 / resolution,
lat_cells = (bbox[4] - bbox[2]) * 111000 / resolution
))
# Create grid
grid <- expand.grid(
x = 1:grid_size[1],
y = 1:grid_size[2]
) %>%
mutate(
lon = bbox[1] + (x - 0.5) * (bbox[3] - bbox[1]) / grid_size[1],
lat = bbox[2] + (y - 0.5) * (bbox[4] - bbox[2]) / grid_size[2],
FI = 0,
FC = 0,
elevation = 0
)
# 2. Load real elevation data
elevation_data <- get_elevation_data(bbox, resolution)
for(i in 1:nrow(grid)) {
# Find closest elevation point
grid$elevation[i] <- get_nearest_elevation(elevation_data, grid$lon[i], grid$lat[i])
}
# 3. Load LANDFIRE fuel data
landfire_data <- load_landfire_data(bbox, resolution)
# Assign fuel types to grid
for(i in 1:nrow(grid)) {
# Extract LANDFIRE value at this point
lf_value <- extract(landfire_data$raster, cbind(grid$lon[i], grid$lat[i]))
# Look up fuel type
if(!is.na(lf_value)) {
fuel_idx <- which(landfire_data$lookup$value == lf_value)
if(length(fuel_idx) > 0) {
grid$fuel_type[i] <- landfire_data$lookup$fuel_type[fuel_idx]
} else {
grid$fuel_type[i] <- "grass" # Default
}
} else {
grid$fuel_type[i] <- "grass" # Default
}
}
# 4. Load weather forecast for center of area
center_lat <- (bbox[2] + bbox[4]) / 2
center_lon <- (bbox[1] + bbox[3]) / 2
weather_forecast <- load_noaa_weather_forecast(center_lat, center_lon)
# 5. Load satellite hotspot data for initial ignition
hotspots <- integrate_satellite_hotspots(grid, start_date)
# Set initial fire at hotspot locations
for(i in 1:nrow(hotspots)) {
grid_idx <- which(grid$x == hotspots$x[i] & grid$y == hotspots$y[i])
if(length(grid_idx) > 0) {
grid$FI[grid_idx] <- hotspots$confidence[i] / 100 # Scale to 0-1
}
}
# 6. Create model configuration
config <- advanced_config
config$grid_size <- list(rows = grid_size[1], cols = grid_size[2])
config$cell_size <- resolution
config$time_steps <- duration_hours
config$dt <- 1 # 1 hour time steps
config$start_time <- as.POSIXct(start_date)
config$grid_coords <- list(
lat_min = bbox[2],
lon_min = bbox[1],
lat_step = (bbox[4] - bbox[2]) / grid_size[2],
lon_step = (bbox[3] - bbox[1]) / grid_size[1]
)
# 7. Run the simulation
results <- list()
results[[1]] <- grid
for(t in 2:(duration_hours+1)) {
cat("Simulating hour", t-1, "of", duration_hours, "\n")
# Get current time
current_time <- config$start_time + (t-2) * 3600 # in seconds
# Update weather
grid_with_weather <- update_grid_weather(results[[t-1]], weather_forecast, current_time)
# Simulate fire spread with all advanced components
results[[t]] <- simulate_fire_atmosphere_step(grid_with_weather, t-1, config)
}
return(results)
}
# Helper function to update grid weather based on forecast
update_grid_weather <- function(grid, weather_forecast, current_time) {
# Find closest forecast time
time_idx <- which.min(abs(as.numeric(weather_forecast$time - current_time)))
if(length(time_idx) > 0) {
# Get weather for this time
temp <- weather_forecast$temperature[time_idx]
wind_speed <- weather_forecast$wind_speed[time_idx]
wind_dir_deg <- weather_forecast$wind_direction_deg[time_idx]
humidity <- weather_forecast$humidity[time_idx]
precip <- weather_forecast$precipitation[time_idx]
# Convert to components
wind_x <- wind_speed * cos(wind_dir_deg * pi/180)
wind_y <- wind_speed * sin(wind_dir_deg * pi/180)
# Apply to all grid cells (simplified - in reality would have spatial variation)
grid$temperature <- temp
grid$wind_x <- wind_x
grid$wind_y <- wind_y
grid$relative_humidity <- humidity
grid$precipitation <- ifelse(is.na(precip), 0, precip)
}
return(grid)
}
# Function to get elevation data
get_elevation_data <- function(bbox, resolution) {
# In a real implementation, would query elevation service
# For this example, create synthetic terrain
# Create lat/lon grid
lon_seq <- seq(bbox[1], bbox[3], length.out = ceiling((bbox[3]-bbox[1])*111000/resolution))
lat_seq <- seq(bbox[2], bbox[4], length.out = ceiling((bbox[4]-bbox[2])*111000/resolution))
coords <- expand.grid(lon = lon_seq, lat = lat_seq)
# Generate synthetic terrain (hills and valleys)
coords$elevation <- 0
# Add some hills
for(i in 1:5) {
center_lon <- runif(1, bbox[1], bbox[3])
center_lat <- runif(1, bbox[2], bbox[4])
height <- runif(1, 100, 500)
width <- runif(1, 0.05, 0.2)
coords$elevation <- coords$elevation +
height * exp(-((coords$lon - center_lon)^2 + (coords$lat - center_lat)^2) / (2 * width^2))
}
# Add a ridgeline
ridge_lon <- runif(1, bbox[1], bbox[3])
ridge_angle <- runif(1, 0, pi)
ridge_height <- runif(1, 300, 800)
ridge_width <- runif(1, 0.01, 0.05)
for(i in 1:nrow(coords)) {
# Distance to ridge
lon_diff <- coords$lon[i] - ridge_lon
lat_diff <- coords$lat[i] - bbox[2] - (coords$lon[i] - bbox[1]) * tan(ridge_angle)
dist <- abs(lat_diff * cos(ridge_angle) - lon_diff * sin(ridge_angle))
# Add elevation from ridge
coords$elevation[i] <- coords$elevation[i] +
ridge_height * exp(-(dist^2) / (2 * ridge_width^2))
}
return(coords)
}
# Helper to get nearest elevation point
get_nearest_elevation <- function(elevation_data, lon, lat) {
dists <- sqrt((elevation_data$lon - lon)^2 + (elevation_data$lat - lat)^2)
return(elevation_data$elevation[which.min(dists)])
}4.8 Running a Full Integrated Simulation
This function runs the complete integrated SDE model simulation of a real-world scenario.
# Example usage of the integrated model for a real case study
run_case_study <- function() {
# Define area of interest (e.g., part of California)
bbox <- c(-120.2, 38.5, -119.8, 38.9) # lon1, lat1, lon2, lat2
# Define time period
start_date <- as.POSIXct("2023-07-15 12:00:00", tz = "UTC")
duration_hours <- 48 # 2 days
# Run the integrated model
results <- run_integrated_wildfire_model(bbox, start_date, duration_hours)
# Visualize results at different time points
time_points <- c(1, 12, 24, 36, 48)
for(t in time_points) {
cat("Creating visualization for hour", t, "\n")
# 2D visualization
png(paste0("wildfire_hour_", t, ".png"), width = 1200, height = 800)
visualize_advanced_fire_model(results, t)
dev.off()
# 3D visualization
p <- visualize_fire_atmosphere_3D(results, t)
htmlwidgets::saveWidget(as_widget(p), paste0("wildfire_3d_hour_", t, ".html"))
}
# Create animation
create_advanced_fire_animation(results, "california_wildfire_simulation.gif")
# Return results for further analysis
return(results)
}
# Function to analyze and summarize simulation results
analyze_wildfire_results <- function(results) {
# Extract key metrics over time
time_steps <- length(results)
metrics <- data.frame(
time = 1:time_steps,
total_burned_area = numeric(time_steps),
avg_intensity = numeric(time_steps),
max_intensity = numeric(time_steps),
perimeter_length = numeric(time_steps),
spread_rate = numeric(time_steps)
)
# Calculate metrics for each time step
for(t in 1:time_steps) {
grid <- results[[t]]
# Burned area (cells with FI > 0.1)
metrics$total_burned_area[t] <- sum(grid$FI > 0.1)
# Average and max intensity
fire_cells <- grid$FI > 0.1
if(sum(fire_cells) > 0) {
metrics$avg_intensity[t] <- mean(grid$FI[fire_cells])
metrics$max_intensity[t] <- max(grid$FI)
}
# Perimeter calculation
perimeter_cells <- 0
for(i in which(fire_cells)) {
x <- grid$x[i]
y <- grid$y[i]
# Check if cell has at least one non-burning neighbor
neighbors <- which(
grid$x >= x-1 & grid$x <= x+1 &
grid$y >= y-1 & grid$y <= y+1 &
!(grid$x == x & grid$y == y)
)
if(any(grid$FI[neighbors] <= 0.1)) {
perimeter_cells <- perimeter_cells + 1
}
}
metrics$perimeter_length[t] <- perimeter_cells
# Calculate spread rate (change in burned area)
if(t > 1) {
metrics$spread_rate[t] <- metrics$total_burned_area[t] - metrics$total_burned_area[t-1]
}
}
# Create visualizations
library(plotly)
# 1. Fire growth over time
growth_plot <- plot_ly(metrics, x = ~time) %>%
add_trace(y = ~total_burned_area, name = "Burned Area",
type = "scatter", mode = "lines", line = list(color = "orange", width = 3)) %>%
layout(
title = "Wildfire Growth Over Time",
xaxis = list(title = "Time (hours)"),
yaxis = list(title = "Burned Area (cells)")
)
# 2. Fire intensity and spread rate
intensity_plot <- plot_ly(metrics, x = ~time) %>%
add_trace(y = ~avg_intensity, name = "Avg Intensity",
type = "scatter", mode = "lines", line = list(color = "red")) %>%
add_trace(y = ~max_intensity, name = "Max Intensity",
type = "scatter", mode = "lines", line = list(color = "darkred")) %>%
add_trace(y = ~spread_rate/10, name = "Spread Rate",
type = "scatter", mode = "lines", yaxis = "y2", line = list(color = "blue")) %>%
layout(
title = "Fire Intensity and Spread Rate",
xaxis = list(title = "Time (hours)"),
yaxis = list(title = "Fire Intensity",
titlefont = list(color = "red"),
tickfont = list(color = "red")),
yaxis2 = list(title = "Spread Rate (cells/hour)",
overlaying = "y", side = "right",
titlefont = list(color = "blue"),
tickfont = list(color = "blue"))
)
# 3. Perimeter analysis
perimeter_plot <- plot_ly(metrics, x = ~time) %>%
add_trace(y = ~perimeter_length, name = "Fire Perimeter",
type = "scatter", mode = "lines", line = list(color = "purple", width = 3)) %>%
add_trace(y = ~total_burned_area/10, name = "Burned Area/10",
type = "scatter", mode = "lines", line = list(color = "orange", dash = "dot")) %>%
layout(
title = "Fire Perimeter Length Over Time",
xaxis = list(title = "Time (hours)"),
yaxis = list(title = "Perimeter Length (cells)")
)
# 4. Create figures showing fire shape characteristics
# Calculate shape metrics
metrics$shape_index <- metrics$perimeter_length / (2 * sqrt(pi * metrics$total_burned_area))
shape_plot <- plot_ly(metrics, x = ~time) %>%
add_trace(y = ~shape_index, name = "Shape Index",
type = "scatter", mode = "lines", line = list(color = "green", width = 3)) %>%
layout(
title = "Fire Shape Complexity Over Time",
xaxis = list(title = "Time (hours)"),
yaxis = list(title = "Shape Index (1.0 = circle)")
)
# Return analysis results
return(list(
metrics = metrics,
plots = list(
growth = growth_plot,
intensity = intensity_plot,
perimeter = perimeter_plot,
shape = shape_plot
)
))
}
# Function to generate a comprehensive report
generate_wildfire_report <- function(results, output_dir = "./wildfire_report") {
# Create output directory if it doesn't exist
if(!dir.exists(output_dir)) {
dir.create(output_dir, recursive = TRUE)
}
# Run analysis
analysis <- analyze_wildfire_results(results)
# Save plots
htmlwidgets::saveWidget(as_widget(analysis$plots$growth),
file.path(output_dir, "fire_growth.html"))
htmlwidgets::saveWidget(as_widget(analysis$plots$intensity),
file.path(output_dir, "fire_intensity.html"))
htmlwidgets::saveWidget(as_widget(analysis$plots$perimeter),
file.path(output_dir, "fire_perimeter.html"))
htmlwidgets::saveWidget(as_widget(analysis$plots$shape),
file.path(output_dir, "fire_shape.html"))
# Generate snapshots for key time points
time_points <- c(1, round(length(results)/4), round(length(results)/2),
round(3*length(results)/4), length(results))
for(t in time_points) {
# 2D visualization
png(file.path(output_dir, paste0("wildfire_hour_", t, ".png")),
width = 1200, height = 800)
visualize_advanced_fire_model(results, t)
dev.off()
# 3D visualization
p <- visualize_fire_atmosphere_3D(results, t)
htmlwidgets::saveWidget(as_widget(p),
file.path(output_dir, paste0("wildfire_3d_hour_", t, ".html")))
}
# Create animation
create_advanced_fire_animation(results,
file.path(output_dir, "wildfire_animation.gif"))
# Generate tabular summary
write.csv(analysis$metrics, file.path(output_dir, "metrics.csv"), row.names = FALSE)
# Generate HTML report
report_html <- file.path(output_dir, "wildfire_report.html")
# Create report content
report_content <- '
<!DOCTYPE html>
<html>
<head>
<title>Advanced Wildfire Simulation Report</title>
<style>
body { font-family: Arial, sans-serif; max-width: 1200px; margin: 0 auto; padding: 20px; }
h1, h2, h3 { color: #333; }
.section { margin-bottom: 30px; }
.flex-container { display: flex; flex-wrap: wrap; }
.flex-item { flex: 1; min-width: 45%; margin: 10px; }
img { max-width: 100%; }
table { border-collapse: collapse; width: 100%; }
th, td { border: 1px solid #ddd; padding: 8px; text-align: left; }
th { background-color: #f2f2f2; }
tr:nth-child(even) { background-color: #f9f9f9; }
</style>
</head>
<body>
<h1>Advanced Wildfire Simulation Report</h1>
<div class="section">
<h2>Simulation Overview</h2>
<p>This report presents the results of an advanced stochastic differential equation (SDE) model
for wildfire propagation, incorporating weather conditions, fuel types, spotting behavior,
and fire-atmosphere interactions.</p>
</div>
<div class="section">
<h2>Fire Progression</h2>
<p>The animation below shows the fire progression over time:</p>
<img src="wildfire_animation.gif" alt="Wildfire Animation">
</div>
<div class="section">
<h2>Key Snapshots</h2>
<div class="flex-container">
<div class="flex-item">
<h3>Initial State</h3>
<img src="wildfire_hour_1.png" alt="Initial State">
<p>3D View: <a href="wildfire_3d_hour_1.html" target="_blank">Open 3D Visualization</a></p>
</div>
<div class="flex-item">
<h3>Halfway Point</h3>
<img src="wildfire_hour_' + round(length(results)/2) + '.png" alt="Halfway Point">
<p>3D View: <a href="wildfire_3d_hour_' + round(length(results)/2) + '.html" target="_blank">Open 3D Visualization</a></p>
</div>
<div class="flex-item">
<h3>Final State</h3>
<img src="wildfire_hour_' + length(results) + '.png" alt="Final State">
<p>3D View: <a href="wildfire_3d_hour_' + length(results) + '.html" target="_blank">Open 3D Visualization</a></p>
</div>
</div>
</div>
<div class="section">
<h2>Analysis Results</h2>
<div class="flex-container">
<div class="flex-item">
<h3>Fire Growth</h3>
<iframe src="fire_growth.html" width="100%" height="400px" frameborder="0"></iframe>
</div>
<div class="flex-item">
<h3>Fire Intensity</h3>
<iframe src="fire_intensity.html" width="100%" height="400px" frameborder="0"></iframe>
</div>
<div class="flex-item">
<h3>Fire Perimeter</h3>
<iframe src="fire_perimeter.html" width="100%" height="400px" frameborder="0"></iframe>
</div>
<div class="flex-item">
<h3>Fire Shape Complexity</h3>
<iframe src="fire_shape.html" width="100%" height="400px" frameborder="0"></iframe>
</div>
</div>
</div>
<div class="section">
<h2>Summary Metrics</h2>
<p>Download full metrics: <a href="metrics.csv">metrics.csv</a></p>
<table id="metrics-table">
<tr>
<th>Time (hours)</th>
<th>Burned Area (cells)</th>
<th>Average Intensity</th>
<th>Max Intensity</th>
<th>Spread Rate (cells/hour)</th>
</tr>
<!-- Table content will be filled by JavaScript -->
</table>
</div>
<script>
// Simple script to populate the table with key metrics
fetch("metrics.csv")
.then(response => response.text())
.then(data => {
const rows = data.split("\\n");
const headers = rows[0].split(",");
// Select key time points
const totalRows = rows.length - 1; // Minus header
const keyPoints = [
1, // Start
Math.round(totalRows/4),
Math.round(totalRows/2),
Math.round(3*totalRows/4),
totalRows // End
];
const table = document.getElementById("metrics-table");
keyPoints.forEach(rowIndex => {
if (rowIndex < totalRows && rowIndex > 0) {
const cells = rows[rowIndex].split(",");
const row = table.insertRow(-1);
// Add selected columns (time, burned area, avg intensity, max intensity, spread rate)
const colIndices = [0, 1, 2, 3, 5]; // Indices from the CSV
colIndices.forEach(colIndex => {
const cell = row.insertCell(-1);
let value = cells[colIndex];
// Format numbers
if (colIndex > 0 && !isNaN(value)) {
value = parseFloat(value).toFixed(2);
}
cell.textContent = value;
});
}
});
});
</script>
</body>
</html>
'
# Write HTML report
writeLines(report_content, report_html)
cat("Wildfire simulation report generated at:", output_dir, "\n")
return(output_dir)
}
enhanced wildfire simulation
wildfire fuel type analysis
wildfire dynamic - longitudinal tracking
metrics
The enhanced wildfire simulation incorporates several key improvements that create a more dynamic and non-linear behavior throughout the entire time period.
Complex Terrain: The terrain now includes a ridge, multiple hills, and a river valley, creating more realistic topography that influences fire spread.
Dynamic Weather Patterns: Weather conditions shift multiple times throughout the simulation
- Initial southwesterly winds
- A weather front with stronger westerly winds and cooler, moister conditions
- A shift to northwesterly winds, and
- A transition to hot, dry, northerly winds.
Multiple Ignition Mechanisms:
- Initial ignition point
- Enhanced ember spotting that creates new fires across barriers
- A delayed ignition point that starts across the river at \(t=25\), and
- Potential firebreak failures in extreme conditions.
More Sophisticated Fire Spread:
- Directional spread influenced by wind direction
- Stronger upslope effects on steeper terrain
- Variable rate of spread based on fuel type and moisture, and
- Fire-atmosphere coupling that creates local weather effects.
Variable Fuel Properties:
- Moisture gradient based on elevation and distance from river
- Multiple fuel types with different characteristics, and
- Fuel consumption rates that vary with fire intensity.
As well as, Stochastic Elements:
- Random variations in fire spread and intensity
- Probabilistic spotting behavior, and
- Chance-based firebreak failures.
The final (static) diagnostic plots showcase the temporal dynamics, e.g., how fire behavior changes over time and varies by fuel type. The metrics track important aspects like spread rate, fire intensity, perimeter length, and their relationships with changing weather conditions.
4.9 Final Case Study Application
We can also construct a complete case study demonstrating all the integrated SDE model components with a simulation for the Lake Tahoe region.
# Define area for a real case study (Northern California)
wildfire_case_study <- function() {
# Area of interest (near Lake Tahoe, CA)
bbox <- c(-120.2, 38.5, -119.8, 38.9) # lon1, lat1, lon2, lat2
# Define time period for simulation
start_date <- as.POSIXct("2023-07-15 12:00:00", tz = "UTC")
duration_hours <- 72 # 3 days
cat("Setting up simulation for Lake Tahoe region...\n")
# Create custom configuration for this region
config <- list(
# Grid parameters
grid_size = list(rows = 40, cols = 40),
cell_size = 100, # meters
grid_coords = list(
lat_min = bbox[2],
lon_min = bbox[1],
lat_step = (bbox[4] - bbox[2]) / 40,
lon_step = (bbox[3] - bbox[1]) / 40
),
# Fire dynamics parameters
alphaP = 0.15, # Base fire spread rate
betaP = 0.2, # Fuel influence factor
gammaP = 0.25, # Slope influence factor
deltaP = 0.3, # Heat influence
epsilonP = 0.12, # Fuel consumption rate
# Weather effect parameters
gamma_wind = 0.2, # Wind influence
delta_temp = 0.05, # Temperature influence
rh_effect = 0.8, # Humidity influence
rain_suppression = 0.5, # Rain suppression effect
temp_effect = 0.03, # Temperature effect on fuel
humidity_effect = 0.7, # Humidity effect on fuel
precipitation_effect = 0.9, # Precipitation effect
# Stochastic parameters
sigma_FI = 0.04, # Fire intensity noise
sigma_FC = 0.06, # Fuel consumption noise
# Simulation parameters
time_steps = duration_hours,
dt = 1, # 1 hour time steps
start_time = start_date,
# Ember spotting parameters
spotting_params = list(
ember_prob_threshold = 0.6,
max_spotting_distance = 2000
),
# Fire-atmosphere parameters
atmosphere_params = list(
plume_height_coef = 100,
air_temp_rise_coef = 5,
convection_strength_coef = 0.2
),
# Firefighting parameters
firefighting_start = 12, # Start firefighting after 12 hours
firefighting_params = list(
strategy = "perimeter",
resources = 10,
effectiveness = 0.3
)
)
# Generate synthetic terrain for the study area
cat("Generating terrain and fuel data...\n")
terrain_data <- get_elevation_data(bbox, 100)
# Create initial grid
grid <- expand.grid(
x = 1:config$grid_size$rows,
y = 1:config$grid_size$cols
) %>%
mutate(
lon = bbox[1] + (x - 0.5) * config$grid_coords$lon_step,
lat = bbox[2] + (y - 0.5) * config$grid_coords$lat_step,
FI = 0,
FC = 0,
elevation = 0
)
# Assign elevation
for(i in 1:nrow(grid)) {
grid$elevation[i] <- get_nearest_elevation(terrain_data, grid$lon[i], grid$lat[i]) / 1000 # Scale to 0-1
}
# Assign fuel types - simplified version for demonstration
set.seed(123) # For reproducibility
# Create realistic fuel patterns
# Forest in higher elevations, grass in lower areas, shrubs in between
grid$fuel_type <- "grass" # Default
for(i in 1:nrow(grid)) {
if(grid$elevation[i] > 0.7) {
grid$fuel_type[i] <- "forest"
} else if(grid$elevation[i] > 0.3) {
grid$fuel_type[i] <- "shrub"
}
# Add some randomness to make it more realistic
if(runif(1) < 0.2) {
if(grid$fuel_type[i] == "grass") {
grid$fuel_type[i] <- "shrub"
} else if(grid$fuel_type[i] == "shrub") {
if(runif(1) < 0.5) {
grid$fuel_type[i] <- "forest"
} else {
grid$fuel_type[i] <- "grass"
}
}
}
}
# Add fuel properties based on fuel type
fuel_types <- define_fuel_types()
for(i in 1:nrow(grid)) {
ft <- fuel_types[[grid$fuel_type[i]]]
grid$fuel_load[i] <- ft$load
grid$moisture_drying[i] <- ft$moisture_drying
grid$moisture_threshold[i] <- ft$moisture_threshold
grid$heat_content[i] <- ft$heat_content
grid$spread_rate_coef[i] <- ft$spread_rate_coef
grid$intensity_coef[i] <- ft$intensity_coef
grid$sustainability[i] <- ft$sustainability
grid$ignition_temp[i] <- ft$ignition_temp
# Initial fuel consumption at full capacity
grid$FC[i] <- 1.0
}
# Synthetic weather data for simulation period
cat("Generating weather data for simulation period...\n")
# Create synthetic diurnal pattern with a weather front passing through
hours <- 0:(duration_hours-1)
# Base diurnal temperature pattern
temps <- 25 + 10 * sin(2*pi*(hours-10)/24) # Peak at 10 hours after start
# Weather front approaches at hour 24, passes through by hour 48
# Cooler temperatures, higher humidity, shifting winds
front_effect <- rep(0, length(hours))
front_effect[24:48] <- -5 * dnorm(1:25, mean=12, sd=6)
temps <- temps + front_effect
# Relative humidity - inverse relationship with temperature
rh <- 100 - 3 * (temps - min(temps))
rh <- pmin(pmax(rh, 20), 95) # Bound between 20-95%
# Wind patterns - rotation during front passage
wind_speed <- 3 + 2 * sin(2*pi*hours/24) # Daily cycle
wind_speed[24:48] <- wind_speed[24:48] + 3 * dnorm(1:25, mean=12, sd=6) # Stronger during front
# Wind direction shifts during front passage
wind_dir <- rep(45, length(hours)) # Start from NE
wind_dir[24:48] <- 45 + 180 * (1:25)/25 # Rotate to SW during front passage
# Precipitation - chance of rain with the front
precip <- rep(0, length(hours))
precip[30:40] <- 0.3 * dnorm(1:11, mean=5, sd=2) # Peak rain at hour 35
# Combine into weather data frame
weather_data <- data.frame(
time = start_date + hours * 3600,
temperature = temps,
relative_humidity = rh,
wind_speed = wind_speed,
wind_direction = wind_dir,
precipitation = precip
)
# Add wind components
weather_data$wind_x <- weather_data$wind_speed * cos(weather_data$wind_direction * pi/180)
weather_data$wind_y <- weather_data$wind_speed * sin(weather_data$wind_direction * pi/180)
# Set initial ignition - lightning strike in forest area
cat("Setting initial ignition conditions...\n")
forest_cells <- which(grid$fuel_type == "forest")
ignition_idx <- sample(forest_cells, 1)
grid$FI[ignition_idx] <- 0.7 # Strong initial ignition
# Print ignition location
cat("Ignition at:", grid$x[ignition_idx], grid$y[ignition_idx], "\n")
cat("Fuel type at ignition:", grid$fuel_type[ignition_idx], "\n")
# Run the simulation
cat("Running wildfire simulation...\n")
results <- list()
results[[1]] <- grid
for(t in 2:(duration_hours+1)) {
cat("Simulating hour", t-1, "of", duration_hours, "\n")
# Get current time
current_time <- start_date + (t-2) * 3600 # in seconds
# Get weather for this time step
current_weather <- weather_data[t-1,]
# Update grid with current weather
grid_with_weather <- results[[t-1]]
grid_with_weather$temperature <- current_weather$temperature
grid_with_weather$relative_humidity <- current_weather$relative_humidity
grid_with_weather$wind_x <- current_weather$wind_x
grid_with_weather$wind_y <- current_weather$wind_y
grid_with_weather$precipitation <- current_weather$precipitation
# Simulate fire spread with all advanced components
next_grid <- simulate_fire_atmosphere_step(grid_with_weather, t-1, config)
# Apply firefighting if active
if(t-1 >= config$firefighting_start) {
next_grid <- apply_firefighting(next_grid, config$firefighting_params)
}
results[[t]] <- next_grid
}
# Generate comprehensive report
cat("Generating wildfire analysis report...\n")
report_path <- generate_wildfire_report(results, "./tahoe_wildfire_report")
cat("\nCase study complete! Report available at:", report_path, "\n")
return(list(
config = config,
results = results,
weather = weather_data,
report_path = report_path
))
}The advanced SDE model integrates detailed weather forecasts with spatial and temporal variation, including temperature, humidity, wind fields, and precipitation effects on fire behavior. Different vegetation types (grass, shrub, forest, slash) are now represented with unique properties that affect fire spread and intensity. Long-distance ember transport is modeled, allowing for realistic jumping of fire across barriers and accelerated spread. The model also captures how fires modify their local weather, creating convection currents and altering humidity and temperature.
The integrated model depicts how different physical processes interact to create complex wildfire behavior that matches observations from real fires. The SDE framework allows for appropriate representation of uncertainty in fire spread, reflecting the chaotic nature of wildfire dynamics.
5 Summary
The stochastic differential equation model supports deeper understanding and capacity for prediction of wildfire propagation. By combining deterministic fire physics with stochastic processes, the model captures the complex and often unpredictable nature of wildfire behavior.
This approach extends the SOCR SDE cancer treatment model, adapting it to wildfire modeling. The integration of spatial components, environmental factors, and multi-core optimization techniques enhances both the accuracy and computational efficiency of the simulations. This SDE model can serve as a foundation for developing wildfire risk assessment, firefighting resource allocation, and evacuation planning, ultimately contributing to improved wildfire management and public safety.