Untitled
unknown
plain_text
3 years ago
10 kB
5
Indexable
import random
import numpy as np
import math
import sys
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib import colors
from matplotlib.animation import FuncAnimation
EMPTY, FLAMMABLE, BURNING, BURNED = 0,1,2,3
def get_states(my_grid):
state_grid = np.empty((GRID_HEIGHT, GRID_WIDTH), dtype=int)
#state_grid = [map(myfunc, x) for x in my_grid]
for y in range(0, GRID_HEIGHT ):
for x in range(0, GRID_WIDTH ):
state_grid[y][x] = my_grid[y][x].state
return state_grid
def get_time_to_burn(my_grid):
time_to_burn_grid = np.empty((GRID_HEIGHT, GRID_WIDTH), dtype=float)
for y in range(0, GRID_HEIGHT ):
for x in range(0, GRID_WIDTH ):
time_to_burn_grid[y][x] = my_grid[y][x].time_to_burn
return time_to_burn_grid
class Fuel:
def __init__(self, Wo, Rho_p, Stau, Mx, Se, h, sigma, delta):
self.Wo = Wo
self.Rho_p = Rho_p
self.Stau = Stau
self.Mx = Mx
self.Se = Se
self.h = h
self.sigma = sigma
self.delta = delta
# w_o Oven-dry fuel load
# ρ_p Fuel particle density
# S_τ Total mineral content
# M_x Moisture of extinction
# S_e Effective mineral content
# h Heat content of the fuel
# σ Surface-area-to-volume-ratio
# δ Fuel bed depth
grass = Fuel(1.347, 512, 0.0555, 25, 0.01, 18622, 57, 0.76)
shrub = Fuel(2.35, 512, 0.0555, 26, 0.01, 18622, 55, 0.80)
conifer = Fuel(2.4, 512, 0.0555, 30, 0.01, 18622, 53, 0.72)
# M_f Fuel moisture
# U Wind speed
# θ Slope angle at midflame height
class Cell:
def __init__(self, Mf, U, pi, fuel=None):
self.Mf = Mf
self.U = U
self.pi = pi
self.fuel = fuel
self.time_to_burn = sys.maxsize #change to better
self.has_burning_neighbour = False
if self.fuel is None:
self.state = EMPTY # let's write burning, burned...
self.ROS = 0
self.time_to_burn = None
else:
self.state = FLAMMABLE
self.ROS = self.calc_ROS()
def calc_ROS(self):
Gm_max = ((30.48 * self.fuel.sigma) ** 1.5) / ( 495 + .0594 * (30.48 * self.fuel.sigma) ** 1.5) # FROM FIRESTATION =/= Paper
Rho_b = self.fuel.Wo / self.fuel.delta #kg/m
b = Rho_b / self.fuel.Rho_p # m^2
bop = 0.20395 * self.fuel.sigma ** (-0.8189) # CORRIGIDO #cm^-1
A = 8.9033 * self.fuel.sigma ** (-0.7913) # CORRIGIDO #cm^-1
Gm = Gm_max * ((b / bop) ** A) * np.exp(A * (1 - b / bop)) # FROM FIRESTATION =/= Paper
Wn = self.fuel.Wo * (1 - self.fuel.Stau) # CORRIGIDO #kg/m^2
Tau_m = self.Mf / self.fuel.Mx
Eta_m = 1 - 2.59 * Tau_m + 5.11 * (Tau_m ** 2) - 3.52 * (Tau_m ** 3)
Eta_s = 0.174 * self.fuel.Se ** (-0.19)
Ir = Gm * Wn * self.fuel.h * Eta_m * Eta_s / 60 # OUTPUT R in m/s
Xi = (192 + 7.9095 * self.fuel.sigma) ** (-1) * np.exp((0.792 + 3.7597 * self.fuel.sigma ** 0.5) * (b + 0.1)) # (42)CORRIGIDO
C = 7.47 * np.exp(-0.8711 * (self.fuel.sigma ** 0.55)) # CORRIGIDO
B = 0.15988 * (self.fuel.sigma ** 0.54) # CORRIGIDO
E = 0.715 * np.exp(-0.01094 * self.fuel.sigma) # CORRIGIDO
Pi_w = C * (3.281 * self.U ** B) * ((b / bop) ** (-E)) # 47 CORRIGIDO
Pi_s = 5.275 * (b ** (-0.3)) * (np.tan(self.pi) ** (2)) # 51 CORRIGIDO
EE = np.exp(-4.528 / self.fuel.sigma) # 14 CORRIGIDO
Qig = 581 + 2594 * self.Mf
R = Ir * Xi * (1 + Pi_w + Pi_s) / (Rho_b * EE * Qig)
R = round(R, 6)
return R
def Harmonic_ROS(self, neighbor_cell):
return (2 * self.ROS * neighbor_cell.ROS) / (self.ROS + neighbor_cell.ROS)
def Harmonic_ROS_time_adjacent(self, neighbor_cell):
if self.Harmonic_ROS(neighbor_cell) == 0:
return sys.maxsize #make this better
return round((LENGTH_EDGE / self.Harmonic_ROS(neighbor_cell)) / 60, 6) # min
def Harmonic_ROS_time_sub_adjacent(self, neighbor_cell):
if self.Harmonic_ROS(neighbor_cell) == 0:
return sys.maxsize #make this better
return round(((math.sqrt(2) * LENGTH_EDGE) / self.Harmonic_ROS(neighbor_cell)) / 60, 6) # min
GRID_HEIGHT = 8
GRID_WIDTH = 8
TOTAL_TIME = 5
LENGTH_EDGE = 25 # m I should do differently. cells' edge should equal to this number
adjacent_neighbours = ((-1, 0), (0, -1), (0, 1), (1, 0))
sub_adjacent_neighbours = ((-1, -1), (-1, 1), (1, -1), (1, 1))
all_neighbours = ((-1, 0), (0, -1), (0, 1), (1, 0), (-1, -1), (-1, 1), (1, -1), (1, 1))
grid = np.empty((GRID_HEIGHT, GRID_WIDTH), dtype=Cell)
state_grid = np.empty((GRID_HEIGHT, GRID_WIDTH), dtype=int)
fuel_grid = np.empty((GRID_HEIGHT, GRID_WIDTH), dtype=str)
results = []
for y in range(0, GRID_HEIGHT):
for x in range(0, GRID_WIDTH):
new_random_number = random.randint(0, 100)
if new_random_number <= 1:
# water cell
new_cell = Cell(0, 0, 0)
fuel_grid[y][x] = 'EMPTY'
elif new_random_number > 1 and new_random_number <= 15:
# grass cell
new_cell = Cell(12.8, 2.82, 0, grass)
fuel_grid[y][x] = 'GRASS'
elif new_random_number > 15 and new_random_number <= 65:
# shrub cell
new_cell = Cell(12.8, 2.82, 0, shrub)
fuel_grid[y][x] = 'SHRUB'
else:
# conifer cell
new_cell = Cell(12.8, 2.82, 0, conifer)
fuel_grid[y][x] = 'CONIFER'
grid[y][x] = new_cell
state_grid[y][x] = new_cell.state
grid[2][3].state = BURNING
grid[4][3].state = BURNING
print('STATES at beginning ')
print(get_states(grid))
print()
for time in range(0, TOTAL_TIME):
time_to_burn_compare = []
for y in range(1, GRID_HEIGHT-1):
for x in range(1, GRID_WIDTH-1):
# if grid[y][x].time_to_burn <= 0:
#
# grid[y][x].time_to_burn = sys.maxsize
if grid[y][x].state == BURNED:
grid[y][x].time_to_burn = sys.maxsize
grid[y][x].has_burning_neighbour = False # maybe we should not put here?
if grid[y][x].state == FLAMMABLE: #checking flammable cells' neighbors
current_shortest_time_pos = (None, None)
Harmonic_ROS_times_to_compare = []
for (dx, dy) in all_neighbours:
# if grid[y+dy][x+dx].state == BURNED or grid[y+dy][x+dx].state == EMPTY or grid[y+dy][x+dx].state == FLAMMABLE:
# #maybe I need to say like if == 0,1,3 don't calculate H.ROS but how?
# continue
if grid[y+dy][x+dx].state == BURNING:
grid[y][x].has_burning_neighbour = True
if (dx, dy) in adjacent_neighbours:
Harmonic_ROS_time = grid[y][x].Harmonic_ROS_time_adjacent(grid[y+dy][x+dx])
if (dx, dy) in sub_adjacent_neighbours:
Harmonic_ROS_time = grid[y][x].Harmonic_ROS_time_sub_adjacent(grid[y+dy][x+dx])
Harmonic_ROS_times_to_compare.append(Harmonic_ROS_time)
if grid[y][x].has_burning_neighbour:
Harmonic_ROS_times_to_compare.append(grid[y][x].time_to_burn)
print('Harmonic_ROS_times_to_compare:')
print(Harmonic_ROS_times_to_compare)
grid[y][x].time_to_burn = min(Harmonic_ROS_times_to_compare)
time_to_burn_compare.append(grid[y][x].time_to_burn)
else:
grid[y][x].time_to_burn = sys.maxsize
print('cell time_to_burn(min):')
print(grid[y][x].time_to_burn)
print()
print('time_to_burn_compare: ')
print(time_to_burn_compare)
if min(time_to_burn_compare) == sys.maxsize:
current_shortest_time = 0
print('current shortest time is 0, there is a problem?')
else:
current_shortest_time = min(time_to_burn_compare)
print()
print('current_shortest_time(min):')
print(current_shortest_time)
for y in range(1, GRID_HEIGHT - 1):
for x in range(1, GRID_WIDTH - 1):
if grid[y][x].time_to_burn != None:
#if grid[y][x].has_burning_neighbour:
grid[y][x].time_to_burn = grid[y][x].time_to_burn - current_shortest_time
if grid[y][x].state == BURNING:
grid[y][x].state = BURNED
grid[y][x].time_to_burn = sys.maxsize #i should make this better
elif grid[y][x].state == FLAMMABLE and grid[y][x].time_to_burn <= 0 and grid[y][x].time_to_burn != None:
grid[y][x].state = BURNING
#grid[y][x].time_to_burn = sys.maxsize #i should make this better because i can not track which cell has 0 timetoburn
print()
print('time_to_burn_after_reducing_CurrentShortestTime(min):')
print(grid[y][x].time_to_burn)
results.append(get_states(grid))
print()
print('STATES at end of TİME: ' + str(time))
print(get_states(grid))
# print()
# print('TTB at end of TİME: ' + str(time))
# print(get_time_to_burn(grid))
# print()
print(fuel_grid)
fig, ax = plt.subplots(figsize=(5, 8))
cmap = colors.ListedColormap(['mediumblue','limegreen', 'red', 'darkred'])
def update(i):
im_normed = results[i]
ax.imshow(im_normed, cmap=cmap)
ax.set_title("Time step: {}".format(i), fontsize=20)
ax.set_axis_off()
anim = animation.FuncAnimation(fig, update, frames=np.arange(0, 5), interval=10)
anim.save('colour_rotation.gif', dpi=80, writer='imagemagick')
plt.close()
print("finished")
Editor is loading...