Untitled

mail@pastecode.io avatar
unknown
python
2 years ago
10 kB
2
Indexable
Never
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")