Merge pull request #90 from Deltares/feature/raster_risk

GridModel: Risk

.testdata/create_test_data.py

@@ -385,7 +385,7 @@ def create_settings_grid():
             "keep_temp_files": True,
         },
         "output": {
-            "path": "output/event",
+            "path": "output/event_grid",
             "grid": {"name": "output.nc"},
         },
         "hazard": {
@@ -415,6 +415,16 @@ def create_settings_grid():
     with open(Path(p, "settings_grid.toml"), "wb") as f:
         tomli_w.dump(doc, f)
 
+    doc_r = copy.deepcopy(doc)
+    doc_r["output"]["path"] = "output/risk_grid"
+    doc_r["hazard"]["file"] = "hazard/risk_map.nc"
+    doc_r["hazard"]["return_periods"] = [2, 5, 10, 25]
+    doc_r["hazard"]["settings"] = {"var_as_band": True}
+    doc_r["hazard"]["risk"] = True
+
+    with open(Path(p, "settings_grid_risk.toml"), "wb") as f:
+        tomli_w.dump(doc_r, f)
+
 
 def create_settings_risk():
     """_summary_."""

src/fiat/cfg.py

@@ -57,6 +57,10 @@ def __init__(
         # Create the hidden temporary folder
         self._create_temp_dir()
 
+        # Create risk directory if needed
+        if self.get("hazard.risk"):
+            self._create_risk_dir()
+
         # Set the cache size per GDAL object
         _cache_size = self.get("global.gdal_cache")
         if _cache_size is not None:
@@ -109,6 +113,14 @@ def _create_output_dir(
         generic_folder_check(_p)
         self["output.path"] = _p
 
+    def _create_risk_dir(
+        self,
+    ):
+        """_summary_."""
+        _ph = Path(self["output.path"], "rp_damages")
+        generic_folder_check(_ph)
+        self["output.path.risk"] = _ph
+
     def _create_temp_dir(
         self,
     ):

src/fiat/check.py

@@ -210,10 +210,11 @@ def check_hazard_rp_iden(
     if deter_type(bn_str, l - 1) != 3:
         return [float(n) for n in bnames]
 
-    if len(rp_cfg) == len(bnames):
-        rp_str = "\n".join([str(n) for n in rp_cfg]).encode()
-        if deter_type(rp_str, l - 1) != 3:
-            return rp_cfg
+    if rp_cfg is not None:
+        if len(rp_cfg) == len(bnames):
+            rp_str = "\n".join([str(n) for n in rp_cfg]).encode()
+            if deter_type(rp_str, l - 1) != 3:
+                return rp_cfg
 
     logger.error(
         f"'{path.name}': cannot determine the return periods for the risk calculation"

src/fiat/io.py

@@ -6,7 +6,6 @@
 import weakref
 from abc import ABCMeta, abstractmethod
 from io import BufferedReader, BufferedWriter, FileIO, TextIOWrapper
-from itertools import product
 from math import floor, log10
 from typing import Any
 
@@ -558,8 +557,6 @@ def __init__(
         else:
             ValueError("")
 
-        self.create_windows()
-
     def __iter__(self):
         self.flush()
         self._reset_chunking()
@@ -637,24 +634,6 @@ def shape(self):
         """_summary_."""
         return self._x, self._y
 
-    def create_windows(self):
-        """_summary_."""
-        _lu = tuple(
-            product(
-                range(0, self._x, self._chunk[0]),
-                range(0, self._y, self._chunk[1]),
-            ),
-        )
-        for _l, _u in _lu:
-            w = min(self._chunk[0], self._x - _l)
-            h = min(self._chunk[1], self._y - _u)
-            yield (
-                _l,
-                _u,
-                w,
-                h,
-            )
-
     def get_metadata_item(
         self,
         entry: str,
@@ -1269,6 +1248,13 @@ def get_srs(self):
         """_summary_."""
         return self.src.GetSpatialRef()
 
+    def set_chunk_size(
+        self,
+        chunk: tuple,
+    ):
+        """_summary_."""
+        self._chunk = chunk
+
     @_BaseIO._check_mode
     def set_geotransform(self, affine: tuple):
         """_summary_."""

src/fiat/models/base.py

@@ -3,7 +3,6 @@
 from abc import ABCMeta, abstractmethod
 from multiprocessing import Manager
 from os import cpu_count
-from pathlib import Path
 
 from osgeo import osr
 
@@ -129,11 +128,6 @@ def _read_hazard_grid(self):
             # Directly calculate the coefficients
             rp_coef = calc_rp_coef(rp)
             self.cfg["hazard.rp_coefficients"] = rp_coef
-            # Set the risk folder for raster calculations
-            self.cfg["output.path.risk"] = Path(
-                self.cfg["output.path"],
-                "rp_damages",
-            )
 
         # Information for output
         ns = check_hazard_band_names(

src/fiat/models/calc.py

@@ -46,11 +46,11 @@ def calc_rp_coef(
     # Step 1: Compute frequencies associated with T-values.
     _rp = sorted(rp)
     idxs = [_rp.index(n) for n in rp]
-    rp.sort()
-    rp_l = len(rp)
+    rp_u = sorted(rp)
+    rp_l = len(rp_u)
 
-    f = [1 / n for n in rp]
-    lf = [math.log(1 / n) for n in rp]
+    f = [1 / n for n in rp_u]
+    lf = [math.log(1 / n) for n in rp_u]
 
     if rp_l == 1:
         return f
@@ -79,9 +79,9 @@ def calc_rp_coef(
         b[0]
         if idx == 0
         else f[idx] + a[idx - 1]
-        if idx == len(rp) - 1
+        if idx == rp_l - 1
         else a[idx - 1] + b[idx]
-        for idx in range(len(rp))
+        for idx in range(rp_l)
     ]
 
     return [alpha[idx] for idx in idxs]

src/fiat/models/grid.py

@@ -11,7 +11,7 @@
 from fiat.io import open_grid
 from fiat.log import spawn_logger
 from fiat.models.base import BaseModel
-from fiat.models.util import grid_worker_exact
+from fiat.models.util import grid_worker_exact, grid_worker_risk
 
 logger = spawn_logger("fiat.model.grid")
 
@@ -59,9 +59,19 @@ def _read_exposure_grid(self):
 
         self.exposure_grid = data
 
-    def resolve():
+    def resolve(self):
         """_summary_."""
-        pass
+        if self.cfg.get("hazard.risk"):
+            logger.info("Setting up risk calculations..")
+
+            # Time the function
+            _s = time.time()
+            grid_worker_risk(
+                self.cfg,
+                self.exposure_grid.chunk,
+            )
+            _e = time.time() - _s
+            logger.info(f"Risk calculation time: {round(_e, 2)} seconds")
 
     def run(self):
         """_summary_."""
@@ -87,6 +97,7 @@ def run(self):
                     )
                     futures.append(fs)
             logger.info("Busy...")
+            # Wait for the children to finish their calculations
             wait(futures)
 
         else:
@@ -106,7 +117,7 @@ def run(self):
             logger.info("Busy...")
             p.join()
         _e = time.time() - _s
-
         logger.info(f"Calculations time: {round(_e, 2)} seconds")
+        self.resolve()
         logger.info(f"Output generated in: '{self.cfg['output.path']}'")
         logger.info("Grid calculation are done!")