开启掘金成长之旅！这是我参与「掘金日新计划 · 12 月更文挑战」的第1天，点击查看活动详情
本文介绍了如何使用Python的seaborn和matplotlib库绘制热力图，以仿真草原火灾演化的过程。
背景描述
-
假设你现在经营一块东西方向宽度500单位，南北方向长度1000单位的麦田。
-
初始时刻，在麦田的最北端，随机产生了5处起火点。
-
我们描述火情有以下几种情况：无火→起火→大火→燃尽。各种情况之间可以互相转化，其转化规则为：
-
无火→起火：在“大火”格的下风处，每分钟有80%的概率转化，在“大火”格的风向90度方向位置，每分钟有30%的概率转化，在“大火”格的背风处，每分钟有1%的概率转化。
-
起火→大火：某一单元格“起火”后2分钟，自动变为“大火”。
-
大火→燃尽：某一单元格转化为“大火”后10分钟，自动变为“燃尽”。
-
燃尽的格子不会转化为任何其他格子。
-
-
初始风向为北风（自北向南刮），一定时间后转化为东风，转化的时间为7分钟。
模拟的要求
产生随机的五个“起火”格，然后模拟火情传播的情况，以分钟为单位进行计算。并给出起火后30分钟，60分钟，100分钟的农田火情图，并计算最终完全烧尽需要多长时间（完全烧尽指的是农田中不再有“起火”和“大火”格）。
核心问题 如何绘制热力图
由参考链接，不难得知绘制热力图的函数为sns.heatmap
vmin和vmax参数控制热力图表示的最小值和最大值范围cmap控制做图的颜色cbar控制是否显示图例
相关代码如下
import matplotlib.pyplot as plt
import seaborn as sns
f, ax = plt.subplots(figsize=(5, 10))
sns.heatmap(farmMap, vmin=100, vmax=400, cmap='YlOrRd', cbar=True)
plt.xticks(fontsize=0)
plt.yticks(fontsize=0)
ax.set_title('Field when t=' + timeTile + 'min')
plt.show()
程序设计思路及完整代码
用二维数组表示整个农田的所有单元格，用数组中元素的值代表单元格的状态，根据规则编写状态转换的语句，利用循环来更新每一分钟过去后，新的农田状态。
农田的状态按以下规则表示，这种方法可以满足起火和大火状态时间在1-100秒内波动，算是本方案的亮点。
- 100 无火
- 200-201 起火
- 301-310 大火
- 400 燃尽
完整的程序如下。
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
# 假设1个单元格长度为10个单位长度
sns.set()
# 初始化随机种子
np.random.seed(4)
# 定义仿真图区域
maxX = 100
maxY = 50
# 初始化地图状态
# 100 无火 200-201 起火 301-310 大火 400 燃尽
farmMap = np.ones((maxX, maxY), dtype=int) * 100
# 随机生成5个起火点
start = np.random.randint(0, maxY, 5)
print(start)
for thisStart in start:
farmMap[0, thisStart] = 200
time = 0
# 如果有状态为2xx 3xx findFire = 1
findFire = 1
# 当存在2xx 3xx起火点时
while findFire == 1:
time = time + 1
# t=31 61 101时做图
if time == 31 or time == 61 or time == 101:
f, ax = plt.subplots(figsize=(5, 10))
timeTile = str(time - 1)
sns.heatmap(farmMap, vmin=100, vmax=400, cmap='YlOrRd', cbar=True)
plt.xticks(fontsize=0)
plt.yticks(fontsize=0)
ax.set_title('Field when t=' + timeTile + 'min')
plt.show()
findFire = 0
# 遍历地图上的所有点
for x in range(0, maxX):
for y in range(0, maxY):
# 处理2xx 状态 起火 如果有起火点 findFire 设置为1
if round(farmMap[x, y] / 100) == 2:
findFire = 1
# 燃烧2秒后进行状态转换 -> 大火
if farmMap[x, y] == 202:
farmMap[x, y] = 301
else:
# 不足2秒累计时间
farmMap[x, y] += 1
continue
# 处理3xx 状态 大火 如果有大火点 findFire 设置为1
if round(farmMap[x, y] / 100) == 3:
findFire = 1
# 燃烧10秒后进行状态转换 ->燃尽
if farmMap[x, y] == 310:
farmMap[x, y] = 400
continue
else:
# 对两种传播状态进行处理 t<7 纵向传播为主 北风
if time < 7:
# 对于每一种状态，生成随机种子，判断概率是否满足将周围区域由100状态转为200状态
p = np.random.random(1)
if x + 1 < maxX:
if farmMap[x + 1, y] == 100 and p <= 0.8:
farmMap[x + 1, y] = 200
p = np.random.random(1)
if x - 1 > 0:
if farmMap[x - 1, y] == 100 and p <= 0.01 and x - 1 > 0:
farmMap[x - 1, y] = 200
p = np.random.random(1)
if y + 1 < maxY:
if farmMap[x, y + 1] == 100 and p <= 0.3 and y + 1 < maxY:
farmMap[x, y + 1] = 200
p = np.random.random(1)
if y - 1 > 0:
if farmMap[x, y - 1] == 100 and p <= 0.3 and y - 1 > 0:
farmMap[x, y - 1] = 200
# 对两种传播状态进行处理 t>7时 横向传播为主 东风
else:
p = np.random.random(1)
if y - 1 > 0:
if farmMap[x, y - 1] == 100 and p <= 0.8 and y - 1 > 0:
farmMap[x, y - 1] = 200
p = np.random.random(1)
if y + 1 < maxY:
if farmMap[x, y + 1] == 100 and p <= 0.01 and y + 1 < maxY:
farmMap[x, y + 1] = 200
p = np.random.random(1)
if x - 1 > 0:
if farmMap[x - 1, y] == 100 and p <= 0.3 and x - 1 > 0:
farmMap[x - 1, y] = 200
p = np.random.random(1)
if x + 1 < maxX:
if farmMap[x + 1, y] == 100 and p <= 0.3 and x + 1 < maxX:
farmMap[x + 1, y] = 200
farmMap[x, y] = farmMap[x, y] + 1
continue
# 如果不存在 2xx 3xx 点 程序终止
if findFire == 0:
print('完全烧尽需要'+str(time) + '分钟')
break
效果如下。
t=30min时火场热力图。
t=99min时火场热力图
t=299min时火场热力图