Merge pull request #81 from Deltares/hazard_sfincs

Hazard sfincs into fiat integrator
Deltares · Jun 1, 2023 · e45a36e · e45a36e
2 parents d1b6473 + 9b18b26
commit e45a36e
Show file tree

Hide file tree

Showing 13 changed files with 122 additions and 77 deletions.
diff --git a/hydromt_fiat/fiat.py b/hydromt_fiat/fiat.py
@@ -157,17 +157,17 @@ def setup_exposure_raster(self):
 
     def setup_hazard(
         self,
-        map_fn: str,
-        map_type: str,
-        rp,
-        crs,
-        nodata,
-        var,
-        chunks,
-        risk_output: bool = True,
-        hazard_type: str = "flooding",
+        map_fn:       Union[str, list[str]],
+        map_type:     Union[str, list[str]],
+        rp:           Union[int, list[int], None] = None ,
+        crs:          Union[str, list[str], None] = None,
+        nodata:       Union[int, list[int]] = None,
+        var:          Union[str, list[str], None] = None,
+        chunks:       Union[str, int, list[int]] = 'auto',
+        risk_output:  bool = True,
+        hazard_type:  str = "flooding",
         name_catalog: str = "flood_maps",
-        maps_id: str = "RP",
+        maps_id:      str = "RP",
 
     ):
         hazard = Hazard()
