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.

- A burning cell turns into an empty cell
- A cell occupied by a tree becomes a burning cell if any of its eight neighbouring cells are burning
- 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
- 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.

*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 # 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: 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.001 # 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 = fig.add_subplot(111) 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) plt.show()

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

Share on Twitter Share on Facebook
## Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

## Stafford Baines 3 years, 11 months ago

Looks interesting. I'll have to peruse it a bit more to see what makes it tick.

Link | ReplyCoincidently 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.

## christian 3 years, 11 months ago

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

Link | Reply## Adam Wolniakowski 2 years, 11 months ago

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

Link | ReplyHe used fish vs. sharks populations.

## christian 2 years, 11 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.

Link | Reply## christian 2 years, 11 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.

Link | Reply## Adam Wolniakowski 2 years, 11 months ago

That's really cool!

Link | ReplyIf 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)

## christian 2 years, 11 months ago

I remember Wireworld from some reading I did a long time ago (Martin Gardner?) – your project sounds interesting (and hard...)

Link | ReplyI 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.

## Pedro Henrique Duque 7 months, 3 weeks 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?

Link | Reply## christian 7 months, 2 weeks 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.

Link | Reply## New Comment