Simulating Wildfires in a Forest🔥
Using only numpy and matplotlib🐍
Simulating any natural process is always quite a challenging task, especially when we are talking about weather and climate. In essence, a so called process-based model always relies on some set of assumptions allowing relatively simple description of a certain environment. These assumptions are a necessity since we can’t account for every process or a trajectory of each object/molecule/particle. So under the hood of this kind of models you’ll always find real-world knowledge wrapped into a set of equations/rules which approximate the environment with some accuracy.
I’m not an expert in the field of process-based modeling, however recently I stumbled upon a really captivating article of Christian Hill. I believe that it’s an extremely good and accessible introduction to real-world simulation, which could be understood by anyone holding basic math knowledge. So today I’ll try to give my perspective on the concept described in the aforementioned blog post and build a little bit more complex and realistic model.
As always the code of this article you can find on my GitHub. The code you’ll see below is fully written by me, I adopted only the Christian’s approach.
Cellular automaton
To build a wildfire simulation we’re going to rely on the concept of Cellular Automaton (CA), which is a computational model used to simulate complex systems and processes, and it consists of an N-dimensional grid of cells.
- Each cell has a state and a neighborhood (i.e. all the cells adjacent to that cell).
- The state of each cell changes over discrete time steps according to a set of rules.
- The state of each cell in the next time step is determined by the current state of the cell and the states of its neighboring cells.
- A set of rules defines how the state of a cell is updated.
Right now it may look too abstract, so let’s try to describe the wildfire problem as a CA. To do that I’m going to use a 2D grid. Each cell in this grid will take one of the following three values: tree (1), wildfire (2), burned area (0). Fig. 1 illustrates a randomly generated 3x3 grid filled with cells of all three states.
It’s time to define rules for our little simulation:
- If a cell with a tree has a wildfire in at least one of the neighboring cells, it starts burning;
- If a cell has a tree and there is no neighbor burning, the cell will start burning with an arbitrary probability f;
- If a cell is empty, it might grow a tree with an arbitrary probability p.
I hope that Fig. 1 and 2 has persuaded you that this model is super intuitive and simple, so let’s try to code it using only numpy and matplotlib.
We will need to import the following for this project:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib import colors
from matplotlib.patches import Rectangle
from datetime import timedelta,datetimeSince we work in python, let’s build the CA using object-oriented programming, i.e. create a class. For now our class will be very simple having only three function: initializing the landscape (grid), checking if neighbors have a fire and propagating the fire:
class WFSim:
def __init__(self, f=0.01, p=1e-4, h=16, w=16):
self.f = f #fire proba
self.p = p #tree growing proba
self.h = h #grid height
self.w = w #grid width
self.landscape = np.random.randint(1,2,(self.h,self.w)) #generating grid
def fire_neighbors_check(self, idx, jdx):
check = False
offsets = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
#offsets = coordinates of all the neighbors of a cell (0,0)
for di,dj in offsets: #checking if any of the neighbors have fire
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
if self.landscape[ni, nj] == 2:
check +=True
return check
def step(self, step):
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = 0 #ash after fire
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #growing a tree
if (self.f > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #firing a tree
self.landscape = new_landscape.copy() Now let’s animate the results using matplotlib:
def update(frame):
im.set_data(Sim.landscape)
ax.axis('off')
Sim.step(frame)
return [im]
colors_list = ['black', 'forestgreen', 'orange'] #setting cell colors
cmap = colors.ListedColormap(colors_list)
bounds = [0,1,2,3]
Sim = WFSim(h=64, w=64) #initializing the simulation
fig, ax = plt.subplots(figsize=(16,16))
im = ax.imshow(Sim.landscape, cmap=cmap, norm=norm)
plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0,
hspace = 0, wspace = 0)
ani = animation.FuncAnimation(fig, update, frames=80, interval=20) #frames - the number of steps in the simulation
ani.save('simple.gif', fps=1.5, savefig_kwargs={'pad_inches':0})
plt.show()
Looks beautiful, right? But the problem is this simulation is really far from reality, let’s COMPLICATE😈 it!
More cell states
The first thing I want to do is to add three new cell states: water, bedrock, grass. For the sake of convenience I’ll adjust the notation: water (-3), bedrock (-2), burned (-1), grass (0), tree (1), fire (2). All the three newly added classes will be initialized right from the beginning, and further during simulation water and bedrock cells won’t be able to take any other values (so they’re fixed). The grass cells will obey the following rules:
- a grass cell can’t burn;
- a burned cell can transform only to grass with the probability grass;
- a grass cell can transform to a tree cell with the (old) probability p (before that we transformed a burned cell to a tree);
Thus I’ve tried to make the model more realistic. In nature, after a wildfire the first thing to grow is grass, only then trees, and also water and bedrock are non-burning materials.
So let’s try to embed these new adjustments to the code:
class WFSim:
def __init__(self, f=0.01, p=1e-4, bedrock=0.005, water=0.05, h=16, w=16):
self.f = f #fire proba
self.p = p #tree growing proba
self.h = h #grid height
self.w = w #grid width
self.bedrock = bedrock #bedrock proba
self.water = water #water proba
self.landscape = np.random.randint(0,2,(self.h,self.w)) #0 - grass, 1 - tree
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
coef = 5 if self.surf_neighbors_check(i, j, "B") else 1
if self.bedrock*coef>np.random.random():
self.landscape[i, j] = -2
coef = 10 if self.surf_neighbors_check(i, j, "W") else 0.1
if self.water*coef>np.random.random():
self.landscape[i, j] = -3
def fire_neighbors_check(self, idx, jdx):
check = False
offsets = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
#offsets = coordinates of all the neighbors of a cell (0,0)
for di,dj in offsets: #checking if any of the neighbors have fire
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
if self.landscape[ni, nj] == 2:
check +=True
return check
def step(self, step):
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = 0 #ash after fire
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #growing a tree
if (self.f > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #firing a tree
self.landscape = new_landscape.copy()As you can see the process of the grid initialization has changed. On the first step we randomly create a grid comprising only grass and trees (we don’t really care about the grass/tree ratio). Then iterating through all the cells with certain probability we assign to each cell a value of water or bedrock. But in nature rocks and water basins often exist in clusters. So to encourage the simulation to have big chunks of water/bedrocks I introduced coefficients. If a cell has a water/bedrock neighbor, then the probability of it having water/bedrock increases by 10 or 5 times respectively.
colors_list = ['steelblue', 'grey','black', 'olivedrab', 'forestgreen', 'orange']
cmap = colors.ListedColormap(colors_list)
bounds = [-3, -2, -1, 0,1,2,3]
norm = colors.BoundaryNorm(bounds, cmap.N)
Sim = WFSim(h=32, w=32)
plt.imshow(Sim.landscape, cmap=cmap, norm=norm)
plt.axis(False)
plt.savefig('Landscape.png', pad_inches=0, bbox_inches='tight')
plt.show()
Let’s run it:
After adding new cell states I think we can move further and extend our external rules for this simulation. Surely, I’m talking about the weather. Wildfire intensity is a function of air temperature, wind and precipitation. So let’s try to implement these three features step by step and start with the simplest one — temperature.
Temperature
First of all, temperature depends a lot on the time component. So from now on I will talk about simulation steps as hours, and the first step is 00 h local time. To account for temperature we need to synthetically generate it. It’s not a tricky task to do since the temperature time series for one day has the shape of a bell (with the maximum around noon). Then we can use the sine function:
average_temp = 20
amplitude = 5
noise_level = 2
hours_in_day = 24
hours = np.arange(hours_in_day)
temperatures = average_temp + amplitude * np.sin(2 * np.pi * hours / 24 - np.pi / 2)
temperatures += np.random.normal(0, noise_level, 24) #noise
plt.plot(hours, temperatures, label='Simulated Temperature')
plt.xlabel('Hour of the Day')
plt.ylabel('Temperature (°C)')
plt.title('Simulated Average Daily Hourly Temperature')
plt.xticks(np.arange(0, 25, 3))
plt.legend()
plt.grid(True)
plt.savefig('temperature.png')
plt.show()
So each time we run this code, we’ll be able to get 24 consequent values of air temperature. How do we use them in the simulation? Easy! For example, let’s suggest that if a temperature at a certain step exceeds 25°C, the probability f (a random tree starting burning) will double. So we will replace previous code
def step(self, step):
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = -1 #ash after fire
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #growing a tree
if (self.f > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #firing a tree
self.landscape = new_landscape.copy()with
def step(self, step):
if step%24==0 and step>0:
self.temp = self.temperature()
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = -1 #ash after fire
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #growing a tree
coef = 2 if self.temp[step%24]>25 else 1
if (self.f*coef > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #firing a tree
self.landscape = new_landscape.copy()
Wind
Adding wind to this simulation will take a little bit more effort. In this work, wind will only have a direction component (no speed). The possible wind directions will be S, N, W, E, NE, NW, SE, SW. In general the wind rule sounds like this: fire propagates to the cells along the wind (Fig. 8).
To implement it, we’ll need to edit the offset rule a little:
#....................
self.wind = np.random.choice(['calm', 'S', 'N', 'W', 'E', 'SW', 'SE', 'NW', 'NE'])
self.offsets = {'calm': [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)],
'N': [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1)],
'S': [(0,-1),(0,1),(1,-1),(1,0),(1,1)],
'E': [(-1,0),(-1,1),(0,1),(1,0),(1,1)],
'W': [(-1,-1),(-1,0),(0,-1),(1,-1),(1,0)],
'NE': [(-1,0),(-1,1),(0,1)],
'NW': [(-1,-1),(-1,0),(0,-1)],
'SE': [(0,1),(1,0),(1,1)],
'SW': [(0,-1),(1,-1),(1,0)]}
def fire_neighbors_check(self, idx, jdx):
check = False
offsets = self.offsets[self.wind]
for di,dj in offsets: #checking if any of the neighbors have fire
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
if self.landscape[ni, nj] == 2:
check +=True
return checkNow you can select the wind direction manually or randomly and get something like this:
Precipitation
What’s the source of precipitation? Right, clouds. So let’s introduce another cell type to the model: cloud (3). The rule imposed by the clouds will be the following:
If a cloud covers a burned or burning cell, then at the next step this cell will transform to grass.
Now we have two difficult things to do:
- Realistically generate clouds;
- Move the clouds;
To generate the cloud let’s introduce a new parameter cloud, which is a probability of a cloud being generated at the current step. To generate a cloud we will create a new function inside the class. The idea will be the following:
- Randomly pick a cell on the landscape and assign it the cloud state (3);
- Randomly pick one neighboring cell and assign it the cloud state as well;
- Randomly pick one neighboring cell of the previously picked cell and assign it the cloud state as well;
- Repeat this N times (N is a cloud size);
The algorithm allows creating quite realistic clouds, which have complex shapes (not simple squares or rectangles):
def generate_cloud(self):
size = 16 #arbitrary cloud size
mask = np.zeros((self.h, self.w))
offsets = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
idx_prev, jdx_prev = np.random.randint(0,self.h),np.random.randint(0,self.w)
for _ in range(size):
for di,dj in [offsets[np.random.randint(len(offsets))]]:
ni, nj = idx_prev+di, jdx_prev+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
cell = (np.random.randint(min(0,idx_prev-1),min(self.h,idx_prev+1)),
np.random.randint(min(0,jdx_prev-1),min(self.w,jdx_prev+1)))
mask[ni, nj] = 1
idx_prev, jdx_prev = ni, nj
return mask.astype(bool)To move the cloud we will need to introduce a new variable storing the landscape at the previous stage. Why? Imagine a cloud covering a water cell, a tree cell and a burned cell. When it moves, we need to keep the water and the tree cells the same, but change the burned one to the grass state. So keeping the previous state is inevitable.
Our clouds will be moving along the wind as well, so at each iteration a cloud moves with the pace of one cell along the wind.
If there is no wind, the clouds degrade and disappear after one step.
Let’s code it!
#..........
self.cloud_offsets = {'calm': [],
'N': [(1,0)],
'S': [(-1,0)],
'E': [(0,-1)],
'W': [(0,1)],
'NE': [(1,-1)],
'NW': [(1,1)],
'SE': [(-1,-1)],
'SW': [(-1,1)]}
def cloud_move(self):
offsets = self.cloud_offsets[self.wind]
mask = np.zeros((self.h, self.w))
for idx in range(self.landscape.shape[0]):
for jdx in range(self.landscape.shape[1]):
for di,dj in offsets:
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w and self.landscape[idx, jdx]==3:
mask[ni,nj] = 1
return mask.astype(bool)Now we can add cloud generating and movement to the step function:
def step(self, step):
if step%24==0 and step>0:
self.temp = self.temperature()
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if new_landscape[i, j] == 3:
if self.old_landscape[i,j]==-1 or self.old_landscape[i,j]==2:
new_landscape[i, j] = 0 #grow grass if burned or burning
else:
new_landscape[i, j] = self.old_landscape[i,j] #else skip
if new_landscape[i, j] == 2:
new_landscape[i, j] = -1 #ash
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #grow a tree on a grass cell
coef = 2 if self.temp[step%24]>25 else 1
if (self.f*coef > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #burn
self.old_landscape = new_landscape.copy() #store the updated landscape
if 3 in self.landscape and self.wind!='calm':
new_landscape[self.cloud_move()] = 3 #move a cloud if exists
if (self.cloud > np.random.rand()):
new_landscape[self.generate_cloud()] = 3 #generate a new cloud
self.landscape = new_landscape.copy() #update the main landscape
Metrics
I really love what we have here, but there is a final detail to add — some numbers. Let’s create a little box showing temperature, wind, tree cover, burned pixel ratio and the step:
def update(frame):
im.set_data(Sim.landscape)
info_text = (
f"Date: {initial_date + timedelta(hours=frame)}\n"
f"Wind: {Sim.wind}\n"
f"Temperature: {round(Sim.temp[frame%24],1)}°C\n"
f"Burned area: {Sim.burned_ratio} %\n"
f"Tree cover: {Sim.tree_cover} %"
)
text.set_text(info_text)
ax.axis('off')
Sim.step(frame)
return [im]
Sim = WFSim(h=64, w=64)
initial_date = datetime(2012, 8, 19, 0)
fig, ax = plt.subplots(figsize=(16,16))
im = ax.imshow(Sim.landscape, cmap=cmap, norm=norm)
text = ax.text(2, 8, "", ha='left', fontsize=20, color='white',
bbox=dict(facecolor='black', alpha=0.8, edgecolor='black',boxstyle='round,pad=0.5'))
plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0,
hspace = 0, wspace = 0)
ani = animation.FuncAnimation(fig, update, frames=80, interval=20)
ani.save('final.gif', fps=1.5, savefig_kwargs={'pad_inches':0})
plt.show()
Much better now!
This simulation is already complex, but it definitely can be improved in terms of new external factors/cell states, so it’ll resemble the real world much more. I encourage you to play around with the simulation and share down below in comments what factors you’d add to the model.
Another thing to think about is what natural/social processes can be visualized using Cellular Automaton. Please, feel free to share your ideas!
I hope that this article has inspired you to try to build your own CA and brought you a new perspective on how to model a wildfire. Thanks for your time:)
===========================================
All my publications on Medium are free and open-access, that’s why I’d really appreciate if you followed me here!
P.s. I’m extremely passionate about (Geo)Data Science, ML/AI and Climate Change. So if you want to work together on some project pls contact me in LinkedIn.
🛰️Follow for more🛰️
