diff --git a/bluemath_tk/core/plotting/scatter.py b/bluemath_tk/core/plotting/scatter.py index a5c2438..d5ca57e 100644 --- a/bluemath_tk/core/plotting/scatter.py +++ b/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 diff --git a/bluemath_tk/wrappers/_base_wrappers.py b/bluemath_tk/wrappers/_base_wrappers.py index 8bd9f66..d33d073 100644 --- a/bluemath_tk/wrappers/_base_wrappers.py +++ b/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, diff --git a/bluemath_tk/wrappers/sfincs/sfincs_wrapper.py b/bluemath_tk/wrappers/sfincs/sfincs_wrapper.py index 79408c9..4e74a9a 100644 --- a/bluemath_tk/wrappers/sfincs/sfincs_wrapper.py +++ b/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,11 +222,26 @@ 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, @@ -227,24 +249,31 @@ def build_case(self, case_context: dict, case_dir: str) -> None: 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() diff --git a/bluemath_tk/wrappers/swash/swash_wrapper.py b/bluemath_tk/wrappers/swash/swash_wrapper.py index e2ae5cd..02d3240 100644 --- a/bluemath_tk/wrappers/swash/swash_wrapper.py +++ b/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 diff --git a/bluemath_tk/wrappers/swash/templates/INPUT b/bluemath_tk/wrappers/swash/templates/INPUT index a2569ae..6311b6c 100644 --- a/bluemath_tk/wrappers/swash/templates/INPUT +++ b/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 $