diff --git a/example/example.py b/example/example.py index 08d09a7..fe23b69 100644 --- a/example/example.py +++ b/example/example.py @@ -1,11 +1,14 @@ +from datetime import timedelta from time import time +import matplotlib.pyplot as plt import numpy as np from propagator.core import ( # type: ignore FUEL_SYSTEM_LEGACY, BoundaryConditions, Propagator, + PropagatorOutOfBoundsError, ) veg = np.full((2000, 2000), 5, dtype=np.int32) @@ -26,42 +29,59 @@ ignition_array[center_x, center_y] = 1 -boundary_conditions_list: list[BoundaryConditions] = [ - BoundaryConditions( - time=0, - ignition_mask=ignition_array, # type: ignore - wind_speed=np.ones(dem.shape) * 40, # km/h - wind_dir=np.ones(dem.shape) * 90, # degrees from north - moisture=np.ones(dem.shape) * 0, # percentage - ), -] -for boundary_condition in boundary_conditions_list: - simulator.set_boundary_conditions(boundary_condition) +boundary_condition = BoundaryConditions( + time=0, + ignition_mask=ignition_array, # type: ignore + wind_speed=np.ones(dem.shape) * 40, # km/h + wind_dir=np.ones(dem.shape) * 90, # degrees from north + moisture=np.ones(dem.shape) * 0, # percentage +) +simulator.set_boundary_conditions(boundary_condition) start_time = time() -while simulator.time < 3600 * 60: +step_time_init = time() +while simulator.time < 3600 * 24: next_time = simulator.next_time() if next_time is None: break + try: + simulator.step() + except PropagatorOutOfBoundsError: + print("Fire reached out of bounds area, stopping simulation.") + break + finally: + step_time_end = time() + if simulator.time % 3600 == 0: + print( + f"Time: {timedelta(seconds=int(simulator.time))} | elapsed: {step_time_end - step_time_init} seconds" + ) + + # create a plot of the fire probability + output = simulator.get_output() + fire_prob = output.fire_probability + ros_mean = output.ros_mean - step_time_init = time() - simulator.step() - step_time_end = time() - if simulator.time % 3600 == 0: - print( - f"Time: {simulator.time} | elapsed: {step_time_end - step_time_init} seconds" - ) + plt.figure(figsize=(8, 6)) + plt.imshow(fire_prob, cmap="hot", vmin=0, vmax=1) + plt.colorbar(label="Fire Probability") + plt.title( + f"Fire Probability at time {timedelta(seconds=int(simulator.time))}" + ) + plt.savefig( + f"example/output/fire_probability_{simulator.time}.png" + ) + plt.close() - # create a plot of the fire probability - fire_prob = simulator.compute_fire_probability() - import matplotlib.pyplot as plt + plt.figure(figsize=(8, 6)) + plt.imshow(ros_mean, cmap="hot") + plt.colorbar(label="Rate of Spread (mean)") + plt.title( + f"Rate of Spread (mean) at time {timedelta(seconds=int(simulator.time))}" + ) + plt.savefig(f"example/output/rate_of_spread_{simulator.time}.png") + plt.close() - plt.figure(figsize=(8, 6)) - plt.imshow(fire_prob, cmap="hot", vmin=0, vmax=1) - plt.colorbar(label="Fire Probability") - plt.title(f"Fire Probability at time {simulator.time} seconds") - plt.savefig(f"example/output/fire_probability_{simulator.time}.png") - plt.close() + step_time_init = time() end_time = time() print(f"Simulation completed in {end_time - start_time} seconds.") diff --git a/example/example_stochastic_ignitions.py b/example/example_stochastic_ignitions.py new file mode 100644 index 0000000..fb63877 --- /dev/null +++ b/example/example_stochastic_ignitions.py @@ -0,0 +1,94 @@ +from datetime import timedelta +from random import random +from time import time + +import matplotlib.pyplot as plt +import numpy as np + +from propagator.core import ( # type: ignore + FUEL_SYSTEM_LEGACY, + BoundaryConditions, + Propagator, + PropagatorOutOfBoundsError, +) + +veg = np.full((2000, 2000), 5, dtype=np.int32) +dem = np.zeros((2000, 2000), dtype=np.float32) +N_REALIZATIONS = 30 + +simulator = Propagator( + dem=dem, + veg=veg, + realizations=N_REALIZATIONS, + fuels=FUEL_SYSTEM_LEGACY, + do_spotting=False, + out_of_bounds_mode="ignore", +) + +ignition_array = np.zeros(dem.shape + (N_REALIZATIONS,), dtype=np.uint8) +center_x, center_y = dem.shape[0] // 2, dem.shape[1] // 2 +for r in range(N_REALIZATIONS): + # set central pixel as ignition point + x, y = ( + center_x + int(random() * 400) - 200, + center_y + int(random() * 400) - 200, + ) + ignition_array[x, y, r] = 1 + + +boundary_condition = BoundaryConditions( + time=0, + ignition_mask=ignition_array, # type: ignore + wind_speed=np.ones(dem.shape) * 40, # km/h + wind_dir=np.ones(dem.shape) * 90, # degrees from north + moisture=np.ones(dem.shape) * 0, # percentage +) +simulator.set_boundary_conditions(boundary_condition) + +start_time = time() +step_time_init = time() +while simulator.time < 3600 * 24: + next_time = simulator.next_time() + if next_time is None: + break + try: + simulator.step() + except PropagatorOutOfBoundsError: + print("Fire reached out of bounds area, stopping simulation.") + break + finally: + step_time_end = time() + if simulator.time % 3600 == 0: + print( + f"Time: {timedelta(seconds=int(simulator.time))} | elapsed: {step_time_end - step_time_init} seconds" + ) + + # create a plot of the fire probability + output = simulator.get_output() + fire_prob = output.fire_probability + ros_mean = output.ros_mean + + plt.figure(figsize=(8, 6)) + plt.imshow(fire_prob, cmap="hot", vmin=0, vmax=1) + plt.colorbar(label="Fire Probability") + plt.title( + f"Fire Probability at time {timedelta(seconds=int(simulator.time))}" + ) + plt.savefig( + f"example/output/fire_probability_{simulator.time}.png" + ) + plt.close() + + plt.figure(figsize=(8, 6)) + plt.imshow(ros_mean, cmap="hot") + plt.colorbar(label="Rate of Spread (mean)") + plt.title( + f"Rate of Spread (mean) at time {timedelta(seconds=int(simulator.time))}" + ) + plt.savefig(f"example/output/rate_of_spread_{simulator.time}.png") + plt.close() + + step_time_init = time() + +end_time = time() +print(f"Simulation completed in {end_time - start_time} seconds.") diff --git a/src/propagator/core/models.py b/src/propagator/core/models.py index 9912536..d12a58c 100644 --- a/src/propagator/core/models.py +++ b/src/propagator/core/models.py @@ -173,7 +173,7 @@ class BoundaryConditions: wind_speed : Optional[npt.NDArray[np.floating]] Wind speed map (km/h). ignition_mask : Optional[npt.NDArray[np.bool_]] - Boolean mask of new ignition points. + 2D or 3D Boolean mask of new ignition points. additional_moisture : Optional[npt.NDArray[np.floating]] Extra moisture to add to fuel (%), can be sparse. vegetation_changes : Optional[npt.NDArray[np.floating]] diff --git a/src/propagator/core/propagator.py b/src/propagator/core/propagator.py index cabec0b..fa520e2 100644 --- a/src/propagator/core/propagator.py +++ b/src/propagator/core/propagator.py @@ -296,10 +296,19 @@ def set_boundary_conditions( if boundary_condition.ignition_mask is not None: ign_arr = boundary_condition.ignition_mask - points = np.argwhere(ign_arr > 0) - points_repeated = np.repeat(points, self.realizations, axis=0) - realizations = np.tile(np.arange(self.realizations), len(points)) + + points = np.argwhere(ign_arr) + + # check ignition_mask shape beforehand + if len(boundary_condition.ignition_mask.shape) == 2: + points_repeated = np.repeat(points, self.realizations, axis=0) + realizations = np.tile( + np.arange(self.realizations), len(points) + ) + else: + points_repeated = points + realizations = points[:, 2] fireline_intensity = np.zeros_like( points_repeated[:, 0], dtype=np.float32