The Forest-fire model

(9 comments)

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.

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.

An animation of the forest fire model

Current rating: 3.4

Comments

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

Stafford Baines 3 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.

Link | Reply
Current rating: 5

christian 3 years, 6 months ago

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

Link | Reply
Currently unrated

Adam Wolniakowski 2 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.

Link | Reply
Currently unrated

christian 2 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.

Link | Reply
Currently unrated

christian 2 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.

Link | Reply
Currently unrated

Adam Wolniakowski 2 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)

Link | Reply
Currently unrated

christian 2 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.

Link | Reply
Currently unrated

Pedro Henrique Duque 2 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
Currently unrated

christian 2 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
Currently unrated

New Comment

required

required (not published)

optional

required