@@ -199,18 +199,19 @@ def setup_hazard(
         hazard_maps = []
         for hazard_map in self.maps.keys():
             hazard_maps.append(
-                str(Path("hazard") / (self.maps[hazard_map].name + ".nc"))
+                str(Path("hazard") / (hazard_map + ".nc"))
             )
 
         hazard_settings["grid_file"] = hazard_maps
-
+        hazard_settings["crs"] = hazard.crs
         if not isinstance(rp, list):
             rp = "Event"
         hazard_settings["return_period"] = rp
-
-        hazard_settings["crs"] = hazard.crs
-        hazard_settings["spatial_reference"] = map_type
-        self.config["hazard"] = hazard_settings
+        if map_type == "water_depth":
+            hazard_settings["spatial_reference"] = "DEM"
+        # else:
+        #     hazard_settings["spatial_reference"] = "DATUM"
+        self.config["hazard"]  = hazard_settings
 
 
 

diff --git a/hydromt_fiat/workflows/hazard.py b/hydromt_fiat/workflows/hazard.py
@@ -22,6 +22,8 @@ def checkInputs(
         nodata,
         var,
     ):
+
+
         # Checks the hazard input parameter types.
         # Checks map path list
         self.map_fn_lst = [map_fn] if isinstance(map_fn, (str, Path)) else map_fn
@@ -159,12 +161,23 @@ def readMaps(
                 )
 
             if os.path.exists(da_map_fn):
-                if not region.empty:
-                    da = model_fiat.data_catalog.get_rasterdataset(
-                        da_map_fn, geom=region, **kwargs
-                    )
+                if da_map_fn.stem == "sfincs_map":
+                    #da = model_fiat.data_catalog.get_rasterdataset(da_map_fn) The file cannot be directly read by HydroMT
+
+                    ds_map = xr.open_dataset(da_map_fn)
+                    da = ds_map[kwargs["variables"]].squeeze(dim="timemax").drop_vars("timemax")
+                    da.raster.set_crs(ds_map.crs.epsg_code)  
+                    da.raster.set_nodata(nodata=ds_map.encoding.get("_FillValue"))
+                    da.reset_coords(['spatial_ref'], drop=False)
+                    da.encoding["_FillValue"] = None
+                    #da_name = kwargs["variables"]
                 else:
-                    da = model_fiat.data_catalog.get_rasterdataset(da_map_fn, **kwargs)
+                    if not region.empty:
+                        da = model_fiat.data_catalog.get_rasterdataset(
+                            da_map_fn, geom=region, **kwargs
+                        )
+                    else:
+                        da = model_fiat.data_catalog.get_rasterdataset(da_map_fn, **kwargs)
             else:
                 if not region.empty:
                     da = model_fiat.data_catalog.get_rasterdataset(
@@ -206,9 +219,12 @@ def readMaps(
                 raise ValueError("The hazard map has no nodata value assigned.")
 
             # Correct (if necessary) the grid orientation from the lower to the upper left corner.
-            if da.raster.res[1] > 0:
-                da = da.reindex({da.raster.y_dim: list(reversed(da.raster.ycoords))})
-
+            if da_name != "sfincs_map":            #This check could not be implemented into the sfincs_map outputs
+                if da.raster.res[1] > 0:
+                    da = da.reindex(
+                        {da.raster.y_dim: list(reversed(da.raster.ycoords))}
+                        )
+
             # Check if the obtained hazard map is identical.
             if (
                 model_fiat.staticmaps
@@ -230,6 +246,9 @@ def readMaps(
                 else None
             )
 
+            if risk_output:
+                da = da.expand_dims({'rp': [da_rp]}, axis=0)
+
             if risk_output and da_rp is None:
                 # Get (if possible) the return period from dataset names if the input parameter is None.
                 if "rp" in da_name.lower():
@@ -240,7 +259,7 @@ def fstrip(x):
                     rp_str = "".join(
                         filter(fstrip, da_name.lower().split("rp")[-1])
                     ).lstrip("0")
-
+                    
                     try:
                         assert isinstance(
                             literal_eval(rp_str) if rp_str else None, (int, float)
@@ -308,71 +327,97 @@ def fstrip(x):
             #             "var": None if not hasattr(self, "var_lst") else self.var_lst[idx],
             #             "chunks": "auto" if chunks == "auto" else self.chunks_lst[idx],
             #            }
-
-            da = da.expand_dims({"rp": [da_rp]}, axis=0)
+            
+
 
             model_fiat.set_maps(da, da_name)
             # post = f"(rp {da_rp})" if rp is not None and risk_output else ""
             post = f"(rp {da_rp})" if risk_output else ""
             model_fiat.logger.info(f"Added {hazard_type} hazard map: {da_name} {post}")
 
-            # new_da = xr.concat([new_da, da], dim='rp')
+            # new_da = xr.concat([new_da, da], dim='rp')   
 
             # import numpy as np
             # new_da = xr.DataArray(np.zeros_like(da), dims=da.dims, coords=da.coords).rename('new_dataframe')
 
-        # # model_fiat.maps.to_netcdf("test_hazard_1/test.nc")
 
-        # #model_fiat.set_maps(catalog, "all")
-        maps = model_fiat.maps
-        list_keys = list(maps.keys())
-        maps_0 = maps[list_keys[0]].rename("risk")
-        list_keys.pop(0)
-
-        for idx, x in enumerate(list_keys):
-            key_name = list_keys[idx]
-            layer = maps[key_name]
-            maps_0 = xr.concat([maps_0, layer], dim="rp")
-
-        # new_da = maps_0.to_dataset(name='RISK')
-        # new_da.attrs = {  "returnperiod": list(list_rp),
-        #                   "type":self.map_type_lst,
-        #                   'name':list_names}
-
-        if not risk_output:
-            # new_da = maps_0.drop_dims('rp')
-            new_da = maps_0.to_dataset(name="EVENT")
-            new_da.attrs = {
-                "returnperiod": list(list_rp),
-                "type": self.map_type_lst,
-                "name": list_names,
-                "Analysis": "Event base",
-            }
-
-        else:
+        if risk_output:
+            maps = model_fiat.maps
+            list_keys = list(maps.keys())
+            maps_0 = maps[list_keys[0]].rename('risk')
+            list_keys.pop(0)
+
+            for idx, x in enumerate(list_keys):
+                key_name = list_keys[idx]
+                layer = maps[key_name]
+                maps_0 = xr.concat([maps_0, layer], dim='rp') 
+
             new_da = maps_0.to_dataset(name='RISK')
             new_da.attrs = {  "returnperiod": list(list_rp),
                             "type":self.map_type_lst,
                             'name':list_names,
                             "Analysis": "Risk"}  
 
-        model_fiat.hazard = new_da
-        model_fiat.set_maps(model_fiat.hazard, 'HydroMT_Fiat_hazard')
+            model_fiat.hazard = new_da
+            model_fiat.set_maps(model_fiat.hazard, 'HydroMT_Fiat_hazard')
 
-        list_maps = list(model_fiat.maps.keys())
+            list_maps = list(model_fiat.maps.keys())
 
-        if risk_output:
-            for item in list_maps[:-1]:
-                model_fiat.maps.pop(item)
+            if risk_output:
+                for item in list_maps[:-1]:
+                    model_fiat.maps.pop(item)
 
+        # else:
+        #     model_fiat.hazard = new_da
+
+
+        # # model_fiat.maps.to_netcdf("test_hazard_1/test.nc")
+
+        # #model_fiat.set_maps(catalog, "all")
+        # maps = model_fiat.maps
+        # list_keys = list(maps.keys())
+        # maps_0 = maps[list_keys[0]].rename('risk')
+        # list_keys.pop(0)
+
+        # for idx, x in enumerate(list_keys):
+        #     key_name = list_keys[idx]
+        #     layer = maps[key_name]
+        #     maps_0 = xr.concat([maps_0, layer], dim='rp') 
+
+        # if not risk_output:
+        #     # new_da = maps_0.drop_dims('rp')
+        #     new_da = maps_0.to_dataset(name='EVENT')
+        #     new_da.attrs = {  "returnperiod": list(list_rp),
+        #                     "type":self.map_type_lst,
+        #                     'name':list_names,
+        #                     "Analysis": "Event base"}  
+                                                            #TODO: A netcdf is not required when working with evetns       
+        # else:
+        #     new_da = maps_0.to_dataset(name='RISK')
+        #     new_da.attrs = {  "returnperiod": list(list_rp),
+        #                     "type":self.map_type_lst,
+        #                     'name':list_names,
+        #                     "Analysis": "Risk"}  
+
+        # if risk_output:
+        #     model_fiat.hazard = new_da
+        #     model_fiat.set_maps(model_fiat.hazard, 'HydroMT_Fiat_hazard')
+        # else:
+        #     model_fiat.hazard = new_da
+
+        # list_maps = list(model_fiat.maps.keys())
+
+        # if risk_output:
+        #     for item in list_maps[:-1]:
+        #         model_fiat.maps.pop(item)
+
 
         # ds = xr.Dataset(maps)
         # ds.raster.set_crs(da.raster.crs)
 
         # # new_da= ds.to_array(dim='rp')
 
         # # new_da = new_da.to_dataset(name='Full_stack')
-
         # # new_da.attrs = {  "returnperiod": list_rp,
         # #                   "type":self.map_type_lst,
         # #                   'name':list_names}
@@ -385,8 +430,7 @@ def fstrip(x):
         #     new_da = xr.concat([new_da, layer], dim='rp')
 
         # # new_da = new_da.to_dataset()
-
         # # config = model_fiat.config
         # # new_da.to_netcdf("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database/test_hazard_1/hazard/test_final_v2.nc")
         # #C:\Users\fuentesm\CISNE\HydroMT_sprint_sessions
-        return model_fiat.maps 
+        return model_fiat.maps
diff --git a/tests/test_SVI.py b/tests/test_SVI.py
@@ -4,7 +4,7 @@
 import pytest
 
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 _cases = {
     "fiat_objects": {
         "folder": "test_hazard",

diff --git a/tests/test_exposure.py b/tests/test_exposure.py
@@ -5,7 +5,7 @@
 import shutil
 
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "vulnerability_and_exposure": {

diff --git a/tests/test_hazard.py b/tests/test_hazard.py
@@ -17,7 +17,7 @@
     "fiat_objects": {
         "folder"   : "test_hazard_1",
         "ini"      : "test_hazard.ini",
-        "catalog"  : "fiat_catalog_hazard.yml",
+        "catalog"  : "fiat_catalog.yml",
 
     }
 }

diff --git a/tests/test_integration.py b/tests/test_integration.py
@@ -4,7 +4,7 @@
 import pytest
 import shutil
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "integration": {

diff --git a/tests/test_integrations_hazard.py b/tests/test_integrations_hazard.py
@@ -8,8 +8,8 @@
 
 _cases = {
     "integration": {
-        "data_catalogue": EXAMPLEDIR / "fiat_catalog_hazard.yml",
-        "dir": "test_hazard",
+        "data_catalogue": EXAMPLEDIR / "fiat_catalog.yml",
+        "dir": "test_hazard_integration",
         "ini": EXAMPLEDIR / "test_hazard.ini",
     },
 }

diff --git a/tests/test_raise_ground_floor_height.py b/tests/test_raise_ground_floor_height.py
@@ -3,7 +3,7 @@
 import pytest
 import shutil
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "raise_ground_floor_height": {

diff --git a/tests/test_reader.py b/tests/test_reader.py
@@ -3,7 +3,7 @@
 import pytest
 
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "read": {

diff --git a/tests/test_setup_new_composite_areas.py b/tests/test_setup_new_composite_areas.py
@@ -3,7 +3,7 @@
 import pytest
 import shutil
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "setup_new_composite_area_datum": {

diff --git a/tests/test_truncate_damage_function.py b/tests/test_truncate_damage_function.py
@@ -3,7 +3,7 @@
 import pytest
 import shutil
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "truncate_damage_function": {

diff --git a/tests/test_update_max_potential_damage.py b/tests/test_update_max_potential_damage.py
@@ -3,7 +3,7 @@
 import pytest
 import shutil
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "update_max_potential_damage": {

diff --git a/tests/test_vulnerability.py b/tests/test_vulnerability.py
@@ -4,7 +4,7 @@
 import pytest
 import shutil
 
-EXAMPLEDIR = Path().absolute() / "local_test_database"
+EXAMPLEDIR = Path("P:/11207949-dhs-phaseii-floodadapt/Model-builder/Delft-FIAT/local_test_database")
 
 _cases = {
     "vulnerability": {