# The Forest-fire model

A simple model of a forest fire is defined as a two-dimensional cellular automaton on a grid of cells which take one of three states: empty, occupied by a tree, or burning. The automaton evolves according to the following rules which are executed simultaneously for every cell at a given generation.

1. A burning cell turns into an empty cell
2. A cell occupied by a tree becomes a burning cell if any of its eight neighbouring cells are burning
3. A cell occupied by a tree becomes burning with a probabilty $f$ (even if none of its neighbouring cells are burning), as though struck by lightning
4. An empty cell becomes occupied by a tree with probability $p$. The model is interesting as a simple dynamical system displaying self-organized criticality. The following Python code simulates a forest fire using this model for probabilities $p = 0.05$, $f = 0.001$. A Matplotlib animation is used for visualization.

Update (January 2020): As noted by Casey (comment, below), the simulation isn't very realistic for small $f$: fires in a dense forest tend to expand out in a square pattern because diagonally-adjacent trees catch fire as readily as orthogonally-adjacent trees. This can be improved by assigning a probability of less than 1 for these trees to ignite. It isn't obvious what probability to assign, but I use 0.573 below. This is the ratio $B/A$, where $B$ is the area of overlap of a circle of radius $\frac{3}{2}$ at (0,0) (the tree on fire) with one of radius $\frac{1}{2}$ at $(0,\sqrt{2})$ (blue, the diagonally-adjacent tree), and $A$ is the overlap of this first circle with one of radius $\frac{1}{2}$ at $(0,1)$ (red, the orthogonally-adjacent tree). $A$ is equal to the "area" of this adjacent tree, since its circle is contained in the first. This code is also available on my github page.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib import colors

# Create a forest fire animation based on a simple cellular automaton model.
# The maths behind this code is described in the scipython blog article
# at https://scipython.com/blog/the-forest-fire-model/
# Christian Hill, January 2016.
# Updated January 2020.

# Displacements from a cell to its eight nearest neighbours
neighbourhood = ((-1,-1), (-1,0), (-1,1), (0,-1), (0, 1), (1,-1), (1,0), (1,1))
EMPTY, TREE, FIRE = 0, 1, 2
# Colours for visualization: brown for EMPTY, dark green for TREE and orange
# for FIRE. Note that for the colormap to work, this list and the bounds list
# must be one larger than the number of different values in the array.
colors_list = [(0.2,0,0), (0,0.5,0), (1,0,0), 'orange']
cmap = colors.ListedColormap(colors_list)
bounds = [0,1,2,3]
norm = colors.BoundaryNorm(bounds, cmap.N)

def iterate(X):
"""Iterate the forest according to the forest-fire rules."""

# The boundary of the forest is always empty, so only consider cells
# indexed from 1 to nx-2, 1 to ny-2
X1 = np.zeros((ny, nx))
for ix in range(1,nx-1):
for iy in range(1,ny-1):
if X[iy,ix] == EMPTY and np.random.random() <= p:
X1[iy,ix] = TREE
if X[iy,ix] == TREE:
X1[iy,ix] = TREE
for dx,dy in neighbourhood:
# The diagonally-adjacent trees are further away, so
# only catch fire with a reduced probability:
if abs(dx) == abs(dy) and np.random.random() < 0.573:
continue
if X[iy+dy,ix+dx] == FIRE:
X1[iy,ix] = FIRE
break
else:
if np.random.random() <= f:
X1[iy,ix] = FIRE
return X1

# The initial fraction of the forest occupied by trees.
forest_fraction = 0.2
# Probability of new tree growth per empty cell, and of lightning strike.
p, f = 0.05, 0.0001
# Forest size (number of cells in x and y directions).
nx, ny = 100, 100
# Initialize the forest grid.
X  = np.zeros((ny, nx))
X[1:ny-1, 1:nx-1] = np.random.randint(0, 2, size=(ny-2, nx-2))
X[1:ny-1, 1:nx-1] = np.random.random(size=(ny-2, nx-2)) < forest_fraction

fig = plt.figure(figsize=(25/3, 6.25))
ax.set_axis_off()
im = ax.imshow(X, cmap=cmap, norm=norm)#, interpolation='nearest')

