[ARC] Improve SWASH output handling and scatter matrix visualization

bluemath_tk/core/plotting/scatter.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import pandas as pd
+import matplotlib.pyplot as plt
 from matplotlib.axes import Axes
 from matplotlib.figure import Figure
 from scipy.stats import gaussian_kde, probplot
@@ -124,6 +125,8 @@ def validation_scatter(
 def plot_scatters_in_triangle(
     dataframes: List[pd.DataFrame],
     data_colors: Optional[List[str]] = None,
+    variables: Optional[List[str]] = None,
+    color: Optional[str] = None,
     **kwargs,
 ) -> Tuple[Figure, np.ndarray]:
     """
@@ -132,15 +135,28 @@ def plot_scatters_in_triangle(
     Parameters
     ----------
     dataframes : List[pd.DataFrame]
-        List of dataframes to plot. Each dataframe should contain the same columns.
+        List of dataframes to plot. Each dataframe should contain the variables
+        specified in `variables`.
     data_colors : Optional[List[str]], optional
         List of colors for the dataframes. If None, uses default_colors.
+        Ignored when `color` is provided.
+    variables : Optional[List[str]], optional
+        List of dataframe columns to use as scatter plot axes. If None, all
+        columns from the first dataframe are used.
+    color : Optional[str], optional
+        Name of the dataframe column used to color the scatter points.
+        If provided, a shared colorbar is added to the figure.
+        If None, each dataframe is plotted using its corresponding color
+        from `data_colors`.
     **kwargs : dict, optional
         Additional keyword arguments for the scatter plot. These will be passed to
         matplotlib.pyplot.scatter. Common parameters include:
         - s : float, marker size
         - alpha : float, transparency
         - marker : str, marker style
+        - cmap : str, colormap used when `color` is provided
+        - vmin : float, minimum value for colormap normalization
+        - vmax : float, maximum value for colormap normalization
 
     Returns
     -------
@@ -152,14 +168,18 @@ def plot_scatters_in_triangle(
     Raises
     ------
     ValueError
-        If the variables in the first dataframe are not present in all other dataframes.
+        If the specified variables or color variable are not present in all
+        dataframes.
     """
 
     if data_colors is None:
         data_colors = default_colors
 
-    # Get the number and names of variables from the first dataframe
-    variables_names = list(dataframes[0].columns)
+    # Get the number and names of variables
+    if variables is None:
+        variables = list(dataframes[0].columns)
+
+    variables_names = variables
     num_variables = len(variables_names)
 
     # Check variables names are in all dataframes
@@ -168,6 +188,18 @@ def plot_scatters_in_triangle(
             raise ValueError(
                 f"Variables {variables_names} are not in dataframe {df.columns}."
             )
+        if color is not None and color not in df.columns:
+            raise ValueError(
+                f"Color variable '{color}' is not in dataframe {df.columns}."
+            )
+
+    # Set common color limits when coloring by a dataframe variable
+    if color is not None:
+        color_values = np.concatenate(
+            [df[color].dropna().values for df in dataframes]
+        )
+        kwargs.setdefault("vmin", np.nanmin(color_values))
+        kwargs.setdefault("vmax", np.nanmax(color_values))
 
     # Create figure and axes
     default_static_plot = DefaultStaticPlotting()
@@ -177,6 +209,7 @@ def plot_scatters_in_triangle(
         sharex=False,
         sharey=False,
     )
+
     if isinstance(axes, Axes):
         axes = np.array([[axes]])
 
@@ -187,10 +220,11 @@ def plot_scatters_in_triangle(
                     ax=axes[c2, c1],
                     x=df[v1],
                     y=df[v2],
-                    c=data_colors[idf],
+                    c=df[color] if color is not None else data_colors[idf],
                     alpha=0.6,
                     **kwargs,
                 )
+
             if c1 == c2:
                 axes[c2, c1].set_xlabel(variables_names[c1 + 1])
                 axes[c2, c1].set_ylabel(variables_names[c2])
@@ -200,4 +234,22 @@ def plot_scatters_in_triangle(
             else:
                 fig.delaxes(axes[c2, c1])
 
+    # Add shared colorbar when coloring by a dataframe variable
+    if color is not None:
+        cmap = kwargs.get("cmap", "viridis")
+        norm = plt.Normalize(
+            vmin=kwargs["vmin"],
+            vmax=kwargs["vmax"],
+        )
+        sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
+        sm.set_array([])
+
+        visible_axes = [ax for ax in axes.flat if ax in fig.axes]
+
+        fig.colorbar(
+            sm,
+            ax=visible_axes,
+            label=color,
+        )
+
     return fig, axes

bluemath_tk/wrappers/_base_wrappers.py

@@ -443,18 +443,20 @@ def build_case_and_render_files(
             case_context=case_context,
             case_dir=case_dir,
         )
-        for template_name in self.templates_name:
-            try:
-                self.render_file_from_template(
-                    template_name=template_name,
-                    context=case_context,
-                    output_filename=op.join(case_dir, template_name),
-                )
-            except UnicodeDecodeError as _ude:
-                self.copy_files(
-                    src=op.join(self.templates_dir, template_name),
-                    dst=op.join(case_dir, template_name),
-                )
+
+        # UNCOMMENT!
+        #for template_name in self.templates_name:
+        #    try:
+        #        self.render_file_from_template(
+        #            template_name=template_name,
+        #            context=case_context,
+        #            output_filename=op.join(case_dir, template_name),
+        #        )
+        #    except UnicodeDecodeError as _ude:
+        #        self.copy_files(
+        #            src=op.join(self.templates_dir, template_name),
+        #            dst=op.join(case_dir, template_name),
+        #        )
 
     def build_cases(
         self,

bluemath_tk/wrappers/sfincs/sfincs_wrapper.py

@@ -8,6 +8,8 @@
 
 from .._base_wrappers import BaseModelWrapper
 
+from pyproj import CRS
+import rioxarray
 
 class SfincsModelWrapper(BaseModelWrapper):
     """
@@ -94,15 +96,20 @@ def setup_infiltration(self, sf: SfincsModel, case_context: dict) -> List[dict]:
         Setup the infiltration for the SFINCS model.
         """
 
-        p_infiltration = case_context.get("path_to_inf_tif")
-        ant_moisture = "avg"
+        #ant_moisture = case_context.get("antecedent_moisture_conditions")
+
+        #p_infiltration = case_context.get("path_to_inf_tif") / f"highres_infiltration_oahu_{ant_moisture}.tif"
 
-        dataset_inf = sf.data_catalog.get_rasterdataset(p_infiltration)
+        dataset_inf = sf.data_catalog.get_rasterdataset(case_context.get("path_to_inf_tif"))
+        #dataset_inf = rioxarray.open_rasterio(p_infiltration, masked=True).squeeze(drop=True)
 
-        dataset_inf.name = "cn_{0}".format(ant_moisture)
+        #dataset_inf.raster.set_spatial_dims(x_dim="x", y_dim="y")
+        #dataset_inf.raster.set_crs(CRS.from_epsg(26904)) 
+
+        dataset_inf.name = "cn"#.format(ant_moisture)
 
         sf.setup_cn_infiltration(
-            dataset_inf.compute(), antecedent_moisture="{0}".format(ant_moisture)
+            dataset_inf.compute(), antecedent_moisture="cn"
         )
 
         return dataset_inf
@@ -215,36 +222,58 @@ def build_template_case(self) -> None:
         self.setup_outflow(sf=sf, case_context=self.fixed_parameters)
 
         self.setup_waterlevel_mask(sf=sf, case_context=self.fixed_parameters)
-
+        
         _ = sf.plot_basemap(bmap="sat", zoomlevel=12)
 
+        sf.config["dtout"] = 900
+        sf.config["storemeteo"] = 1
+        sf.config["bzifile"] = 'sfincs.bzi'
+        sf.config["bzsfile"] = 'sfincs.bzs'
+        sf.config["tstop"] = '20100201 000000'
+        sf.config["tstart"] = '20100131 000000'
+        sf.config["tref"] = '20100131 000000'
+        sf.config["stop"] = '20100201 000000'
+        sf.config["netamprfile"] = "precip_2d.nc"
+        sf.config["dtmaxout"] = None
+        sf.config["bndfile"] = "sfincs.bnd"
+
         sf.write()
 
+        self.sf = sf
+
+
     def build_case(self, case_context: dict, case_dir: str) -> None:
         """
         Build the base SFINCS model. This includes setting up the grid,
         depth, friction, mask, outflow and waterlevel mask. It also
         applies the precipitation and waterlevel forcing if specified.
         """
 
-        sf = SfincsModel(root=case_dir, mode="w+")
+        #sf = SfincsModel(root=case_dir, mode="w+")
 
-        sf.setup_grid(
-            x0=case_context["x0"],
-            y0=case_context["y0"],
-            dx=case_context["dx"],
-            dy=case_context["dy"],
-            nmax=case_context["nmax"],
-            mmax=case_context["mmax"],
-            rotation=case_context["rotation"],
-            epsg=case_context["epsg"],
-        )
-        tstart, tstop = self.set_ctimes(case_context=case_context)
+        #sf.setup_grid(
+        #    x0=case_context["x0"],
+        #    y0=case_context["y0"],
+        #    dx=case_context["dx"],
+        #    dy=case_context["dy"],
+        #    nmax=case_context["nmax"],
+        #    mmax=case_context["mmax"],
+        #    rotation=case_context["rotation"],
+        #    epsg=case_context["epsg"],
+        #)
+
+        sf = self.sf
 
-        sf.config["tstop"] = tstop
-        sf.config["tstart"] = tstart
-        sf.config["dtout"] = 60
-        sf.config["storemeteo"] = 1
+        sf.set_root(root=case_dir, mode="r+")
+
+        #datasets_dep = self.setup_dem(sf=sf, case_context=self.fixed_parameters)
+
+        #tstart, tstop = self.set_ctimes(case_context=case_context)
+
+        #sf.config["tstop"] = tstop
+        #sf.config["tstart"] = tstart
+        #sf.config["dtout"] = 60
+        #sf.config["storemeteo"] = 1
 
         if case_context.get("quickly_waterlevel_forcing") is not None:
             """
@@ -259,7 +288,7 @@ def build_case(self, case_context: dict, case_dir: str) -> None:
             )
             sf.write_forcing()
             os.rename(op.join(case_dir, "sfincs.bzs"), op.join(case_dir, "sfincs.bzi"))
-
+        
         if case_context.get("slowly_waterlevel_forcing") is not None:
             sf.setup_waterlevel_forcing(
                 timeseries=case_context.get("slowly_waterlevel_forcing"),
@@ -271,14 +300,22 @@ def build_case(self, case_context: dict, case_dir: str) -> None:
             sf.setup_precip_forcing_from_grid(
                 precip=case_context.get("precipitation_forcing"), aggregate=False
             )
+
+            sf.write_forcing()
 
-        if case_context.get("gdf_crs") is not None:
-            sf.setup_observation_lines(
-                locations=case_context.get("gdf_crs"), merge=False
-            )
+        #if case_context.get("gdf_crs") is not None:
+        #    sf.setup_observation_lines(
+        #        locations=case_context.get("gdf_crs"), merge=False
+        #    )
+
+        #if case_context.get("gdf_obs") is not None:
+        #    sf.setup_observation_points(locations=case_context.get("gdf_obs"))
 
-        if case_context.get("gdf_obs") is not None:
-            sf.setup_observation_points(locations=case_context.get("gdf_obs"))
+        #if case_context.get("antecedent_moisture_conditions") is not None:
+        #    self.setup_infiltration(
+        #        sf=sf,
+        #        case_context=case_context
+        #    )
 
         # if case_context.get("precipitation_forcing") is not None and case_context.get("waterlevel_forcing") is not None:
         #    self.setup_rivers(sf)
@@ -288,4 +325,4 @@ def build_case(self, case_context: dict, case_dir: str) -> None:
 
         # sf.write_forcing()
 
-        sf.write()
+        #sf.write()

bluemath_tk/wrappers/swash/swash_wrapper.py

@@ -254,25 +254,30 @@ def _convert_case_output_files_to_nc(
         """
 
         df_output = self._read_tabfile(file_path=output_path)
-        df_output[["Xp", "Yp", "Tsec"]] = df_output[["Xp", "Yp", "Tsec"]].astype(
-            int
-        )  # TODO: check if this is correct
-        df_output.set_index(
-            ["Xp", "Yp", "Tsec"], inplace=True
-        )  # set index to Xp, Yp and Tsec
+
+        # Set water level to NaN at dry points
+        dry_threshold = 0.05
+        dry_points = df_output["Depth"].round(12) <= dry_threshold
+        df_output.loc[dry_points, "Watlev"] = np.nan
+
+        df_output[["Xp"]] = df_output[["Xp"]].astype(int)
+        df_output["Tsec"] = df_output["Tsec"].round().astype(int)
+
+        df_output.set_index(["Xp", "Tsec"], inplace=True)
         ds_output = df_output.to_xarray()
 
+        # Read runup output
         df_run = self._read_tabfile(file_path=run_path)
-        df_run[["Tsec"]] = df_run[["Tsec"]].astype(
-            int
-        )  # TODO: check if this is correct
+
+        df_run["Tsec"] = df_run["Tsec"].round().astype(int)
+
         df_run.set_index(["Tsec"], inplace=True)
         ds_run = df_run.to_xarray()
 
-        # merge output files to one xarray.Dataset
+        # Merge output files into one xarray.Dataset
         ds = xr.merge([ds_output, ds_run], compat="no_conflicts")
 
-        # assign correct coordinate case_num
+        # Assign case number
         ds.coords["case_num"] = case_num
 
         return ds

bluemath_tk/wrappers/swash/templates/INPUT

@@ -48,7 +48,7 @@ QUANTITY XP hexp=10
 QUANT RUNUP delrp 0.01
 $
 CURVE 'line' 0 0 {{ xlenc }} {{ xlenc }} 0
-TABLE 'line' HEAD 'output.tab' TSEC XP YP WATL OUTPUT 0 1 SEC
+TABLE 'line' HEAD 'output.tab' TSEC XP DEPTH WATL OUTPUT 0 1 SEC
 $
 TABLE 'NOGRID' HEAD 'run.tab' TSEC RUNUP OUTPUT 0 1 SEC
 $