From c07cfe55fb524a7ecdc301405664b4851461dc9e Mon Sep 17 00:00:00 2001 From: Maarten van Ormondt Date: Sun, 2 Mar 2025 09:21:02 +0100 Subject: [PATCH] lots of updates --- hydromt_sfincs/boundary_conditions.py | 695 ++++++++++++++++ hydromt_sfincs/config.py | 35 +- hydromt_sfincs/config_variables.py | 23 + hydromt_sfincs/discharge_points.py | 511 ++++++++++-- hydromt_sfincs/observation_points.py | 6 +- hydromt_sfincs/quadtree.py | 3 +- hydromt_sfincs/quadtree_mask.py | 2 +- hydromt_sfincs/sfincs.py | 48 +- .../snapwave_boundary_conditions.py | 769 ++++++++++++++++++ hydromt_sfincs/snapwave_conditions.py | 72 -- hydromt_sfincs/utils.py | 2 +- hydromt_sfincs/waterlevel_conditions.py | 70 -- hydromt_sfincs/wave_makers.py | 280 +++---- 13 files changed, 2140 insertions(+), 376 deletions(-) create mode 100644 hydromt_sfincs/boundary_conditions.py create mode 100644 hydromt_sfincs/snapwave_boundary_conditions.py delete mode 100644 hydromt_sfincs/snapwave_conditions.py delete mode 100644 hydromt_sfincs/waterlevel_conditions.py diff --git a/hydromt_sfincs/boundary_conditions.py b/hydromt_sfincs/boundary_conditions.py new file mode 100644 index 00000000..569da926 --- /dev/null +++ b/hydromt_sfincs/boundary_conditions.py @@ -0,0 +1,695 @@ +import geopandas as gpd +import numpy as np +import pandas as pd +import xarray as xr +from pathlib import Path +from typing import Union, List +import shapely +from pyproj import Transformer + +# from tabulate import tabulate + +from hydromt.model.components import ModelComponent +from hydromt.model import Model +from hydromt_sfincs import utils + + +class SfincsBoundaryConditions(ModelComponent): + def __init__( + self, + model: Model, + ): + # self._filename: str = "sfincs.bnd" # FIXME - List(str = "sfincs.bnd" and str = "sfincs.bzs" or str = "sfincs_netbndbzsbzi.nc") + self.data = gpd.GeoDataFrame() + super().__init__( + model=model, + ) + + # @property + # def data(self) -> gpd.GeoDataFrame: + # """Water level boundary conditions data. + + # Return pd.GeoDataFrame + # """ + # if self._data is None: + # self._initialize() + # return self._data + + # def _initialize(self) -> None: + # """Initialize boundary conditions data.""" + # if self._data is None: + # self._data = gpd.GeoDataFrame() + + def read(self, format: str = None): + """Read SFINCS boundary conditions (*.bnd, *.bzs, *.bca files) or netcdf file. + + The format of the boundary conditions files can be specified, + otherwise it is determined from the model configuration. + + Parameters + ---------- + format : str, optional + Format of the boundary conditions files, "asc" or "netcdf". + """ + + if format is None: + if self.model.config.get("netbndbzsbzifile"): + format = "netcdf" + else: + format = "asc" + + if format == "asc": + self.read_boundary_points() + # Check if there are any points + if not self.data.empty: + self.read_boundary_conditions_timeseries() + # Read astro if bcafile is defined + if self.model.config.get("bcafile"): + self.read_boundary_conditions_astro() + elif format == "netcdf": + # Read netcdf file + self.read_boundary_conditions_netcdf() + + def read_boundary_points(self, filename: str | Path = None): + """Read SFINCS boundary condition points (*.bnd) file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if crsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "bndfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if bnd file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition points file not found: {abs_file_path}" + ) + + # HydroMT does not have open_vector at the moment ... + # Read bnd file + # gdf = utils.read_xy(abs_file_path, crs=self.model.crs) + # # Add columns for timeseries and astro and add empty DataFrames + # gdf["timeseries"] = pd.DataFrame() + # gdf["astro"] = pd.DataFrame() + # # Add to self.data + # self.data = gdf + + # Read the bnd file + df = pd.read_csv( + abs_file_path, index_col=False, header=None, names=["x", "y"], sep="\s+" + ) + + gdf_list = [] + # Loop through points + for ind in range(len(df.x.values)): + name = str(ind + 1).zfill(4) + x = df.x.values[ind] + y = df.y.values[ind] + point = shapely.geometry.Point(x, y) + d = { + "name": name, + "timeseries": pd.DataFrame(), + "astro": pd.DataFrame(), + "geometry": point, + } + gdf_list.append(d) + self.data = gpd.GeoDataFrame(gdf_list, crs=self.model.crs) + + def read_boundary_conditions_timeseries(self, filename: str | Path = None): + """Read SFINCS boundary condition timeseries (*.bzs) file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if crsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "bzsfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if bzs file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition timeseries file not found: {abs_file_path}" + ) + + # Read bzs file (this creates one DataFrame with all timeseries) + df = utils.read_timeseries(abs_file_path, tref=self.model.config.get("tref")) + + # Now we need to split the timeseries into the different points + for idx, row in self.data.iterrows(): + # Get the timeseries for this point + ts = pd.DataFrame(df.iloc[:, idx]) + # Set the column name to wl + ts.columns = ["wl"] + # # Set the index to time + # ts.index.name = "time" + # Add to the point + self.data.at[idx, "timeseries"] = ts + + def read_boundary_conditions_astro(self, filename: str | Path = None): + """Read SFINCS boundary condition astro (*.bca) file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if bcafile is not None + abs_file_path = self.model.config.get_set_file_variable( + "bcafile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if bca file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition astro file not found: {abs_file_path}" + ) + + # Read bca file, which is actually some sort of toml file + d = IniStruct(filename=abs_file_path) + # Loop through boundary points + for ip, point in self.data.iterrows(): + # Set data in row of gdf + self.data.at[ip, "astro"] = d.section[ip].data + + def read_boundary_conditions_netcdf(self, filename: str | Path = None): + """Read SFINCS boundary conditions netcdf file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if netbndbzsbzifile is not None + abs_file_path = self.model.config.get_set_file_variable( + "netbndbzsbzifile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if netbndbzsbzifile exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition netcdf file not found: {abs_file_path}" + ) + + # Read netcdf file + ds = xr.open_dataset(abs_file_path) + + # Loop through boundary points + # FIXME - we first need to get the points! + for ip, point in self.data.iterrows(): + # Get the timeseries for this point + ts = ds["timeseries"].sel(point=ip).to_dataframe() + # Add to the point + self.data.at[ip, "timeseries"] = ts + + # Get the astro for this point + astro = ds["astro"].sel(point=ip).to_dataframe() + # Add to the point + self.data.at[ip, "astro"] = astro + + def write(self, format: str = None): + """Write SFINCS boundary conditions (*.bnd, *.bzs, *.bca files) or netcdf file. + + The format of the boundary conditions files can be specified, + otherwise it is determined from the model configuration. + + Parameters + ---------- + format : str, optional + Format of the boundary conditions files, "asc" (default), or "netcdf". + """ + + if self.data.empty: + # There are no boundary points + return + + if format is None: + if self.model.config.get("netbndbzsbzifile"): + format = "netcdf" + else: + format = "asc" + + if format == "asc": + self.write_boundary_points() + self.write_boundary_conditions_timeseries() + if self.model.config.get("bcafile"): + self.write_boundary_conditions_astro() + else: + self.write_boundary_conditions_netcdf() + + def write_boundary_points(self, filename: str | Path = None): + """Write SFINCS boundary condition points (*.bnd) file""" + + # Check that write mode is on + self.root._assert_write_mode() + + # Get absolute file name and set it in config if bndfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "bndfile", value=filename, default="sfincs.bnd" + ) + + # Write bnd file + # Change precision of coordinates according to crs + if self.model.crs.is_geographic: + fmt = "%11.6f" + else: + fmt = "%11.1f" + utils.write_xy(abs_file_path, self.data, fmt=fmt) + + def write_boundary_conditions_timeseries(self, filename: str | Path = None): + """Write SFINCS boundary condition timeseries (*.bzs) file""" + + # Check that write mode is on + self.root._assert_write_mode() + + # Get absolute file name and set it in config if bzsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "bzsfile", value=filename, default="sfincs.bzs" + ) + + # Get all timeseries and stick in one DataFrame + df = pd.DataFrame() + for ip, point in self.data.iterrows(): + df = pd.concat([df, point["timeseries"]["wl"]], axis=1) + + # Write to file + # This does NOT work at the moment! + # utils.write_timeseries(abs_file_path, df, self.model.config.get("tref")) + + # For now use 'ugly' to_csv method without control of column width + # Convert time index to datetime64 + time = pd.to_datetime(df.index) + tref = self.model.config.get("tref") + time = (time - tref).total_seconds() + df.index = time + df.to_csv( + abs_file_path, index=True, sep=" ", header=False, float_format="%0.3f" + ) + + # to_fwf(df, abs_file_path) + + def write_boundary_conditions_astro(self, filename: str | Path = None): + """Write SFINCS boundary condition astro (*.bca) file""" + + # Check that write mode is on + self.root._assert_write_mode() + + # Get absolute file name and set it in config if bcafile is not None + abs_file_path = self.model.config.get_set_file_variable( + "bcafile", value=filename, default="sfincs.bca" + ) + + # Write bca file + # Create IniStruct + d = IniStruct() + # Loop through boundary points + for ip, point in self.data.iterrows(): + # Add data to IniStruct + d.section[ip].data = point["astro"] + # Write to file + d.write(abs_file_path) + + def set(self, gdf: gpd.GeoDataFrame, merge: bool = True): + """Set boundary conditions data. + + Parameters + ---------- + gdf : gpd.GeoDataFrame + GeoDataFrame with boundary points. + merge : bool, optional + Merge data with existing data, by default True. + """ + + if merge: + self.data = pd.concat([self.data, gdf], ignore_index=True) + else: + self.data = gdf + + def add_point( + self, + gdf: gpd.GeoDataFrame = None, + x: float = None, + y: float = None, + wl: float = 0.0, + ): + """Add a single point to the boundary conditions data. Either gdf, + or x, y must be provided. + + Parameters + ---------- + gdf : gpd.GeoDataFrame + GeoDataFrame with a single point + x : float + x-coordinate of the point + y : float + y-coordinate of the point + wl : float + Water level of the point + """ + if gdf is not None: + if len(gdf) != 1: + raise ValueError( + "Only GeoDataFrame with a single point in a can be added." + ) + gdf = gdf.to_crs(self.model.crs) + if "timeseries" not in gdf: + gdf["timeseries"] = pd.DataFrame() + if "astro" not in gdf: + gdf["astro"] = pd.DataFrame() + else: + # Create a GeoDataFrame with a single point + if x is None or y is None: + raise ValueError("Either gdf or x, y, and name must be provided.") + point = shapely.geometry.Point(x, y) + gdf = gpd.GeoDataFrame( + [ + { + "timeseries": pd.DataFrame(), + "astro": pd.DataFrame(), + "geometry": point, + } + ], + crs=self.model.crs, + ) + + # Check if there is data in the timeseries + if gdf["timeseries"][0].empty: + # Now add the water level + if not self.data.empty: + # Set water level at same times as first existing point by copying + gdf.at[0, "timeseries"] = self.data.iloc[0]["timeseries"].copy() + gdf.at[0, "timeseries"]["wl"] = wl + else: + # First point, so need to generate df with constant water level + time = [self.model.config.get("tstart"), self.model.config.get("tstop")] + wl = [wl] * 2 + # Create DataFrame with columns time and wl + df = pd.DataFrame() + df["time"] = time + df["wl"] = wl + df = df.set_index("time") + gdf.at[0, "timeseries"] = df + else: + # Check if the timeseries is the same length as the first point + if len(gdf["timeseries"][0]) != len(self.data.iloc[0]["timeseries"]): + raise ValueError( + "Timeseries in gdf must be the same length as the first point in the boundary conditions data." + ) + + # Add to self.data + self.data = pd.concat([self.data, gdf], ignore_index=True) + + def delete(self, index: Union[int, List[int]]): + """Delete a single point from the boundary conditions data. + + Parameters + ---------- + index : int or list of int + Index or list of indices of points to be deleted. + """ + + if self.data.empty: + return + + if not isinstance(index, list): + index = [index] + # Check if indices are within range + if any(x > (len(self.data.index) - 1) for x in index): + raise ValueError("One of the indices exceeds length of index range!") + self.data = self.data.drop(index).reset_index(drop=True) + + if self.data.empty: + self.model.config.set("bndfile", None) + self.model.config.set("bzsfile", None) + self.model.config.set("bcafile", None) + self.model.config.set("netbndbzsbzifile", None) + + def clear(self): + """Clean GeoDataFrame with boundary points.""" + self.data = gpd.GeoDataFrame() + + def set_timeseries( + self, + index: Union[int, List[int]] = None, + shape: str = "constant", + timestep: float = 600.0, + offset: float = 0.0, + amplitude: float = 1.0, + phase: float = 0.0, + period: float = 43200.0, + peak: float = 1.0, + tpeak: float = 86400.0, + duration: float = 43200.0, + ): + """Applies time series boundary conditions for each point + Create numpy datetime64 array for time series with python datetime.datetime objects + + Parameters + ---------- + shape : str + Shape of the time series. Options are "constant", "sine", "gaussian", "astronomical". + timestep : float + Time step [s] + offset : float + Offset of the time series [m] + amplitude : float + Amplitude of the sine wave [m] + phase : float + Phase of the sine wave [degrees] + period : float + Period of the sine wave [s] + peak : float + Peak of the Gaussian wave [m] + tpeak : float + Time of the peak of the Gaussian wave [s] + duration : float + Duration of the Gaussian wave [s] + """ + + if self.data.empty: + return + + if shape == "astronomical": + # Use existing method + self.generate_bzs_from_bca(dt=timestep, offset=offset, write=False) + return + + t0 = np.datetime64(self.model.config.get("tstart")) + t1 = np.datetime64(self.model.config.get("tstop")) + if shape == "constant": + dt = np.timedelta64(int((t1 - t0).astype(float) / 1e6), "s") + else: + dt = np.timedelta64(int(timestep), "s") + time = np.arange(t0, t1 + dt, dt) + dtsec = dt.astype(float) + # Convert time to seconds since tref + tsec = ( + (time - np.datetime64(self.model.config.get("tref"))) + .astype("timedelta64[s]") + .astype(float) + ) + nt = len(tsec) + if shape == "constant": + wl = [offset] * nt + elif shape == "sine": + wl = offset + amplitude * np.sin( + 2 * np.pi * tsec / period + phase * np.pi / 180 + ) + elif shape == "gaussian": + wl = offset + peak * np.exp(-(((tsec - tpeak) / (0.25 * duration)) ** 2)) + elif shape == "astronomical": + # Not implemented + return + + times = pd.date_range( + start=t0, end=t1, freq=pd.tseries.offsets.DateOffset(seconds=dtsec) + ) + + if index is None: + index = list(self.data.index) + elif not isinstance(index, list): + index = [index] + + for i in index: + df = pd.DataFrame() + df["time"] = times + df["wl"] = wl + df = df.set_index("time") + self.data.at[i, "timeseries"] = df + + def generate_bzs_from_bca( + self, dt: float = 600.0, offset: float = 0.0, write_file: bool = True + ): + """Generate bzs file from bca file""" + + if self.data.empty: + return + + if not self.model.input.variables.bzsfile: + self.model.input.variables.bzsfile = "sfincs.bzs" + + times = pd.date_range( + start=self.model.input.variables.tstart, + end=self.model.input.variables.tstop, + freq=pd.tseries.offsets.DateOffset(seconds=dt), + ) + + # Make boundary conditions based on bca file + for icol, point in self.gdf.iterrows(): + v = predict(point.astro, times) + offset + ts = pd.Series(v, index=times) + # Convert this pandas series to a DataFrame + df = pd.DataFrame() + df["time"] = ts.index + df["wl"] = ts.values + df = df.set_index("time") + self.gdf.at[icol, "timeseries"] = df + + if write_file: + self.write_boundary_conditions_timeseries() + + def get_boundary_points_from_mask(self, min_dist=None, bnd_dist=5000.0): + # Should move this to mask? + + if min_dist is None: + # Set minimum distance between to grid boundary points on polyline to 2 * dx + min_dist = self.model.quadtree_grid.data.attrs["dx"] * 2 + + mask = self.model.quadtree_grid.data["mask"] + ibnd = np.where(mask == 2) + xz, yz = self.model.quadtree_grid.face_coordinates() + xp = xz[ibnd] + yp = yz[ibnd] + + # Make boolean array for points that are include in a polyline + used = np.full(xp.shape, False, dtype=bool) + + # Make list of polylines. Each polyline is a list of indices of boundary points. + polylines = [] + + while True: + if np.all(used): + # All boundary grid points have been used. We can stop now. + break + + # Find first the unused points + i1 = np.where(~used)[0][0] + + # Set this point to used + used[i1] = True + + # Start new polyline with index i1 + polyline = [i1] + + while True: + # Compute distances to all points that have not been used + xpunused = xp[~used] + ypunused = yp[~used] + # Get all indices of unused points + unused_indices = np.where(~used)[0] + + dst = np.sqrt((xpunused - xp[i1]) ** 2 + (ypunused - yp[i1]) ** 2) + if np.all(np.isnan(dst)): + break + inear = np.nanargmin(dst) + inearall = unused_indices[inear] + if dst[inear] < min_dist: + # Found next point along polyline + polyline.append(inearall) + used[inearall] = True + i1 = inearall + else: + # Last point found + break + + # Now work the other way + # Start with first point of polyline + i1 = polyline[0] + while True: + if np.all(used): + # All boundary grid points have been used. We can stop now. + break + # Now we go in the other direction + xpunused = xp[~used] + ypunused = yp[~used] + unused_indices = np.where(~used)[0] + dst = np.sqrt((xpunused - xp[i1]) ** 2 + (ypunused - yp[i1]) ** 2) + inear = np.nanargmin(dst) + inearall = unused_indices[inear] + if dst[inear] < min_dist: + # Found next point along polyline + polyline.insert(0, inearall) + used[inearall] = True + # Set index of next point + i1 = inearall + else: + # Last nearby point found + break + + if len(polyline) > 1: + polylines.append(polyline) + + gdf_list = [] + ip = 0 + # Transform to web mercator to get distance in metres + if self.model.crs.is_geographic: + transformer = Transformer.from_crs(self.model.crs, 3857, always_xy=True) + # Loop through polylines + for polyline in polylines: + x = xp[polyline] + y = yp[polyline] + points = [(x, y) for x, y in zip(x.ravel(), y.ravel())] + line = shapely.geometry.LineString(points) + if self.model.crs.is_geographic: + # Line in web mercator (to get length in metres) + xm, ym = transformer.transform(x, y) + pointsm = [(xm, ym) for xm, ym in zip(xm.ravel(), ym.ravel())] + linem = shapely.geometry.LineString(pointsm) + num_points = int(linem.length / bnd_dist) + 2 + else: + num_points = int(line.length / bnd_dist) + 2 + # Interpolate to new points + new_points = [ + line.interpolate(i / float(num_points - 1), normalized=True) + for i in range(num_points) + ] + # Loop through points in polyline + for point in new_points: + name = str(ip + 1).zfill(4) + d = { + "name": name, + "timeseries": pd.DataFrame(), + "astro": pd.DataFrame(), + "geometry": point, + } + gdf_list.append(d) + ip += 1 + + self.data = gpd.GeoDataFrame(gdf_list, crs=self.model.crs) + + +# def to_fwf(df, fname, floatfmt=".3f"): +# indx = df.index.tolist() +# vals = df.values.tolist() +# for it, t in enumerate(vals): +# t.insert(0, indx[it]) +# content = tabulate(vals, [], tablefmt="plain", floatfmt=floatfmt) +# open(fname, "w").write(content) diff --git a/hydromt_sfincs/config.py b/hydromt_sfincs/config.py index fabd523d..67f42663 100644 --- a/hydromt_sfincs/config.py +++ b/hydromt_sfincs/config.py @@ -1,4 +1,5 @@ from datetime import datetime +import time from typing import TYPE_CHECKING, List, Optional, Dict, Any from ast import literal_eval from os.path import abspath, isabs, join, split, exists @@ -6,7 +7,8 @@ from hydromt.model.components import ModelComponent -from hydromt_sfincs.config_variables import SfincsConfigVariables +# from hydromt_sfincs.config_variables import SfincsConfigVariables +from hydromt_sfincs.config_variables import sfincs_config_variables if TYPE_CHECKING: from hydromt_sfincs import SfincsModel @@ -17,14 +19,14 @@ class SfincsConfig(ModelComponent): def __init__(self, model: "SfincsModel"): self._filename = "sfincs.inp" - self._data: SfincsConfigVariables = None + self._data: sfincs_config_variables = None super().__init__(model=model) @property - def data(self) -> SfincsConfigVariables: + def data(self): """Return the SfincsConfig object.""" if self._data is None: - self._data = SfincsConfigVariables() + self._data = sfincs_config_variables return self._data def read(self, filename: str = "sfincs.inp") -> None: @@ -75,7 +77,7 @@ def read(self, filename: str = "sfincs.inp") -> None: inp_dict[name] = val # Convert dictionary to SfincsConfig instance - self._data = SfincsConfigVariables(**inp_dict) + self._data = self._data.copy(update=inp_dict) # Update the grid properties from the configuration # This will either drop the quadtree component or the regular component? @@ -126,21 +128,23 @@ def get(self, key: str, fallback: Any = None, abs_path: bool = False) -> Any: return value - def set(self, key: str, value: Any) -> None: + def set(self, key: str, value: Any, skip_validation=False) -> None: """Set a value with validation using Pydantic's model_copy.""" + if not hasattr(self.data, key): raise KeyError(f"'{key}' is not a valid attribute of SfincsConfig.") - # Validate the new data - # FIXME implement this in a better way - try: - # If key was an extra field, do NOT set it here - if key in self.data.model_fields: - value = SfincsConfigVariables(**{key: value}).__dict__[key] - except Exception as e: - raise TypeError(f"Invalid input type for '{key}'") + if not skip_validation: + # Validate the new data + # FIXME implement this in a better way + # It works, but it is quite slow when all the variables are set in a loop + # Therefore the skip_validation option is added + try: + self.data.model_validate({key: value}) + except Exception as e: + raise TypeError(f"Invalid input type for '{key}'") - self._data = self._data.model_copy(update={key: value}) + self.data.__setattr__(key, value) def update(self, dict: Dict[str, Any]) -> None: """ @@ -167,6 +171,7 @@ def update_grid_from_config(self) -> None: self.model.grid.update_grid_from_config() # drop quadtree component self.model.components.pop("quadtree", None) + # TODO also drop mask and subgrid components? elif self.model.grid_type == "quadtree": # drop regular component self.model.components.pop("grid", None) diff --git a/hydromt_sfincs/config_variables.py b/hydromt_sfincs/config_variables.py index 0f00b396..dd82b254 100644 --- a/hydromt_sfincs/config_variables.py +++ b/hydromt_sfincs/config_variables.py @@ -618,6 +618,26 @@ class Config: le=1, description="Option in integrated SnapWave solver to turn on wind growth process (1: yes, 0: no)", ) + snapwave_bndfile: str | None = Field( + None, + description="Name of the SnapWave boundary points file", + ) + snapwave_bhsfile: str | None = Field( + None, + description="Name of the SnapWave wave height time-series file", + ) + snapwave_btpfile: str | None = Field( + None, + description="Name of the SnapWave wave period time-series file", + ) + snapwave_bwdfile: str | None = Field( + None, + description="Name of the SnapWave wave direction time-series file", + ) + snapwave_bdsfile: str | None = Field( + None, + description="Name of the SnapWave wave spreading time-series file", + ) # # Wind drag # @@ -637,3 +657,6 @@ class Config: # bcafile: str | None = Field(None, description="Name of the calibration file") corfile: str | None = Field(None, description="Name of the correction file") + + +sfincs_config_variables = SfincsConfigVariables() diff --git a/hydromt_sfincs/discharge_points.py b/hydromt_sfincs/discharge_points.py index b5a1c6fa..7d34a936 100644 --- a/hydromt_sfincs/discharge_points.py +++ b/hydromt_sfincs/discharge_points.py @@ -4,6 +4,10 @@ import xarray as xr from pathlib import Path from typing import Union, List +import shapely +from pyproj import Transformer + +# from tabulate import tabulate from hydromt.model.components import ModelComponent from hydromt.model import Model @@ -15,60 +19,465 @@ def __init__( self, model: Model, ): - self._filename: str = "sfincs.src" - self._data: gpd.GeoDataFrame = None + # self._filename: str = "sfincs.dis" # FIXME - List(str = "sfincs.dis" and str = "sfincs.src" or str = "sfincs_netbndbzsbzi.nc") + self.data = gpd.GeoDataFrame() super().__init__( model=model, ) - @property - def data(self) -> gpd.GeoDataFrame: - """Discharge points data. + # @property + # def data(self) -> gpd.GeoDataFrame: + # """Water level discharge conditions data. + + # Return pd.GeoDataFrame + # """ + # if self._data is None: + # self._initialize() + # return self._data + + # def _initialize(self) -> None: + # """Initialize discharge conditions data.""" + # if self._data is None: + # self._data = gpd.GeoDataFrame() + + def read(self, format: str = None): + """Read SFINCS discharge points (*.dis, *.src files) or netcdf file. + + The format of the discharge conditions files can be specified, + otherwise it is determined from the model configuration. + + Parameters + ---------- + format : str, optional + Format of the discharge files, "asc" or "netcdf". + """ + + if format is None: + if self.model.config.get("netsrcdisfile"): + format = "netcdf" + else: + format = "asc" + + if format == "asc": + self.read_discharge_points() + # Check if there are any points + if not self.data.empty: + self.read_discharge_timeseries() + elif format == "netcdf": + # Read netcdf file + self.read_discharges_netcdf() + + def read_discharge_points(self, filename: str | Path = None): + """Read SFINCS discharge points (*.src) file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if crsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "srcfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if src file exists + if not abs_file_path.exists(): + raise FileNotFoundError(f"Discharge points file not found: {abs_file_path}") + + # HydroMT does not have open_vector at the moment ... + # Read bnd file + # gdf = utils.read_xy(abs_file_path, crs=self.model.crs) + # # Add columns for timeseries and astro and add empty DataFrames + # gdf["timeseries"] = pd.DataFrame() + # gdf["astro"] = pd.DataFrame() + # # Add to self.data + # self.data = gdf + + # Read the bnd file + df = pd.read_csv( + abs_file_path, index_col=False, header=None, names=["x", "y"], sep="\s+" + ) + + gdf_list = [] + # Loop through points + for ind in range(len(df.x.values)): + name = str(ind + 1).zfill(4) + x = df.x.values[ind] + y = df.y.values[ind] + point = shapely.geometry.Point(x, y) + d = { + "name": name, + "timeseries": pd.DataFrame(), + "geometry": point, + } + gdf_list.append(d) + self.data = gpd.GeoDataFrame(gdf_list, crs=self.model.crs) + + def read_discharge_timeseries(self, filename: str | Path = None): + """Read SFINCS discharge condition timeseries (*.bzs) file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if crsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "disfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if dis file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Discharge timeseries file not found: {abs_file_path}" + ) + + # Read bzs file (this creates one DataFrame with all timeseries) + df = utils.read_timeseries(abs_file_path, tref=self.model.config.get("tref")) + + # Now we need to split the timeseries into the different points + for idx, row in self.data.iterrows(): + # Get the timeseries for this point + ts = pd.DataFrame(df.iloc[:, idx]) + # Set the column name to wl + ts.columns = ["q"] + # # Set the index to time + # ts.index.name = "time" + # Add to the point + self.data.at[idx, "timeseries"] = ts + + def read_discharge_conditions_netcdf(self, filename: str | Path = None): + """Read SFINCS discharge conditions netcdf file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if netbndbzsbzifile is not None + abs_file_path = self.model.config.get_set_file_variable( + "netbndbzsbzifile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if netbndbzsbzifile exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"discharge condition netcdf file not found: {abs_file_path}" + ) + + # Read netcdf file + ds = xr.open_dataset(abs_file_path) + + # Loop through discharge points + # FIXME - we first need to get the points! + for ip, point in self.data.iterrows(): + # Get the timeseries for this point + ts = ds["timeseries"].sel(point=ip).to_dataframe() + # Add to the point + self.data.at[ip, "timeseries"] = ts + + # Get the astro for this point + astro = ds["astro"].sel(point=ip).to_dataframe() + # Add to the point + self.data.at[ip, "astro"] = astro + + def write(self, format: str = None): + """Write SFINCS discharges (*.src, *.dis files) or netcdf file. + + The format of the discharge files can be specified, + otherwise it is determined from the model configuration. + + Parameters + ---------- + format : str, optional + Format of the discharge files, "asc" (default), or "netcdf". + """ + + if self.data.empty: + # There are no discharge points + return + + if format is None: + if self.model.config.get("netsrcdisfile"): + format = "netcdf" + else: + format = "asc" + + if format == "asc": + self.write_discharge_points() + self.write_discharge_timeseries() + else: + self.write_discharges_netcdf() + + def write_discharge_points(self, filename: str | Path = None): + """Write SFINCS discharge points (*.src) file""" + + # Check that write mode is on + self.root._assert_write_mode() + + # Get absolute file name and set it in config if bndfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "srcfile", value=filename, default="sfincs.src" + ) + + # Write src file + # Change precision of coordinates according to crs + if self.model.crs.is_geographic: + fmt = "%11.6f" + else: + fmt = "%11.1f" + + utils.write_xyn(abs_file_path, self.data, fmt=fmt) + + def write_discharge_timeseries(self, filename: str | Path = None): + """Write SFINCS discharge timeseries (*.dis) file""" + + # Check that write mode is on + self.root._assert_write_mode() + + # Get absolute file name and set it in config if bzsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "disfile", value=filename, default="sfincs.dis" + ) + + # Get all timeseries and stick in one DataFrame + df = pd.DataFrame() + for ip, point in self.data.iterrows(): + df = pd.concat([df, point["timeseries"]["q"]], axis=1) + + # Write to file + # This does NOT work at the moment! + # utils.write_timeseries(abs_file_path, df, self.model.config.get("tref")) + + # For now use 'ugly' to_csv method without control of column width + # Convert time index to datetime64 + time = pd.to_datetime(df.index) + tref = self.model.config.get("tref") + time = (time - tref).total_seconds() + df.index = time + df.to_csv( + abs_file_path, index=True, sep=" ", header=False, float_format="%0.3f" + ) + + # to_fwf(df, abs_file_path) + + def set(self, gdf: gpd.GeoDataFrame, merge: bool = True): + """Set discharge data. + + Parameters + ---------- + gdf : gpd.GeoDataFrame + GeoDataFrame with discharge points. + merge : bool, optional + Merge data with existing data, by default True. + """ + + if merge: + self.data = pd.concat([self.data, gdf], ignore_index=True) + else: + self.data = gdf + + def add_point( + self, + gdf: gpd.GeoDataFrame = None, + name: str = None, + x: float = None, + y: float = None, + q: float = 0.0, + ): + """Add a single point to the discharge data. Either gdf, + or x, y must be provided. + + Parameters + ---------- + gdf : gpd.GeoDataFrame + GeoDataFrame with a single point + name : str + Name of the point + x : float + x-coordinate of the point + y : float + y-coordinate of the point + wl : float + Water level of the point + """ + if gdf is not None: + if len(gdf) != 1: + raise ValueError( + "Only GeoDataFrame with a single point in a can be added." + ) + gdf = gdf.to_crs(self.model.crs) + if "timeseries" not in gdf: + gdf["timeseries"] = pd.DataFrame() + else: + # Create a GeoDataFrame with a single point + if x is None or y is None or name is None: + raise ValueError("Either gdf or x, y, and name must be provided.") + point = shapely.geometry.Point(x, y) + gdf = gpd.GeoDataFrame( + [ + { + "name": name, + "timeseries": pd.DataFrame(), + "geometry": point, + } + ], + crs=self.model.crs, + ) + + # Check if there is data in the timeseries + if gdf["timeseries"][0].empty: + # Now add the water level + if not self.data.empty: + # Set water level at same times as first existing point by copying + gdf.at[0, "timeseries"] = self.data.iloc[0]["timeseries"].copy() + gdf.at[0, "timeseries"]["q"] = q + else: + # First point, so need to generate df with constant water level + time = [self.model.config.get("tstart"), self.model.config.get("tstop")] + q = [q] * 2 + # Create DataFrame with columns time and wl + df = pd.DataFrame() + df["time"] = time + df["q"] = q + df = df.set_index("time") + gdf.at[0, "timeseries"] = df + else: + # Check if the timeseries is the same length as the first point + if len(gdf["timeseries"][0]) != len(self.data.iloc[0]["timeseries"]): + raise ValueError( + "Timeseries in gdf must be the same length as the first point in the discharge conditions data." + ) + + # Add to self.data + self.data = pd.concat([self.data, gdf], ignore_index=True) + + def delete(self, index: Union[int, List[int]]): + """Delete a single point from the discharge data. + + Parameters + ---------- + index : int or list of int + Index or list of indices of points to be deleted. + """ + + if self.data.empty: + return + + if not isinstance(index, list): + index = [index] + # Check if indices are within range + if any(x > (len(self.data.index) - 1) for x in index): + raise ValueError("One of the indices exceeds length of index range!") + self.data = self.data.drop(index).reset_index(drop=True) + + if self.data.empty: + self.model.config.set("srcfile", None) + self.model.config.set("disfile", None) + self.model.config.set("netsrcdisfile", None) + + def clear(self): + """Clean GeoDataFrame with discharge points.""" + self.data = gpd.GeoDataFrame() + + def set_timeseries( + self, + index: Union[int, List[int]] = None, + shape: str = "constant", + timestep: float = 600.0, + offset: float = 0.0, + amplitude: float = 1.0, + phase: float = 0.0, + period: float = 43200.0, + peak: float = 1.0, + tpeak: float = 86400.0, + duration: float = 43200.0, + ): + """Applies time series discharges for each point + Create numpy datetime64 array for time series with python datetime.datetime objects - Return gpd.GeoDataFrame + Parameters + ---------- + shape : str + Shape of the time series. Options are "constant", "sine", or "gaussian". + timestep : float + Time step [s] + offset : float + Offset of the time series [m] + amplitude : float + Amplitude of the sine wave [m] + phase : float + Phase of the sine wave [degrees] + period : float + Period of the sine wave [s] + peak : float + Peak of the Gaussian wave [m] + tpeak : float + Time of the peak of the Gaussian wave [s] + duration : float + Duration of the Gaussian wave [s] """ - if self._data is None: - self._initialize() - return self._data - - # Original HydroMT-SFINCS setup_ functions: - # setup_discharge_forcing - # setup_discharge_forcing_from_grid - - -# %% core HydroMT-SFINCS functions: -# _initialize -# read -# read_discharge_netcdf - netsrcdisfile -# read_discharge_points - srcfile -# read_discharge_timeseries - disfile -# write -# write_snapwave_conditions - netsnapwavefile (or other default name) -# and/or: write_discharge_points - srcfile + write_discharge_conditions_timeseries - disfile> FIXME don't support? -# set -# set_points - gpd.GeoDataFrame with points -# set_timeseries - pd.DataFrame with timeseries (?) -# create: -# create -# create_from_grid -# add -# delete -# clear - -# %% DDB GUI focused additional functions: -# add_point (exception, call add for just one point) -# delete_point (exception, call delete for just one point) - -# to_gdf -# has_open_boundaries - -# get_boundary_points_from_mask -# set_timeseries -# set_timeseries_uniform -# set_conditions_at_point -# read_timeseries_file -# to_fwf -# inpolygon -# -# class SnapWaveBoundaryEnclosure > we are not going to support this > directly in snapwave_mask -# class SnapWaveMask > will be part of Mask class - exactly same functionality as for SFINCS mask + + if self.data.empty: + return + + t0 = np.datetime64(self.model.config.get("tstart")) + t1 = np.datetime64(self.model.config.get("tstop")) + if shape == "constant": + dt = np.timedelta64(int((t1 - t0).astype(float) / 1e6), "s") + else: + dt = np.timedelta64(int(timestep), "s") + time = np.arange(t0, t1 + dt, dt) + dtsec = dt.astype(float) + # Convert time to seconds since tref + tsec = ( + (time - np.datetime64(self.model.config.get("tref"))) + .astype("timedelta64[s]") + .astype(float) + ) + nt = len(tsec) + if shape == "constant": + q = [offset] * nt + elif shape == "sine": + q = offset + amplitude * np.sin( + 2 * np.pi * tsec / period + phase * np.pi / 180 + ) + elif shape == "gaussian": + q = offset + peak * np.exp(-(((tsec - tpeak) / (0.25 * duration)) ** 2)) + else: + # Not implemented + return + + times = pd.date_range( + start=t0, end=t1, freq=pd.tseries.offsets.DateOffset(seconds=dtsec) + ) + + if index is None: + index = list(self.data.index) + elif not isinstance(index, list): + index = [index] + + for i in index: + df = pd.DataFrame() + df["time"] = times + df["q"] = q + df = df.set_index("time") + self.data.at[i, "timeseries"] = df + + +# def to_fwf(df, fname, floatfmt=".3f"): +# indx = df.index.tolist() +# vals = df.values.tolist() +# for it, t in enumerate(vals): +# t.insert(0, indx[it]) +# content = tabulate(vals, [], tablefmt="plain", floatfmt=floatfmt) +# open(fname, "w").write(content) diff --git a/hydromt_sfincs/observation_points.py b/hydromt_sfincs/observation_points.py index a38ef527..0a68d25b 100644 --- a/hydromt_sfincs/observation_points.py +++ b/hydromt_sfincs/observation_points.py @@ -105,11 +105,11 @@ def write(self, filename=None): key="obsfile", value=filename, default_filename="sfincs.obs" ) - # change precision of coordinates according to crs + # Change precision of coordinates according to crs if self.model.crs.is_geographic: - fmt = "%.6f" + fmt = "%11.6f" else: - fmt = "%.1f" + fmt = "%11.1f" utils.write_xyn(file_path, self.data, fmt=fmt) # =utils.py function diff --git a/hydromt_sfincs/quadtree.py b/hydromt_sfincs/quadtree.py index d95da23c..5fc0ef5d 100644 --- a/hydromt_sfincs/quadtree.py +++ b/hydromt_sfincs/quadtree.py @@ -124,7 +124,8 @@ def read( if not isabs(filename): self._filename = join(self.root.path, filename) - dsu = xu_open_dataset(self._filename) + # dsu = xu_open_dataset(self._filename) + dsu = xu.load_dataset(self._filename) # set CRS (not sure if that should be stored in the netcdf in this way) dsu.grid.set_crs(CRS.from_wkt(dsu["crs"].crs_wkt)) diff --git a/hydromt_sfincs/quadtree_mask.py b/hydromt_sfincs/quadtree_mask.py index a2d618c1..df53e61c 100644 --- a/hydromt_sfincs/quadtree_mask.py +++ b/hydromt_sfincs/quadtree_mask.py @@ -65,7 +65,7 @@ def build( print("Building mask ...") mask = np.zeros(self.model.quadtree_grid.nr_cells, dtype=np.int8) - x, y = self.model.quadtree_grid.face_coordinates() + x, y = self.model.quadtree_grid.face_coordinates z = self.model.quadtree_grid.data["z"].values[:] # Indices are 1-based in SFINCS so subtract 1 for python 0-based indexing diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py index 432c69a8..4b1bbd43 100644 --- a/hydromt_sfincs/sfincs.py +++ b/hydromt_sfincs/sfincs.py @@ -30,24 +30,25 @@ # %% Import model components from hydromt.model import Model -# input types: +# Input types from hydromt_sfincs.config import SfincsConfig -# grid types: -from hydromt_sfincs.quadtree import QuadtreeGrid -from hydromt_sfincs.quadtree_mask import QuadtreeMask +# Grid types from hydromt_sfincs.regulargrid import RegularGrid -from hydromt_sfincs.subgrid import SubgridTableRegular +from hydromt_sfincs.quadtree import QuadtreeGrid -# map types: +# Map types # from hydromt_sfincs.mask import SfincsMask +from hydromt_sfincs.quadtree_mask import QuadtreeMask + # from hydromt_sfincs.bathymetry import SfincsBathymetry +from hydromt_sfincs.subgrid import SubgridTableRegular from hydromt_sfincs.infiltration import SfincsInfiltration from hydromt_sfincs.manning_roughness import SfincsManningRoughness from hydromt_sfincs.initial_conditions import SfincsInitialConditions from hydromt_sfincs.storage_volume import SfincsStorageVolume -# geoms types: +# Geoms types from hydromt_sfincs.observation_points import SfincsObservationPoints from hydromt_sfincs.cross_sections import SfincsCrossSections from hydromt_sfincs.weirs import SfincsWeirs @@ -56,12 +57,10 @@ from hydromt_sfincs.drainage_structures import SfincsDrainageStructures from hydromt_sfincs.rivers import SfincsRivers -# forcing types: +# Forcing types +from hydromt_sfincs.boundary_conditions import SfincsBoundaryConditions from hydromt_sfincs.discharge_points import SfincsDischargePoints - -# from hydromt_sfincs.boundary_conditions import SfincsBoundaryConditions #/ -from hydromt_sfincs.waterlevel_conditions import SfincsWaterlevelConditions -from hydromt_sfincs.snapwave_conditions import SfincsSnapWaveConditions +from hydromt_sfincs.snapwave_boundary_conditions import SfincsSnapWaveBoundaryConditions # from hydromt_sfincs.meteo import SfincsMeteo from hydromt_sfincs.meteo import SfincsPrecipitation, SfincsPressure, SfincsWind @@ -126,34 +125,37 @@ def __init__( self.grid_type = "regular" self.add_component("config", SfincsConfig(self)) - # grid types: + + # Grid types self.add_component("grid", RegularGrid(self)) self.add_component("quadtree_grid", QuadtreeGrid(self)) - self.add_component("quadtree_mask", QuadtreeMask(self)) - # self.add_component("subgrid", SubgridTableRegular(self)) - # self.add_component("subgrid", SubgridTableRegular(self)) - # map types: + # Map types + self.add_component("quadtree_mask", QuadtreeMask(self)) # self.add_component("mask", SfincsMask(self)) # self.add_component("bathymetry", SfincsBathymetry(self)) # self.add_component("infiltration", SfincsInfiltration(self)) # self.add_component("manning_roughness", SfincsManningRoughness(self)) # self.add_component("initial_conditions", SfincsInitialConditions(self)) # self.add_component("storage_volume", SfincsStorageVolume(self)) + # self.add_component("subgrid", SubgridTableRegular(self)) + # self.add_component("quadtree_subgrid", SubgridTableQuadtree(self)) - # geoms types: + # Geoms types self.add_component("observation_points", SfincsObservationPoints(self)) self.add_component("cross_sections", SfincsCrossSections(self)) self.add_component("thin_dams", SfincsThinDams(self)) # self.add_component("weirs", SfincsWeirs(self)) - # self.add_component("wave_makers", SfincsWaveMakers(self)) + self.add_component("wave_makers", SfincsWaveMakers(self)) # self.add_component("drainage_structures", SfincsDrainageStructures(self)) # self.add_component("rivers", SfincsRivers(self)) - # forcing types: - # self.add_component("discharge_points", SfincsDischargePoints(self)) - # self.add_component("waterlevel_conditions", SfincsWaterlevelConditions(self)) - # self.add_component("snapwave_conditions", SfincsSnapWaveConditions(self)) + # Forcing types + self.add_component("boundary_conditions", SfincsBoundaryConditions(self)) + self.add_component("discharge_points", SfincsDischargePoints(self)) + self.add_component( + "snapwave_boundary_conditions", SfincsSnapWaveBoundaryConditions(self) + ) # self.add_component("meteo", SfincsMeteo(self)) # self.add_component("precipitation", SfincsPrecipitation(self)) # self.add_component("pressure", SfincsPressure(self)) diff --git a/hydromt_sfincs/snapwave_boundary_conditions.py b/hydromt_sfincs/snapwave_boundary_conditions.py new file mode 100644 index 00000000..e3ea106d --- /dev/null +++ b/hydromt_sfincs/snapwave_boundary_conditions.py @@ -0,0 +1,769 @@ +import geopandas as gpd +import numpy as np +import pandas as pd +import xarray as xr +from pathlib import Path +from typing import Union, List +import shapely +from pyproj import Transformer + +# from tabulate import tabulate + +from hydromt.model.components import ModelComponent +from hydromt.model import Model +from hydromt_sfincs import utils + + +class SfincsSnapWaveBoundaryConditions(ModelComponent): + def __init__( + self, + model: Model, + ): + self.data = gpd.GeoDataFrame() + super().__init__( + model=model, + ) + + # @property + # def data(self) -> gpd.GeoDataFrame: + # """Water level boundary conditions data. + + # Return pd.GeoDataFrame + # """ + # if self._data is None: + # self._initialize() + # return self._data + + # def _initialize(self) -> None: + # """Initialize boundary conditions data.""" + # if self._data is None: + # self._data = gpd.GeoDataFrame() + + def read(self, format: str = None): + """Read SFINCS boundary conditions (*.bnd, *.bzs, *.bca files) or netcdf file. + + The format of the boundary conditions files can be specified, + otherwise it is determined from the model configuration. + + Parameters + ---------- + format : str, optional + Format of the boundary conditions files, "asc" or "netcdf". + """ + + if format is None: + if self.model.config.get("netbndbzsbzifile"): + format = "netcdf" + else: + format = "asc" + + if format == "asc": + self.read_boundary_points() + # Check if there are any points + if not self.data.empty: + self.read_boundary_conditions_timeseries() + # Read astro if bcafile is defined + if self.model.config.get("bcafile"): + self.read_boundary_conditions_astro() + elif format == "netcdf": + # Read netcdf file + self.read_boundary_conditions_netcdf() + + def read_boundary_points(self, filename: str | Path = None): + """Read SnapWave boundary condition points (*.bnd) file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if crsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bndfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if bnd file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition points file not found: {abs_file_path}" + ) + + # HydroMT does not have open_vector at the moment ... + # Read bnd file + # gdf = utils.read_xy(abs_file_path, crs=self.model.crs) + # # Add columns for timeseries and astro and add empty DataFrames + # gdf["timeseries"] = pd.DataFrame() + # gdf["astro"] = pd.DataFrame() + # # Add to self.data + # self.data = gdf + + # Read the bnd file + df = pd.read_csv( + abs_file_path, index_col=False, header=None, names=["x", "y"], sep="\s+" + ) + + gdf_list = [] + # Loop through points + for ind in range(len(df.x.values)): + name = str(ind + 1).zfill(4) + x = df.x.values[ind] + y = df.y.values[ind] + point = shapely.geometry.Point(x, y) + d = { + "name": name, + "timeseries": pd.DataFrame(), + "geometry": point, + } + gdf_list.append(d) + self.data = gpd.GeoDataFrame(gdf_list, crs=self.model.crs) + + def read_boundary_conditions_timeseries(self, filename: str | Path = None): + """Read SFINCS boundary condition timeseries (*.bhs, *.btp, *.bwd, *.bds) files""" + + # Check that read mode is on + self.root._assert_read_mode() + + ### BHS file + + # Get absolute file name and set it in config if snapwave_bhsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bhsfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if bzs file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition timeseries file not found: {abs_file_path}" + ) + + # Read bhs file (this creates one DataFrame with all timeseries) + dfhs = utils.read_timeseries(abs_file_path, tref=self.model.config.get("tref")) + + ### BTP file + + # Get absolute file name and set it in config if snapwave_bhsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_btpfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if bzs file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition timeseries file not found: {abs_file_path}" + ) + + # Read btp file (this creates one DataFrame with all timeseries) + dftp = utils.read_timeseries(abs_file_path, tref=self.model.config.get("tref")) + + ### BWD file + + # Get absolute file name and set it in config if snapwave_bwdfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bwdfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if snapwave_bwd file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition timeseries file not found: {abs_file_path}" + ) + + # Read bwd file (this creates one DataFrame with all timeseries) + dfwd = utils.read_timeseries(abs_file_path, tref=self.model.config.get("tref")) + + ### BDS file + + # Get absolute file name and set it in config if snapwave_bdsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bdsfile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if snapwave_bds file exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition timeseries file not found: {abs_file_path}" + ) + + # Read bds file (this creates one DataFrame with all timeseries) + dfds = utils.read_timeseries(abs_file_path, tref=self.model.config.get("tref")) + + # Now we need to split the timeseries into the different points + for idx, row in self.data.iterrows(): + # Get the timeseries for this point + tshs = pd.DataFrame(dfhs.iloc[:, idx]) + tstp = pd.DataFrame(dftp.iloc[:, idx]) + tswd = pd.DataFrame(dfwd.iloc[:, idx]) + tsds = pd.DataFrame(dfds.iloc[:, idx]) + # Set the column name to hs + tshs.columns = ["hs"] + tstp.columns = ["tp"] + tswd.columns = ["wd"] + tsds.columns = ["ds"] + # Concatenate the DataFrames + ts = pd.concat([tshs, tstp, tswd, tsds], axis=1) + # Add to the point + self.data.at[idx, "timeseries"] = ts + + def read_boundary_conditions_netcdf(self, filename: str | Path = None): + """Read SFINCS boundary conditions netcdf file""" + + # Check that read mode is on + self.root._assert_read_mode() + + # Get absolute file name and set it in config if netbndbzsbzifile is not None + abs_file_path = self.model.config.get_set_file_variable( + "netbndbzsbzifile", value=filename + ) + + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined + return + + # Check if netbndbzsbzifile exists + if not abs_file_path.exists(): + raise FileNotFoundError( + f"Boundary condition netcdf file not found: {abs_file_path}" + ) + + # Read netcdf file + ds = xr.open_dataset(abs_file_path) + + # Loop through boundary points + # FIXME - we first need to get the points! + for ip, point in self.data.iterrows(): + # Get the timeseries for this point + ts = ds["timeseries"].sel(point=ip).to_dataframe() + # Add to the point + self.data.at[ip, "timeseries"] = ts + + ds.close() + + def write(self, format: str = None): + """Write SnapWave boundary conditions (*.bnd, *.bhs, *.btp, *.bwd, *.bds files) or netcdf file. + + The format of the boundary conditions files can be specified, + otherwise it is determined from the model configuration. + + Parameters + ---------- + format : str, optional + Format of the boundary conditions files, "asc" (default), or "netcdf". + """ + + if self.data.empty: + # There are no boundary points + return + + if format is None: + if self.model.config.get("netbndsnapwavefile"): + format = "netcdf" + else: + format = "asc" + + if format == "asc": + self.write_boundary_points() + self.write_boundary_conditions_timeseries() + else: + self.write_boundary_conditions_netcdf() + + def write_boundary_points(self, filename: str | Path = None): + """Write SnapWave boundary condition points (*.bnd) file""" + + # Check that write mode is on + self.root._assert_write_mode() + + # Get absolute file name and set it in config if bndfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bndfile", value=filename, default="snapwave.bnd" + ) + + # Write bnd file + # Change precision of coordinates according to crs + if self.model.crs.is_geographic: + fmt = "%11.6f" + else: + fmt = "%11.1f" + utils.write_xy(abs_file_path, self.data, fmt=fmt) + + def write_boundary_conditions_timeseries(self, filename: str | Path = None): + """Write SnapWave boundary condition timeseries (*.bhs, *.btp, *.bwd, *.bds) file""" + + # Check that write mode is on + self.root._assert_write_mode() + + # BHS + + # Get absolute file name and set it in config if snapwave_bhsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bhsfile", value=filename, default="snapwave.bhs" + ) + + # Get all timeseries and stick in one DataFrame + df = pd.DataFrame() + for ip, point in self.data.iterrows(): + df = pd.concat([df, point["timeseries"]["hs"]], axis=1) + + # Write to file + # This does NOT work at the moment! + # utils.write_timeseries(abs_file_path, df, self.model.config.get("tref")) + # For now use 'ugly' to_csv method without control of column width + # Convert time index to datetime64 + time = pd.to_datetime(df.index) + tref = self.model.config.get("tref") + time = (time - tref).total_seconds() + df.index = time + df.to_csv( + abs_file_path, index=True, sep=" ", header=False, float_format="%0.3f" + ) + + # BTP + + # Get absolute file name and set it in config if snapwave_btpfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_btpfile", value=filename, default="snapwave.btp" + ) + + # Get all timeseries and stick in one DataFrame + df = pd.DataFrame() + for ip, point in self.data.iterrows(): + df = pd.concat([df, point["timeseries"]["tp"]], axis=1) + + # Write to file + # This does NOT work at the moment! + # utils.write_timeseries(abs_file_path, df, self.model.config.get("tref")) + # For now use 'ugly' to_csv method without control of column width + # Convert time index to datetime64 + time = pd.to_datetime(df.index) + tref = self.model.config.get("tref") + time = (time - tref).total_seconds() + df.index = time + df.to_csv( + abs_file_path, index=True, sep=" ", header=False, float_format="%0.3f" + ) + + # BWD + + # Get absolute file name and set it in config if snapwave_bwdfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bwdfile", value=filename, default="snapwave.bwd" + ) + + # Get all timeseries and stick in one DataFrame + df = pd.DataFrame() + for ip, point in self.data.iterrows(): + df = pd.concat([df, point["timeseries"]["wd"]], axis=1) + + # Write to file + # This does NOT work at the moment! + # utils.write_timeseries(abs_file_path, df, self.model.config.get("tref")) + # For now use 'ugly' to_csv method without control of column width + # Convert time index to datetime64 + time = pd.to_datetime(df.index) + tref = self.model.config.get("tref") + time = (time - tref).total_seconds() + df.index = time + df.to_csv( + abs_file_path, index=True, sep=" ", header=False, float_format="%0.3f" + ) + + # BDS + + # Get absolute file name and set it in config if snapwave_bdsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "snapwave_bdsfile", value=filename, default="snapwave.bds" + ) + + # Get all timeseries and stick in one DataFrame + # This does NOT work at the moment! + # utils.write_timeseries(abs_file_path, df, self.model.config.get("tref")) + # For now use 'ugly' to_csv method without control of column width + # Convert time index to datetime64 + time = pd.to_datetime(df.index) + tref = self.model.config.get("tref") + time = (time - tref).total_seconds() + df.index = time + df.to_csv( + abs_file_path, index=True, sep=" ", header=False, float_format="%0.3f" + ) + + def set(self, gdf: gpd.GeoDataFrame, merge: bool = True): + """Set SnapWave boundary conditions data. + + Parameters + ---------- + gdf : gpd.GeoDataFrame + GeoDataFrame with boundary points. + merge : bool, optional + Merge data with existing data, by default True. + """ + + if merge: + self.data = pd.concat([self.data, gdf], ignore_index=True) + else: + self.data = gdf + + def add_point( + self, + gdf: gpd.GeoDataFrame = None, + x: float = None, + y: float = None, + hs: float = 1.0, + tp: float = 10.0, + wd: float = 270.0, + ds: float = 20.0, + ): + """Add a single point to the boundary conditions data. Either gdf, + or x, y must be provided. + + Parameters + ---------- + gdf : gpd.GeoDataFrame + GeoDataFrame with a single point + x : float + x-coordinate of the point + y : float + y-coordinate of the point + hs : float + Wave height of the point + tp : float + Peak period of the point + wd : float + Wave direction of the point + ds : float + Directional spread of the point + """ + if gdf is not None: + if len(gdf) != 1: + raise ValueError( + "Only GeoDataFrame with a single point in a can be added." + ) + gdf = gdf.to_crs(self.model.crs) + if "timeseries" not in gdf: + gdf["timeseries"] = pd.DataFrame() + else: + # Create a GeoDataFrame with a single point + if x is None or y is None: + raise ValueError("Either gdf or x, y, and name must be provided.") + point = shapely.geometry.Point(x, y) + gdf = gpd.GeoDataFrame( + [ + { + "timeseries": pd.DataFrame(), + "geometry": point, + } + ], + crs=self.model.crs, + ) + + # Check if there is data in the timeseries + if gdf["timeseries"][0].empty: + # Now add the water level + if not self.data.empty: + # Set water level at same times as first existing point by copying + gdf.at[0, "timeseries"] = self.data.iloc[0]["timeseries"].copy() + gdf.at[0, "timeseries"]["hs"] = hs + gdf.at[0, "timeseries"]["tp"] = tp + gdf.at[0, "timeseries"]["wd"] = wd + gdf.at[0, "timeseries"]["ds"] = ds + else: + # First point, so need to generate df with constant water level + time = [self.model.config.get("tstart"), self.model.config.get("tstop")] + hs = [hs] * 2 + tp = [tp] * 2 + wd = [wd] * 2 + ds = [ds] * 2 + # Create DataFrame with columns time and wl + df = pd.DataFrame() + df["time"] = time + df["hs"] = hs + df["tp"] = tp + df["wd"] = wd + df["ds"] = ds + df = df.set_index("time") + gdf.at[0, "timeseries"] = df + else: + # Check if the timeseries is the same length as the first point + if len(gdf["timeseries"][0]) != len(self.data.iloc[0]["timeseries"]): + raise ValueError( + "Timeseries in gdf must be the same length as the first point in the boundary conditions data." + ) + + # Add to self.data + self.data = pd.concat([self.data, gdf], ignore_index=True) + + def delete(self, index: Union[int, List[int]]): + """Delete a single point from the SnapWave boundary conditions data. + + Parameters + ---------- + index : int or list of int + Index or list of indices of points to be deleted. + """ + + if self.data.empty: + return + + if not isinstance(index, list): + index = [index] + # Check if indices are within range + if any(x > (len(self.data.index) - 1) for x in index): + raise ValueError("One of the indices exceeds length of index range!") + self.data = self.data.drop(index).reset_index(drop=True) + + if self.data.empty: + self.model.config.set("snapwave_bndfile", None) + self.model.config.set("snapwave_btpfile", None) + self.model.config.set("snapwave_bwdfile", None) + self.model.config.set("snapwave_bdsfile", None) + # self.model.config.set("netbndsnapwavefile", None) + + def clear(self): + """Clean GeoDataFrame with boundary points.""" + self.data = gpd.GeoDataFrame() + + def set_timeseries( + self, + index: Union[int, List[int]] = None, + shape: str = "constant", + timestep: float = 600.0, + hs: float = 1.0, + tp: float = 10.0, + wd: float = 270.0, + ds: float = 20.0, + tpeak: float = 86400.0, + duration: float = 43200.0, + ): + """Applies time series boundary conditions for each point + Create numpy datetime64 array for time series with python datetime.datetime objects + + Parameters + ---------- + shape : str + Shape of the time series. Options are "constant" or "gaussian". + timestep : float + Time step [s] + hs : float + Wave height [m] + tp : float + Peak period [s] + wd : float + Wave direction [degrees] + ds : float + Directional spread [degrees] + tpeak : float + Time of the peak of the Gaussian wave [s] + duration : float + Duration of the Gaussian wave [s] + """ + + if self.data.empty: + return + + t0 = np.datetime64(self.model.config.get("tstart")) + t1 = np.datetime64(self.model.config.get("tstop")) + if shape == "constant": + dt = np.timedelta64(int((t1 - t0).astype(float) / 1e6), "s") + else: + dt = np.timedelta64(int(timestep), "s") + time = np.arange(t0, t1 + dt, dt) + dtsec = dt.astype(float) + # Convert time to seconds since tref + tsec = ( + (time - np.datetime64(self.model.config.get("tref"))) + .astype("timedelta64[s]") + .astype(float) + ) + nt = len(tsec) + if shape == "constant": + hs = [hs] * nt + tp = [tp] * nt + wd = [wd] * nt + ds = [ds] * nt + elif shape == "gaussian": + hs = hs * np.exp(-(((tsec - tpeak) / (0.25 * duration)) ** 2)) + tp = [tp] * nt + wd = [wd] * nt + ds = [ds] * nt + else: + # Not implemented + raise ValueError( + f"Shape {shape} not implemented for SnapWave boundary conditions!" + ) + + times = pd.date_range( + start=t0, end=t1, freq=pd.tseries.offsets.DateOffset(seconds=dtsec) + ) + + if index is None: + index = list(self.data.index) + elif not isinstance(index, list): + index = [index] + + for i in index: + df = pd.DataFrame() + df["time"] = times + df["hs"] = hs + df["tp"] = tp + df["wd"] = wd + df["ds"] = ds + df = df.set_index("time") + self.data.at[i, "timeseries"] = df + + def get_boundary_points_from_mask(self, min_dist=None, bnd_dist=5000.0): + # Should move this to mask? Yes. + if min_dist is None: + # Set minimum distance between to grid boundary points on polyline to 2 * dx + min_dist = self.model.quadtree_grid.data.attrs["dx"] * 2 + + mask = self.model.quadtree_grid.data["snapwave_mask"] + ibnd = np.where(mask == 2) + xz, yz = self.model.quadtree_grid.face_coordinates() + xp = xz[ibnd] + yp = yz[ibnd] + + # Make boolean array for points that are include in a polyline + used = np.full(xp.shape, False, dtype=bool) + + # Make list of polylines. Each polyline is a list of indices of boundary points. + polylines = [] + + while True: + if np.all(used): + # All boundary grid points have been used. We can stop now. + break + + # Find first the unused points + i1 = np.where(~used)[0][0] + + # Set this point to used + used[i1] = True + + # Start new polyline with index i1 + polyline = [i1] + + while True: + # Compute distances to all points that have not been used + xpunused = xp[~used] + ypunused = yp[~used] + # Get all indices of unused points + unused_indices = np.where(~used)[0] + + dst = np.sqrt((xpunused - xp[i1]) ** 2 + (ypunused - yp[i1]) ** 2) + if np.all(np.isnan(dst)): + break + inear = np.nanargmin(dst) + inearall = unused_indices[inear] + if dst[inear] < min_dist: + # Found next point along polyline + polyline.append(inearall) + used[inearall] = True + i1 = inearall + else: + # Last point found + break + + # Now work the other way + # Start with first point of polyline + i1 = polyline[0] + while True: + if np.all(used): + # All boundary grid points have been used. We can stop now. + break + # Now we go in the other direction + xpunused = xp[~used] + ypunused = yp[~used] + unused_indices = np.where(~used)[0] + dst = np.sqrt((xpunused - xp[i1]) ** 2 + (ypunused - yp[i1]) ** 2) + inear = np.nanargmin(dst) + inearall = unused_indices[inear] + if dst[inear] < min_dist: + # Found next point along polyline + polyline.insert(0, inearall) + used[inearall] = True + # Set index of next point + i1 = inearall + else: + # Last nearby point found + break + + if len(polyline) > 1: + polylines.append(polyline) + + gdf_list = [] + ip = 0 + # Transform to web mercator to get distance in metres + if self.model.crs.is_geographic: + transformer = Transformer.from_crs(self.model.crs, 3857, always_xy=True) + # Loop through polylines + for polyline in polylines: + x = xp[polyline] + y = yp[polyline] + points = [(x, y) for x, y in zip(x.ravel(), y.ravel())] + line = shapely.geometry.LineString(points) + if self.model.crs.is_geographic: + # Line in web mercator (to get length in metres) + xm, ym = transformer.transform(x, y) + pointsm = [(xm, ym) for xm, ym in zip(xm.ravel(), ym.ravel())] + linem = shapely.geometry.LineString(pointsm) + num_points = int(linem.length / bnd_dist) + 2 + else: + num_points = int(line.length / bnd_dist) + 2 + # Interpolate to new points + new_points = [ + line.interpolate(i / float(num_points - 1), normalized=True) + for i in range(num_points) + ] + # Loop through points in polyline + for point in new_points: + name = str(ip + 1).zfill(4) + d = { + "name": name, + "timeseries": pd.DataFrame(), + "geometry": point, + } + gdf_list.append(d) + ip += 1 + + self.data = gpd.GeoDataFrame(gdf_list, crs=self.model.crs) + + self.set_timeseries( + shape="constant", + timestep=600.0, + hs=1.0, + tp=10.0, + wd=270.0, + ds=20.0, + ) + + +# def to_fwf(df, fname, floatfmt=".3f"): +# indx = df.index.tolist() +# vals = df.values.tolist() +# for it, t in enumerate(vals): +# t.insert(0, indx[it]) +# content = tabulate(vals, [], tablefmt="plain", floatfmt=floatfmt) +# open(fname, "w").write(content) diff --git a/hydromt_sfincs/snapwave_conditions.py b/hydromt_sfincs/snapwave_conditions.py deleted file mode 100644 index 4a4b648d..00000000 --- a/hydromt_sfincs/snapwave_conditions.py +++ /dev/null @@ -1,72 +0,0 @@ -import geopandas as gpd -import numpy as np -import pandas as pd -import xarray as xr -from pathlib import Path -from typing import Union, List - -from hydromt.model.components import ModelComponent -from hydromt.model import Model -from hydromt_sfincs import utils - - -class SfincsSnapWaveConditions(ModelComponent): - def __init__( - self, - model: Model, - ): - self._filename: str = "snapwave.nc" - self._data: xr.Dataset = None - super().__init__( - model=model, - ) - - @property - def data(self) -> xr.Dataset: - """SnapWaveboundary conditions data. - - Return xr.Dataset - """ - if self._data is None: - self._initialize() - return self._data - - # Original HydroMT-SFINCS setup_ functions: - # setup_snapwave_forcing - - -# %% core HydroMT-SFINCS functions: -# _initialize -# read -# read_snapwave_netcdf - netsnapwavefile -# read_snapwave_points - snapwave_bndfile > FIXME don't support? -# read_snapwave_timeseries - bhsfile/btpfile/bwdfile/bdsfile> FIXME don't support? -# write -# write_snapwave_conditions - netsnapwavefile (or other default name) -# and/or: write_boundary_points - snapwave_bndfile + write_boundary_conditions_timeseries - bhsfile/btpfile/bwdfile/bdsfile> FIXME don't support? -# set -# set_points - gpd.GeoDataFrame with points -# set_timeseries - pd.DataFrame with timeseries (?) -# create: -# create_snapwave_timeseries -# add -# delete -# clear - -# %% DDB GUI focused additional functions: -# add_point (exception, call add for just one point) -# delete_point (exception, call delete for just one point) - -# to_gdf -# has_open_boundaries - -# get_boundary_points_from_mask -# set_timeseries -# set_timeseries_uniform -# set_conditions_at_point -# read_timeseries_file -# to_fwf -# inpolygon -# -# class SnapWaveBoundaryEnclosure > we are not going to support this > directly in snapwave_mask -# class SnapWaveMask > will be part of Mask class - exactly same functionality as for SFINCS mask diff --git a/hydromt_sfincs/utils.py b/hydromt_sfincs/utils.py index 5c43768a..8a69d98c 100644 --- a/hydromt_sfincs/utils.py +++ b/hydromt_sfincs/utils.py @@ -292,7 +292,7 @@ def write_timeseries( fn: Union[str, Path], df: Union[pd.DataFrame, pd.Series], tref: Union[str, datetime], - fmt: str = "%7.2f", + fmt: str = "%7.3f", ) -> None: """Write pandas.DataFrame to fixed width ascii timeseries files such as sfincs.bzs, sfincs.dis and sfincs.precip. The output time index is given in diff --git a/hydromt_sfincs/waterlevel_conditions.py b/hydromt_sfincs/waterlevel_conditions.py deleted file mode 100644 index 5f9a8925..00000000 --- a/hydromt_sfincs/waterlevel_conditions.py +++ /dev/null @@ -1,70 +0,0 @@ -import geopandas as gpd -import numpy as np -import pandas as pd -import xarray as xr -from pathlib import Path -from typing import Union, List - -from hydromt.model.components import ModelComponent -from hydromt.model import Model -from hydromt_sfincs import utils - - -class SfincsWaterlevelConditions(ModelComponent): - def __init__( - self, - model: Model, - ): - self._filename: str = "sfincs.bnd" # FIXME - List(str = "sfincs.bnd" and str = "sfincs.bzs" or str = "sfincs_netbndbzsbzi.nc") - self._data: gpd.GeoDataFrame = None - super().__init__( - model=model, - ) - - @property - def data(self) -> gpd.GeoDataFrame: - """Water level boundary conditions data. - - Return pd.GeoDataFrame - """ - if self._data is None: - self._initialize() - return self._data - - # Original HydroMT-SFINCS setup_ functions: - # setup_waterlevel_forcing - # setup_waterlevel_bnd_from_mask - # - # To add: setup_tidal_forcing - - -# %% core HydroMT-SFINCS functions: -# _initialize -# read -# read_waterlevel_points - bndfile -# read_waterlevel_timeseries - bzsfile -# read_waterlevel_astro - bcafile -# read_waterlevel_netcdf - netbndbzsbzifile --> ds = GeoDataset.from_netcdf(fn, crs=self.crs, chunks="auto") -# write -# write_boundary_points - bndfile + write_boundary_conditions_timeseries - bzsfile -# and/or: write_boundary_conditions - netbndbzsbzifile -# set -# set_points - gpd.GeoDataFrame with points -# set_timeseries - pd.DataFrame with timeseries (?) -# create: -# create_waterlevel_timeseries -# create_astro_timeseries -# add -# delete -# clear - -# %% DDB GUI focused additional functions: -# add_point (exception, call add for just one point) -# delete_point (exception, call delete for just one point) -# read_boundary_conditions_astro -# write_boundary_conditions_astro -# generate_bzs_from_bca -# get_boundary_points_from_mask -# set_timeseries -# read_timeseries_file -# to_fwf diff --git a/hydromt_sfincs/wave_makers.py b/hydromt_sfincs/wave_makers.py index 33935198..dfe183eb 100644 --- a/hydromt_sfincs/wave_makers.py +++ b/hydromt_sfincs/wave_makers.py @@ -1,34 +1,42 @@ +import logging import geopandas as gpd -import shapely + +# import shapely import pandas as pd from pathlib import Path -from typing import Union +from typing import TYPE_CHECKING, Union, List + +# import os +from os.path import abspath, join, exists from hydromt.model.components import ModelComponent from hydromt.model import Model from hydromt_sfincs import utils +if TYPE_CHECKING: + from hydromt_sfincs.sfincs import SfincsModel +logger = logging.getLogger(__name__) + class SfincsWaveMakers(ModelComponent): def __init__( self, - model: Model, + model: "SfincsModel", ): - self._filename: str = "sfincs.wvm" - self._data: gpd.GeoDataFrame = None + self.data = gpd.GeoDataFrame() super().__init__( model=model, ) - @property - def data(self) -> gpd.GeoDataFrame: - """Wavemakers lines data. + # @property + # def data(self) -> gpd.GeoDataFrame: + # """Cross-section lines data. - Return geopandas.GeoDataFrame - """ - if self._data is None: - self._initialize() - return self._data + # Return geopandas.GeoDataFrame + # """ + # if self._data is None: + # self._initialize() + # return self._data # %% core HydroMT-SFINCS functions: # _initialize @@ -36,168 +44,162 @@ def data(self) -> gpd.GeoDataFrame: # write # set # create - # add # delete # clear - def _initialize(self, skip_read=False) -> None: - """Initialize wavemakers lines.""" - if self._data is None: - self._data = gpd.GeoDataFrame() # FIXME - right? - if self.root.is_reading_mode() and not skip_read: - self.read() + # def _initialize(self, skip_read=False) -> None: + # """Initialize cross-section lines.""" + # if self._data is None: + # # self._data = dict() + # self._data = gpd.GeoDataFrame() # FIXME - right? + # if self.root.is_reading_mode() and not skip_read: + # self.read() - def read(self): - """Read in all wavemakers lines.""" - # Read input file: - struct = utils.read_geoms(self._filename) # =utils.py function - gdf = utils.linestring2gdf(struct, crs=self.model.crs) # =utils.py function + def read(self, filename: str | Path = None): + """Read SFINCS wave makers (*.wvm) file""" - self.set(gdf, merge=False) # Add to self._data + # Check that read mode is on + self.root._assert_read_mode() - def write(self, filename=None): # TODO - TL: filename=None - still needed? - """Write wvmfile.""" - # change precision of coordinates according to crs - if self.model.crs.is_geographic: - fmt = "%.6f" - else: - fmt = "%.1f" + # Get absolute file name and set it in config if crsfile is not None + abs_file_path = self.model.config.get_set_file_variable( + "wvmfile", value=filename + ) - # #TODO add to config - - # # If filename is not None: - # self.config.XXX - # self._filename = XXX - struct = utils.gdf2linestring(self.data) - utils.write_geoms( - self._filename, struct, stype="wvm", fmt=fmt - ) # =utils.py function + # Check if abs_file_path is None + if abs_file_path is None: + # File name not defined, so no wave makers in this model + return - # TODO - write also as geojson - TL: at what level do we want to do that? - # if self._write_gis: - # self.write_vector(variables=["geoms"]) + # Check if wvm file exists + if not abs_file_path.exists(): + raise FileNotFoundError(f"wave makers file not found: {abs_file_path}") - def set(self, gdf: gpd.GeoDataFrame, merge: bool = True): - """Set wavemakers lines. + # Read wvm file + struct = utils.read_geoms(abs_file_path) + gdf = utils.linestring2gdf(struct, crs=self.model.crs) - Arguments - --------- - gdf: geopandas.GeoDataFrame - Set GeoDataFrame with wavemakers lines to self.data - name: str - Geometry name. - """ - if not gdf.geometry.type.isin(["LineString"]).all(): - raise ValueError("wavemakers must be of type LineString.") - - # Clip points outside of model region: - within = gdf.within(self.model.region) # same as 'inpolygon' function - if within.all() == False: - raise ValueError("None of wavemakers fall within model domain.") - elif within.any() == False: - gdf = gdf[~within] - self.logger.info( - "Some of wavemakers fall out of model domain. Removing lines." - ) + # Add to self.data + self.set(gdf, merge=False) - if merge and self.data is not None: - gdf0 = self.data - gdf = gpd.GeoDataFrame(pd.concat([gdf, gdf0], ignore_index=True)) - self.logger.info("Adding new wavemakers to existing ones.") + def write(self, filename: str | Path = None): + """Write SFINCS wave makers (*.wvm) file, + and set wvmfile in config (if it was not already set)""" - # set gdf in self.data - self._data = gdf + # Check that data is not empty + if self.data.empty: + logger.info("No wave makers available to write.") + return - def create( - self, - locations: Union[str, Path, gpd.GeoDataFrame], - merge: bool = True, - **kwargs, - ): - """Create model wavemakers lines. - (old name: -) + # Set file name and get absolute path + abs_file_path = self.model.config.get_set_file_variable( + "wvmfile", + value=filename, + default="sfincs.wvm", + ) - Adds model layers: + # Change precision of coordinates according to crs + if self.model.crs.is_geographic: + fmt = "%11.6f" + else: + fmt = "%11.1f" - * **wvm** geom: wavemakers lines + # Get linestring geometries from gdf + struct = utils.gdf2linestring(self.data) + # Write to wvm file + utils.write_geoms(abs_file_path, struct, stype="wvm", fmt=fmt) + + # TODO - write also as geojson - TL: at what level do we want to do that? + # if self._write_gis: + # self.write_vector(variables=["geoms"]) + + def set(self, gdf: Union[gpd.GeoDataFrame, str, Path], merge: bool = True): + """Set SFINCS wave makers. Arguments --------- - locations: str, Path, gpd.GeoDataFrame, optional - Path, data source name, or geopandas object for wavemakers lines. - merge: bool, optional - If True, merge the new wavemakers lines with the existing ones. By default True. + str, Path, gpd.GeoDataFrame : + data source name, Path, or geopandas object with LineString geometries. + merge: bool + Merge with existing wave makers. If False, overwrite existing wave makers. """ - # FIXME ensure the catalog is loaded before adding any new entries - # self.data_catalog.sources TODO: check if still needed - gdf = self.data_catalog.get_geodataframe( - locations, geom=self.model.region, assert_gtype=None, **kwargs - ).to_crs(self.model.crs) + # Check if gdf is a string or Path. If so, read the file. + if isinstance(gdf, (str, Path)): + gdf = self.data_catalog.get_geodataframe( + gdf, geom=self.model.region, assert_gtype="LineString" + ).to_crs(self.model.crs) - # make sure MultiLineString are converted to LineString - gdf = gdf.explode(index_parts=True).reset_index(drop=True) - - self.set(gdf, merge) - - # TODO - add to config: self.model.config - # self.model.config(f"{name}file", f"sfincs.{name}") - # self.set_config(f"{name}file", f"sfincs.{name}") + if not gdf.geometry.type.isin(["LineString"]).all(): + raise ValueError("wave makers must be of type LineString.") + + # Check if any of the cross sections fall completely outside the model domain + # If so, give a warning and remove these lines + outside = gdf.disjoint(self.model.region) + if outside.any(): + logger.warning( + "Some wave makers fall outside model domain. Removing these lines." + ) + gdf = gdf[~outside] - def add( - self, - gdf: gpd.GeoDataFrame, - ): - """Add multiple lines to wavemakers. + # Check if there are any cross sections left + if gdf.empty: + logger.warning("All wave makers fall outside model domain!") + return - Arguments - --------- - gdf: gpd.GeoDataFrame - GeoDataFrame with locations and names of wavemakers lines to be added. - **NOTE** - coordinates of wavemakers in GeoDataFrame need to be in the same CRS as SFINCS model. - """ - self.set(gdf, merge=True) + if merge: + self.data = pd.concat([self.gdf, gdf], ignore_index=True) + logger.info("Adding new wave makers to existing ones") + else: + self.data = gdf + logger.info("Setting new wave makers") def delete( self, - index: int, # FIXME - should this be List(int)? + index: Union[list, int], ): - """Remove (multiple) line(s) from wavemakers. + """Remove one or more wave makers. Arguments --------- - index: int - Specify indices (int) of wavemakers(s) to be dropped from GeoDataFrame of wavemakers. + index: list, int + Specify wave makers to be dropped from GeoDataFrame. + If int, drop a single wave maker based on index. + If list, drop multiple wave makers based on index. """ - if index.any() > (len(self.data.index) - 1): # TODO - check if this is correct + # Turn int or str into list + if type(index) == int: + index = [index] + + # Check that any integer in list is not larger than the number of lines + if max(index) > (len(self.data) - 1) or min(index) < 0: raise ValueError("One of the indices exceeds length of index range!") - self.data.drop(index).reset_index(drop=True) - self.logger.info("Dropping line(s) from wavemakers") + # Drop lines from GeoDataFrame + self.data = self.data.drop(index).reset_index(drop=True) + logger.info("Dropping line(s) from wave makers") + + # Check if any cross sections are left + if self.data.empty: + logger.warning("All wave makers have been removed!") + # Set crsfile to None + self.model.config.set("wvmfile", None) def clear(self): - """Clean GeoDataFrame with wavemakerss.""" + """Clean GeoDataFrame with wave makers.""" self.data = gpd.GeoDataFrame() - - # %% DDB GUI focused additional functions: - # delete_polyline - # list_names - - def delete_polyline( - self, - index: int, - ): - """Remove single polyline from wavemakers. - This function sets the wanted index, after which the generic delete function is called. - - Arguments - --------- - name_or_index: str, int - Specify either name (str) or index (int) of cross-section to be dropped from GeoDataFrame. - """ - self.delete(index) # calls the generic delete function - return + # Set crsfile to None + self.model.config.set("wvmfile", None) def list_names(self): - """Give list of names of cross sections.""" - names = list(self.data.name) + """Give list of names of wave makers.""" + # The wave makers do not really have names, + # but we can use the index and turn into strings + names = [str(i + 1) for i in self.data.index] return names + + def snap_to_grid(self): + """Returns GeoDataFrame with wave makers snapped to model grid.""" + # TODO - this probably only works for quadtree grids for now + snap_gdf = self.model.grid.snap_to_grid(self.data) + return snap_gdf