# The animation function: called to produce a frame for each generation.
def animate(i):
im.set_data(animate.X)
animate.X = iterate(animate.X)
# Bind our grid to the identifier X in the animate function's namespace.
animate.X = X

# Interval between frames (ms).
interval = 100
anim = animation.FuncAnimation(fig, animate, interval=interval, frames=200)
plt.show()


Here is a short animation of the forest fire created by this program.

Current rating: 3.5 #### Stafford Baines 5 years, 6 months ago

Looks interesting. I'll have to peruse it a bit more to see what makes it tick.
Coincidently I'm working on a vaguely similar idea which I call 'deer park'. N deer are introduced onto a park and move move around the park randomly. The brown deer
eat the green grass which turns black and remains black for n steps after which it 'regrows'. Each deer has a weight index W which can increase to a maximum M if grass is available and decreases otherwise. If W reaches zero the deer is no more. The idea is to see if by varying the various parameters a stable deer population can found. Later additions may include predators and possibly the affect of weather - droughts slowing the rate of growth of the grass.

Current rating: 5 #### christian 5 years, 6 months ago

Sounds like an interesting project. Let me know what you find!

Currently unrated #### Adam Wolniakowski 4 years, 6 months ago

Reminds me of old Wa-Tor simulator designed by A. K. Dewdney: https://en.wikipedia.org/wiki/Wa-Tor

He used fish vs. sharks populations.

Currently unrated #### christian 4 years, 6 months ago

That's really interesting – I had not seen it before. I might see if I can adapt this code to model the Wa-Tor system.

Currently unrated #### christian 4 years, 6 months ago

Here is my implementation of the Wa-Tor population dynamics model: http://scipython.com/blog/wa-tor-world/ Thanks for the tip-off about it.

Current rating: 5 #### Adam Wolniakowski 4 years, 6 months ago

That's really cool!
If you're into cellular automata, you could also take a look into Wireworld: https://en.wikipedia.org/wiki/Wireworld

I am currently trying to make an interactive compiler to work with the Wireworld computer (https://www.quinapalus.com/wires11.html)

Current rating: 5 #### christian 4 years, 6 months ago

I remember Wireworld from some reading I did a long time ago (Martin Gardner?) – your project sounds interesting (and hard...)
I like your JavaScript programs on the Forest Fire model and Langton's Ant as well, by the way. When I get the time, I might try to implement something in Python for this site as well.

Currently unrated #### Pedro Henrique Duque 2 years, 2 months ago

Is the code complete? I am trying to run it, but I only get the initial image, without fire. Where do you use the probabilities of new tree growth and lightning strike?

Currently unrated #### christian 2 years, 2 months ago

The code runs OK under my environment (Python 3.5 / Anaconda). The probabilities p and f are used in the function iterate() which picks them up in global scope. If you let me know what Python version and libraries you have installed, I might be able to help.

Currently unrated #### Sara 3 months, 1 week ago

I have the same problem , I only get the initial image and I work in python 3 / JupyterLab

Currently unrated #### christian 3 months, 1 week ago

I can reproduce this in Jupyter Notebook. I do see the animation if I use

%matplotlib notebook

Currently unrated #### Casey 1 year, 3 months ago

I'm interested in this problem for a project I'm working on. I'm primarily interested in finding an extremely fast way to run this algorithm once to completion for a pre-defined "forest" with pre-defined fire locations and no random generation of trees or lightning strikes. I would be really interested to hear your thoughts on this. And cool code, thanks for sharing!

Current rating: 5 #### christian 1 year, 3 months ago

Hi Casey,
You're right, of course: the orthogonally- and diagonally-adjacent cells are not treated quite correctly. I guess a fix might be to assign a probability < 1 of the fire spreading to a diagonal cell when a given tree is on fire. For example, some factor related to the ratio A/B of A=the overlap of a circle of radius 2r at (0,0) with one of radius r at (0,r) and B=the overlap of this first circle with one of radius r at (r, 1.414r).

In any case, I'm not sure how quantitative this model can really be and how applicable it is to real life (there's no wind, terrain, all the trees are the same, etc.); but I'd be interested to hear how you get on with your project.
Cheers, Christian