Как стать автором
Поиск
Написать публикацию
Обновить

Моделирование лесных пожаров: теория, клеточный автомат на Python

Уровень сложностиСредний
Время на прочтение10 мин
Количество просмотров6.9K
Python*Алгоритмы*Читальный залНаучно-популярноеЭкология
Обзор
Перевод
Всего голосов 23: ↑23 и ↓0+23
Комментарии15

Комментарии 15

Дополнение: описанная симуляция не очень реалистична для небольших пожаров – пожары в густом лесу в такой модели имеют тенденцию распространяться в форме квадрата, потому что деревья, расположенные по диагонали, загораются так же легко, как и деревья, расположенные по ортогональной оси

А если попробовать не клеточки, а замощение правильными шестиугольниками (сотами)? Кажется, там всё "равномернее"

Всего голосов 3: ↑3 и ↓0+3

Всё так) Это одна из вариаций клеточного автомата, с ней не будет надобности в этом дополнении

Всего голосов 1: ↑1 и ↓0+1

Так у вас пожар распространяется до бесконечности

Добавьте что дерево соседствующее с пожапром загорится с некоторой вероятностью

А еще лучше набор вероятностей в зависимости от от количества соседних пожаров

Всего голосов 2: ↑2 и ↓0+2

Можно считать, что мы добавили характеристику повышенной влажности. Пожалуйста)

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

# https://scipython.com/blog/the-forest-fire-model/
# Christian Hill, январь 2016.
# Обновлено в январе 2020.
# Дополнения: ветер, водные преграды by Jen Ciarochi:
# https://triplebyte.com/blog/how-fire-spreads-mathematical-models-and-simulators

NY = (-1, -1, -1, 0, 0, 1, 1, 1)
NX = (-1, 0, 1, -1, 1, -1, 0, 1)
NZ = (0.1, 0.1, 0.1, 0.1, 1, 0.1, 1, 1)
EMPTY, TREE, FIRE, WATER = 0, 1, 2, 3
colors_list = [(0.2, 0, 0), (0, 0.5, 0), (1, 0, 0), 'orange', 'blue']
cmap = colors.ListedColormap(colors_list)
bounds = [0, 1, 2, 3, 4]
norm = colors.BoundaryNorm(bounds, cmap.N)


def iterate(X):
    """Итерации леса в соответствии с четырьмя правилами лесного пожара"""

    X1 = np.zeros((ny, nx))
    for ix in range(1, nx - 1):
        for iy in range(1, ny - 1):
            if X[iy, ix] == WATER:
                X1[iy, ix] = WATER
            elif X[iy, ix] == EMPTY and np.random.random() <= p:
                X1[iy, ix] = TREE
            elif X[iy, ix] == TREE:
                X1[iy, ix] = TREE
                burning_neighbours = 0
                for dx, dy, dz in zip(NY, NX, NZ):
                    if abs(dx) == abs(dy) and np.random.random() < 0.573:
                        continue
                    if X[iy+dy, ix+dx] == FIRE and np.random.random() <= dz:
                        # X1[iy, ix] = FIRE
                        burning_neighbours += 1
                if np.random.random() <= spread_chances[burning_neighbours]:
                    X1[iy, ix] = FIRE
    return X1


forest_fraction = 0.2
p = 0.03
spread_chances = [0.0001] + [0.5 + 0.06*x for x in range(8)]
nx, ny = 100, 100
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
X[10:90, 10:15] = WATER
X[10:90, 40:45] = WATER
X[10:90, 60:65] = WATER
X[10:90, 80:85] = WATER

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 animate(i):
    im.set_data(animate.X)
    animate.X = iterate(animate.X)


animate.X = X

interval = 50
anim = animation.FuncAnimation(fig, animate, interval=interval, frames=200)
plt.show()

Жаль, тут в блоках кода нельзя вручную выделять цветом/жиром.

Изменения:
33 строка – добавили счётчик горящих соседей (учитывает 0.573 шанс на диагональных);
37-39 – теперь не объявляем дерево горящим при первом же огненном соседе с break, а просто увеличиваем счётчик;
40-41 – убираем else, он уже не нужен. По spread_chances от количества горящих соседей бросаем кубик на продолжение пожара;
45 – чуть снизили вероятность появления нового дерева;
46 – f убрали, это теперь spread_chances[0];
47 – добавили вероятность продолжения пожара, от количества горящих соседей. [0.0001, 0.5, 0.56, 0.62, 0.68, 0.74, 0.8, 0.86, 0.92]. Немного поигрался со значениями, эти выглядят неплохо.

upd: эх, при редактировании номера строк видны, а при прочтении нет :с

Всего голосов 1: ↑1 и ↓0+1

Есть одно упущение в моделях: деревья не растут примерно так же быстро, как они горят.

Всего голосов 3: ↑3 и ↓0+3

Опустим шутку про рост деревьев в реальном времени.

Модель без ветра и воды; вероятность появления нового дерева в два раза ниже вероятности первичного возгорания; оставим низкую влажность (0.8+0.03*x), при высокой совсем наглядность пропадёт; начало с 60% леса:

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

