Feat/3d ignition mask

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.")

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.")

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]]

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