O Fire Model em questão considera uma floresta onde toda lateral esquerda está pegando fogo, e esse fogo pode se espalhar nas quatro direções cardinais, bastando que haja uma árvore em uma dessas posições.
Dessa forma, a variável chave para analisar o quanto o fogo irá se espalhar é a densidade de árvores nessa floresta.
Podemos aplicar o Fire Model usando o seguinte código:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
# Displacements from a cell to its Four nearest neighbours
neighbourhood = ((-1,0), (0,-1), (0, 1), (1,0))
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)
# The initial fraction of the forest occupied by trees.
forest_fraction = 0.58
# Forest size (number of cells in x and y directions).
nx, ny = 100, 100
# Initialize the forest grid.
X = np.zeros((ny, nx))
X[0:100, 0] = 2
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
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] == TREE:
X1[iy,ix] = TREE
for dx,dy in neighbourhood:
if X[iy+dy,ix+dx] == FIRE:
X1[iy,ix] = FIRE
break
return X1
cjm = [iterate(X)]
while (np.array_equal(iterate(cjm[len(cjm)-1]), cjm[len(cjm)-1]) == False):
ai = iterate(cjm[len(cjm)-1])
cjm.append(ai)
len(cjm)
def grd(X):
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')
grd(cjm[len(cjm)-1])
Fire model com densidade fixada em 0.58:

Nota-se que a um ponto crítico nesse modelo, a partir do qual a floresta começa a ser incendiada em proporção quase total. Acredita-se que esse ponto seja quando a densidade é igual a 0.59 ( podemos testar isso reaplicando o código alterando a variável forest_fraction):
Duas simulações com densidade em 0.57:


Duas simulações com densidade em 0.61:


Também podemos adicionar extensões a esse modelo. Um bom exemplo é adicionar uma probabilidade para o incêndio de árvores vizinhas, levando o incêndio a ser bem menor e mais focado, conforme código abaixo:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
# Displacements from a cell to its four nearest neighbours
neighbourhood = ((-1,0), (0,-1), (0, 1),(1,0))
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 grd(X):
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')
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] == TREE:
X1[iy,ix] = TREE
for dx,dy in neighbourhood:
if X[iy+dy,ix+dx] == FIRE and np.random.random() <= fs:
X1[iy,ix] = FIRE
break
return X1
# The initial fraction of the forest occupied by trees.
forest_fraction = 0.72
# Probability of fire spread.
fs = 0.75
# Forest size (number of cells in x and y directions).
nx, ny = 100, 100
# Initialize the forest grid.
X = np.zeros((ny, nx))
X[0:100, 0] = 2
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
cjm = [iterate(X)]
while (np.array_equal(iterate(cjm[len(cjm)-1]), cjm[len(cjm)-1]) == False):
ai = iterate(cjm[len(cjm)-1])
cjm.append(ai)
grd(cjm[len(cjm)-1])
Simulação do fire model com densidade em 0.72 e probabilidade de spread em 0.75:

Também podemos adicionar o vento, considerando que ele reduz a probabilidade de spread do fogo que se propaga na direção do vento e diminui a do fogo que se propaga contra o vento, protegendo mais a parte da floresta que está mais perto da direção do vento:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
# southwind force
swind = 0.25
#northwind force
nwind = -swind
#westwind force
wwind = 0.25
#eastwind force
ewind = -wwind
# Displacements from a cell to its eight nearest neighbours
nbnw = ((-1,0, wwind), (0,-1, nwind), (0, 1, swind),(1,0, ewind))
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 grd(X):
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')
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] == TREE:
X1[iy,ix] = TREE
for dx,dy,dw in nbnw:
if X[iy+dy,ix+dx] == FIRE and np.random.random() <= fs+dw :
X1[iy,ix] = FIRE
break
return X1
# The initial fraction of the forest occupied by trees.
forest_fraction = 0.82
# Probability of fire spread.
fs = 0.57
# Forest size (number of cells in x and y directions).
nx, ny = 100, 100
# Initialize the forest grid.
X = np.zeros((ny, nx))
X[0:100, 0] = 2
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
cjm = [iterate(X)]
while (np.array_equal(iterate(cjm[len(cjm)-1]), cjm[len(cjm)-1]) == False):
ai = iterate(cjm[len(cjm)-1])
cjm.append(ai)
grd(cjm[len(cjm)-1])
Simulação do fire model com vento na direção sudoeste, força de 0.25, densidade da floresta de 0.82 e probabilidade de spread de 0.57:

Logo, como o vento está vindo do sudoeste, a região mais ao sul da floresta fica preservada (não notamos o mesmo efeito na região mais a oeste devido ao incêndio ser oriundo dela). Por fim, podemos permitir ao fogo que de saltos longos por meio de faíscas na direção do vento. Isso resulta em faixas estreitas de floresta que ficam preservadas, cercadas por áreas devastadas. Saltos longos são incluídos pelo seguinte código:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
# southwind force
swind = 0.25
#northwind force
nwind = -swind
#westwind force
wwind = 0.25
#eastwind force
ewind = -wwind
# Displacements from a cell to its four nearest neighbours
nbnw = ((-1,0, wwind), (0,-1, nwind), (0, 1, swind),(1,0, ewind))
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 grd(X):
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')
ns = list()
for i in range(0, len(nbnw)):
n = (nbnw[i][0]*nbnw[i][8], nbnw[i][9]*nbnw[i][10])
ns.append(n)
amp = 20/2
ni = (ns[0][0]+ns[1][0]+ns[2][0]+ns[3][0], ns[0][11]+ns[1][12]+ns[2][13]+ns[3][14])
wd = (round(amp*ni[0]), round(amp*ni[1]))
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] == TREE:
X1[iy,ix] = TREE
for dx,dy,dw in nbnw:
if X[iy+dy,ix+dx] == FIRE and np.random.random() <= fs+dw :
X1[iy,ix] = FIRE
if iy+wd[1] <= 99 and ix+wd[0] <= 99 and X[iy+wd[1], ix+wd[0]] == FIRE and np.random.random() <= lfs:
X1[iy,ix] = FIRE
break
return X1
# The initial fraction of the forest occupied by trees.
forest_fraction = 0.72
# Probability of fire spread.
fs = 0.39
#probability of long-fire spread.
lfs = fs
# Forest size (number of cells in x and y directions).
nx, ny = 100, 100
# Initialize the forest grid.
X = np.zeros((ny, nx))
X[0:100, 0] = 2
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
cjm = [iterate(X)]
while (np.array_equal(iterate(cjm[len(cjm)-1]), cjm[len(cjm)-1]) == False):
ai = iterate(cjm[len(cjm)-1])
cjm.append(ai)
grd(cjm[len(cjm)-1])
Simulação do fire model com saltos longos e vento na direção sudoeste, força de 0.25, densidade da floresta de 0.72 e probabilidade de spread de 0.39:

Uma referência importante para esse trabalho foi o Fire Model aplicado a trovões, presente nesse link.