# https://scipython.com/blog/the-forest-fire-model/
# Christian Hill, январь 2016.
# Обновлено в январе 2020.

NY = (-1, -1, -1, 0, 0, 1, 1, 1)
NX = (-1, 0, 1, -1, 1, -1, 0, 1)
EMPTY, TREE, FIRE = 0, 1, 2
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):
    """Итерации леса в соответствии с четырьмя правилами лесного пожара"""

    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
            elif X[iy, ix] == TREE:
                X1[iy, ix] = TREE
                burning_neighbours = 0
                for dx, dy in zip(NY, NX):
                    if abs(dx) == abs(dy) and np.random.random() < 0.573:
                        continue
                    if X[iy+dy, ix+dx] == FIRE:
                        burning_neighbours += 1
                if np.random.random() <= spread_chances[burning_neighbours]:
                    X1[iy, ix] = FIRE
    return X1


forest_fraction = 0.6
p = 0.0001
spread_chances = [0.0002] + [0.8 + 0.03*x for x in range(8)]
nx, ny = 100, 100
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')

def animate(i):
    im.set_data(animate.X)
    animate.X = iterate(animate.X)


animate.X = X

interval = 5
anim = animation.FuncAnimation(fig, animate, interval=interval, frames=200)
plt.show()

Всего голосов 5: ↑5 и ↓0+5

Этот вариант выглядит как-то странно. Сначала всё выгорело, потом какие-то отдельно стоящие деревья, которые вспыхивают...

Имхо, сейчас в модели главный недостаток - деревья вырастают на пустых клетках "из ничего". Нет леса. Обычно они вырастают из семян/шишек/etc, которые падают на землю недалеко от имеющихся деревьев. Потому надо вероятность вырастания дерева считать по количеству деревьев в соседних (или не слишком отдаленных) зеленых клеток.

Всего голосов 1: ↑1 и ↓0+1

Окей

Добавили новые состояния клеток – семена и ростки. У всех свои шансы на сгорание (с потолка), в зависимости от количества огня вокруг. Деревья разбрасывают семена, аналогично распространению огня, но в зависимости от количества взрослых деревьев вокруг, с шансом 0.0001 появления семечка из ничего. Семечки (коричневые) за 1 шаг превращаются в ростки, ростки (светло-зелёные) – в деревья. Модель без ветра, с водой.

Выглядит хаотично, надо со значениями играться.

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

# https://scipython.com/blog/the-forest-fire-model/
# Christian Hill, январь 2016.
# Обновлено в январе 2020.
# Дополнения: ветер, водные преграды by Jen Ciarochi:
# https://triplebyte.com/blog/how-fire-spreads-mathematical-models-and-simulators

NY = (-1, -1, -1, 0, 0, 1, 1, 1)
NX = (-1, 0, 1, -1, 1, -1, 0, 1)
NZ = (0.1, 0.1, 0.1, 0.1, 1, 0.1, 1, 1)
EMPTY, TREE, FIRE, WATER, SEED, SPROUT = 0, 1, 2, 3, 4, 5
colors_list = [(0.2, 0, 0), (0, 0.5, 0), (1, 0, 0), 'blue', 'brown', 'lightgreen']
cmap = colors.ListedColormap(colors_list)
bounds = [0, 1, 2, 3, 4, 5, 6]
norm = colors.BoundaryNorm(bounds, cmap.N)


def iterate(X):
    """Итерации леса в соответствии с четырьмя правилами лесного пожара"""

    X1 = np.zeros((ny, nx))
    for ix in range(1, nx - 1):
        for iy in range(1, ny - 1):
            cell = X[iy, ix]
            if cell == WATER:
                X1[iy, ix] = WATER
                continue
            elif cell == FIRE:
                X1[iy, ix] = EMPTY
                continue
            burning_neighbours = 0
            fertile_neighbours = 0
            for dy, dx in zip(NY, NX):
                if abs(dy) == abs(dx) and np.random.random() < 0.573:
                    continue
                if X[iy+dy, ix+dx] == FIRE:
                    burning_neighbours += 1
                elif X[iy+dy, ix+dx] == TREE:
                    fertile_neighbours += 1
            if np.random.random() <= fire_spread_chances[X[iy, ix]][burning_neighbours]:
                X1[iy, ix] = FIRE
            elif cell == SEED:
                X1[iy, ix] = SPROUT
            elif cell == SPROUT:
                X1[iy, ix] = TREE
            elif cell == TREE:
                X1[iy, ix] = TREE
            elif np.random.random() <= forest_spread_chances[fertile_neighbours]:
                X1[iy, ix] = SEED
    return X1


forest_fraction = 0.3
p = 0.0001
fire_spread_chances = {
        0: [0] * 9,  # EMPTY
        1: [0.0002] + [0.8 + 0.03*x for x in range(8)],  # TREE
        4: [0.0000] + [0.2 + 0.05*x for x in range(8)],  # SEED
        5: [0.0001] + [0.3 + 0.07*x for x in range(8)],  # SPROUT
    }
forest_spread_chances = [0.0001] + [0.1 + 0.1*x for x in range(8)]
nx, ny = 100, 100
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
X[10:50, 10:12] = WATER
X[30:31, 10:30] = WATER
X[10:50, 28:30] = WATER

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 animate(i):
    im.set_data(animate.X)
    animate.X = iterate(animate.X)


animate.X = X

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

Всего голосов 3: ↑3 и ↓0+3
Обычно они вырастают из семян/шишек/etc
Спящие почки же…
Комментарий пока не оценивали0
НЛО прилетело и опубликовало эту надпись здесь

Лаадно

На выгоревшей земле (серый) 5 ходов ничего не растёт

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

# https://scipython.com/blog/the-forest-fire-model/
# Christian Hill, январь 2016.
# Обновлено в январе 2020.
# Дополнения: водные преграды by Jen Ciarochi:
# https://triplebyte.com/blog/how-fire-spreads-mathematical-models-and-simulators

NY = (-1, -1, -1, 0, 0, 1, 1, 1)
NX = (-1, 0, 1, -1, 1, -1, 0, 1)
NZ = (0.1, 0.1, 0.1, 0.1, 1, 0.1, 1, 1)
EMPTY, TREE, FIRE, WATER, SEED, SPROUT, BURNT = 0, 1, 2, 3, 4, 5, 6
colors_list = [(0.2, 0, 0), (0, 0.5, 0), (1, 0, 0), 'blue', 'brown', 'lightgreen', 'darkgray']
cmap = colors.ListedColormap(colors_list)
bounds = tuple(range(8))
norm = colors.BoundaryNorm(bounds, cmap.N)
burnt_cells = {}


def iterate(X):
    """Итерации леса в соответствии с четырьмя правилами лесного пожара"""

    X1 = np.zeros((ny, nx))
    for ix in range(1, nx - 1):
        for iy in range(1, ny - 1):
            cell = X[iy, ix]
            if cell == WATER:
                X1[iy, ix] = WATER
                continue
            elif cell == FIRE:
                burnt_cells[(iy, ix)] = 5
                X1[iy, ix] = BURNT
                continue
            elif cell == BURNT:
                if burnt_cells[(iy, ix)]:
                    burnt_cells[(iy, ix)] -= 1
                    X1[iy, ix] = BURNT
                else:
                    del burnt_cells[(iy, ix)]
                    X1[iy, ix] = EMPTY
                continue
            burning_neighbours = 0
            fertile_neighbours = 0
            for dy, dx in zip(NY, NX):
                if abs(dy) == abs(dx) and np.random.random() < 0.573:
                    continue
                if X[iy+dy, ix+dx] == FIRE:
                    burning_neighbours += 1
                elif X[iy+dy, ix+dx] == TREE:
                    fertile_neighbours += 1
            if np.random.random() <= fire_spread_chances[X[iy, ix]][burning_neighbours]:
                X1[iy, ix] = FIRE
            elif cell == SEED:
                X1[iy, ix] = SPROUT
            elif cell == SPROUT:
                X1[iy, ix] = TREE
            elif cell == TREE:
                X1[iy, ix] = TREE
            elif np.random.random() <= forest_spread_chances[fertile_neighbours]:
                X1[iy, ix] = SEED
    return X1


forest_fraction = 0.3
p = 0.0001
fire_spread_chances = {
        0: [0] * 9,  # EMPTY
        1: [0.0002] + [0.8 + 0.03*x for x in range(8)],  # TREE
        4: [0.0000] + [0.2 + 0.05*x for x in range(8)],  # SEED
        5: [0.0001] + [0.3 + 0.07*x for x in range(8)],  # SEED
    }
forest_spread_chances = [0.0001] + [0.1 + 0.1*x for x in range(8)]
nx, ny = 100, 100
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
X[10:50, 10:12] = WATER
X[30:31, 10:30] = WATER
X[10:50, 28:30] = WATER
X[11:40, 30:32] = WATER
X[10:40, 43:45] = WATER
X[25:26, 30:44] = WATER
X[10:11, 30:44] = WATER

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 animate(i):
    im.set_data(animate.X)
    animate.X = iterate(animate.X)


animate.X = X

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

Значения из прошлого комментария, с припиской "выглядит хаотично. С ними ещё нужно играться"

Комментарий пока не оценивали0

а какая версия питона и матплотлиба ?
у меня python3.9 c matplotlib 3.5.1 пишет

raise ValueError("BoundaryNorm is not invertible")

Комментарий пока не оценивали0

Python 3.11, matplotlib 3.6.3

Всего голосов 1: ↑1 и ↓0+1
НЛО прилетело и опубликовало эту надпись здесь

Спасибо за перевод этой статьи, переписал на C#, буду использовать в своей игре, как симуляцию песчаной бури.

Всего голосов 1: ↑1 и ↓0+1
Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации

Показать лучшие за всё время

Истории

Ваш аккаунт

Разделы

Информация

Услуги

Техническая поддержка
© 2006–2024, Habr