diff --git a/bluemath_tk/topo_bathy/swan_grid.py b/bluemath_tk/topo_bathy/swan_grid.py index 1aa1857..219234e 100644 --- a/bluemath_tk/topo_bathy/swan_grid.py +++ b/bluemath_tk/topo_bathy/swan_grid.py @@ -2,47 +2,141 @@ import xarray as xr -def generate_grid_parameters(bathy_data: xr.DataArray) -> dict: +def generate_grid_parameters( + bathy_data: xr.DataArray, + alpc: float = 0, + xpc: float = None, + ypc: float = None, + xlenc: float = None, + ylenc: float = None, + buffer_distance: float = None, +) -> dict: """ - Generate the grid parameters for the SWAN model. + Generate grid parameters for the SWAN model based on bathymetry. Parameters ---------- bathy_data : xr.DataArray - Bathymetry data. + Bathymetry data. Must have the following dimensions: - - lon: longitude - - lat: latitude + - lon/x: longitude or x coordinate + - lat/y: latitude or y coordinate + alpc: float + Computational Grid Rotation angle in degrees. + xpc: float + X origin. + ypc: float + Y origin. Returns ------- dict - Grid parameters for the SWAN model. - - Contact @bellidog on GitHub for more information. + Dictionary with grid configuration for SWAN input. """ - return { - "xpc": int(np.nanmin(bathy_data.lon)), # x origin - "ypc": int(np.nanmin(bathy_data.lat)), # y origin - "alpc": 0, # x-axis direction - "xlenc": int( - np.nanmax(bathy_data.lon) - np.nanmin(bathy_data.lon) - ), # grid length in x - "ylenc": int( - np.nanmax(bathy_data.lat) - np.nanmin(bathy_data.lat) - ), # grid length in y - "mxc": len(bathy_data.lon) - 1, # number mesh x, una menos pq si no SWAN peta - "myc": len(bathy_data.lat) - 1, # number mesh y, una menos pq si no SWAN peta - "xpinp": np.nanmin(bathy_data.lon), # x origin - "ypinp": np.nanmin(bathy_data.lat), # y origin - "alpinp": 0, # x-axis direction - "mxinp": len(bathy_data.lon) - 1, # number mesh x - "myinp": len(bathy_data.lat) - 1, # number mesh y - "dxinp": abs( - bathy_data.lon[1].values - bathy_data.lon[0].values - ), # size mesh x (resolution in x) - "dyinp": abs( - bathy_data.lat[1].values - bathy_data.lat[0].values - ), # size mesh y (resolution in y) - } + # Determine coordinate system based on coordinate names + coord_names = list(bathy_data.coords) + + # Get coordinate variables + if any(name in ["lon", "longitude"] for name in coord_names): + x_coord = next(name for name in coord_names if name in ["lon", "longitude"]) + y_coord = next(name for name in coord_names if name in ["lat", "latitude"]) + else: + x_coord = next( + name for name in coord_names if name in ["x", "X", "cx", "easting"] + ) + y_coord = next( + name for name in coord_names if name in ["y", "Y", "cy", "northing"] + ) + + # Get resolution from cropped data + grid_resolution_x = abs( + bathy_data[x_coord][1].values - bathy_data[x_coord][0].values + ) + grid_resolution_y = abs( + bathy_data[y_coord][1].values - bathy_data[y_coord][0].values + ) + + if alpc != 0: + angle_rad = np.radians(alpc) + # Create rotation matrix + R = np.array( + [ + [np.cos(angle_rad), -np.sin(angle_rad)], + [np.sin(angle_rad), np.cos(angle_rad)], + ] + ) + + # Create unrotated rectangle corners + dx = np.array([0, xlenc, xlenc, 0, 0]) + dy = np.array([0, 0, ylenc, ylenc, 0]) + points = np.column_stack([dx, dy]) + + # Rotate points + rotated = np.dot(points, R.T) + + # Translate to corner position + x = rotated[:, 0] + xpc + y = rotated[:, 1] + ypc + corners = np.column_stack([x, y]) + + x_min = np.min(x) - buffer_distance + x_max = np.max(x) + buffer_distance + y_min = np.min(y) - buffer_distance + y_max = np.max(y) + buffer_distance + + print(f"Cropping bounds: x=[{x_min}, {x_max}], y=[{y_min}, {y_max}]") + + # Crop bathymetry + cropped = bathy_data.sel( + { + x_coord: slice(x_min, x_max), + y_coord: slice(y_max, y_min), + } # Note: slice from max to min for descending coordinates + ) + + grid_parameters = { + "xpc": xpc, + "ypc": ypc, + "alpc": alpc, + "xlenc": xlenc, + "ylenc": ylenc, + "mxc": int(np.round(xlenc / grid_resolution_x) - 1), + "myc": int(np.round(ylenc / grid_resolution_y) - 1), + "xpinp": np.nanmin(cropped[x_coord]), # x origin from cropped data + "ypinp": np.nanmin(cropped[y_coord]), # y origin from cropped data + "alpinp": 0, # x-axis direction + "mxinp": len(cropped[x_coord]) - 1, # number mesh x from cropped data + "myinp": len(cropped[y_coord]) - 1, # number mesh y from cropped data + "dxinp": grid_resolution_x, # resolution from cropped data + "dyinp": grid_resolution_y, # resolution from cropped data + } + return grid_parameters, cropped, corners + + else: + # Compute parameters from full bathymetry + grid_parameters = { + "xpc": float(np.nanmin(bathy_data[x_coord])), # origin x + "ypc": float(np.nanmin(bathy_data[y_coord])), # origin y + "alpc": alpc, # x-axis direction + "xlenc": float( + np.nanmax(bathy_data[x_coord]) - np.nanmin(bathy_data[x_coord]) + ), # grid length x + "ylenc": float( + np.nanmax(bathy_data[y_coord]) - np.nanmin(bathy_data[y_coord]) + ), # grid length y + "mxc": len(bathy_data[x_coord]) - 1, # num mesh x + "myc": len(bathy_data[y_coord]) - 1, # num mesh y + "xpinp": float(np.nanmin(bathy_data[x_coord])), # origin x + "ypinp": float(np.nanmin(bathy_data[y_coord])), # origin y + "alpinp": 0, # x-axis direction + "mxinp": len(bathy_data[x_coord]) - 1, # num mesh x + "myinp": len(bathy_data[y_coord]) - 1, # num mesh y + "dxinp": float( + abs(bathy_data[x_coord][1].values - bathy_data[x_coord][0].values) + ), # resolution x + "dyinp": float( + abs(bathy_data[y_coord][1].values - bathy_data[y_coord][0].values) + ), # resolution y + } + return grid_parameters diff --git a/bluemath_tk/waves/binwaves.py b/bluemath_tk/waves/binwaves.py index 163bb0f..743890f 100644 --- a/bluemath_tk/waves/binwaves.py +++ b/bluemath_tk/waves/binwaves.py @@ -16,18 +16,33 @@ def generate_swan_cases( - frequencies_array: np.ndarray, - directions_array: np.ndarray, + frequencies_array: np.ndarray = None, + directions_array: np.ndarray = None, + direction_range: tuple = (0, 360), + direction_divisions: int = 24, + direction_sector: tuple = None, + frequency_range: tuple = (0.035, 0.5), + frequency_divisions: int = 29, + gamma: float = 50, + spr: float = 2, ) -> xr.Dataset: """ Generate the SWAN cases monocromatic wave parameters. Parameters ---------- - directions_array : np.ndarray - The directions array. - frequencies_array : np.ndarray - The frequencies array. + frequencies_array : np.ndarray, optional + The frequencies array. If None, it is generated using frequency_range and frequency_divisions. + directions_array : np.ndarray, optional + The directions array. If None, it is generated using direction_range and direction_divisions. + direction_range : tuple + (min, max) range for directions in degrees. + direction_divisions : int + Number of directional divisions. + frequency_range : tuple + (min, max) range for frequencies in Hz. + frequency_divisions : int + Number of frequency divisions. Returns ------- @@ -35,9 +50,33 @@ def generate_swan_cases( The SWAN monocromatic cases Dataset with coordinates freq and dir. """ - # Wave parameters - gamma = 50 # waves gamma - spr = 2 # waves directional spread + # Auto-generate directions if not provided + if directions_array is None: + step = (direction_range[1] - direction_range[0]) / direction_divisions + directions_array = np.arange( + direction_range[0] + step / 2, direction_range[1], step + ) + + if direction_sector is not None: + start, end = direction_sector + if start < end: + directions_array = directions_array[ + (directions_array >= start) & (directions_array <= end) + ] + else: # caso circular, ej. 270–90 + directions_array = directions_array[ + (directions_array >= start) | (directions_array <= end) + ] + + # Auto-generate frequencies if not provided + if frequencies_array is None: + frequencies_array = np.geomspace( + frequency_range[0], frequency_range[1], frequency_divisions + ) + + # Constants for SWAN + gamma = gamma # waves gamma + spr = spr # waves directional spread # Initialize data arrays for each variable hs = np.zeros((len(directions_array), len(frequencies_array))) @@ -54,11 +93,11 @@ def generate_swan_cases( # Create xarray Dataset ds = xr.Dataset( - data_vars={ - "hs": (["dir", "freq"], hs), - "tp": (["dir", "freq"], tp), - "spr": (["dir", "freq"], spr_arr), - "gamma": (["dir", "freq"], gamma_arr), + { + "hs": (("dir", "freq"), hs), + "tp": (("dir", "freq"), tp), + "spr": (("dir", "freq"), spr_arr), + "gamma": (("dir", "freq"), gamma_arr), }, coords={ "dir": directions_array, @@ -115,6 +154,7 @@ def process_kp_coefficients( concatened_kp = output_kp_list[0] for file in output_kp_list[1:]: concatened_kp = xr.concat([concatened_kp, file], dim="case_num") + return concatened_kp.fillna(0.0).sortby("freq").sortby("dir") diff --git a/bluemath_tk/waves/calibration.py b/bluemath_tk/waves/calibration.py index acd61b8..ec884e8 100644 --- a/bluemath_tk/waves/calibration.py +++ b/bluemath_tk/waves/calibration.py @@ -348,7 +348,11 @@ def _create_vec_direc(self, waves: np.ndarray, direcs: np.ndarray) -> np.ndarray if direcs[i] < 0: direcs[i] = direcs[i] + 360 if direcs[i] > 0 and waves[i] > 0: - bin_idx = int(direcs[i] / self.direction_bin_size) + # Handle direction = 360° case by mapping to the first bin (0-22.5°) + if direcs[i] >= 360: + bin_idx = 0 + else: + bin_idx = int(direcs[i] / self.direction_bin_size) data[i, bin_idx] = waves[i] return data @@ -558,6 +562,8 @@ def correct( [ self.calibration_params["sea_correction"][ int(peak_direction / self.direction_bin_size) + if peak_direction < 360 + else 0 # TODO: Check if this with Javi ] for peak_direction in peak_directions.isel( part=n_part @@ -569,6 +575,8 @@ def correct( [ self.calibration_params["swell_correction"][ int(peak_direction / self.direction_bin_size) + if peak_direction < 360 + else 0 # TODO: Check if this with Javi ] for peak_direction in peak_directions.isel( part=n_part @@ -595,6 +603,8 @@ def correct( [ self.calibration_params["sea_correction"][ int(peak_direction / self.direction_bin_size) + if peak_direction < 360 + else 0 ] for peak_direction in corrected_data["Dirsea"] ] @@ -609,6 +619,8 @@ def correct( [ self.calibration_params["swell_correction"][ int(peak_direction / self.direction_bin_size) + if peak_direction < 360 + else 0 ] for peak_direction in corrected_data[f"Dirswell{n_part}"] ] @@ -753,11 +765,16 @@ def plot_calibration_results(self) -> Tuple[Figure, list]: ) # Plot sea wave climate - x, y, z = density_scatter( - self._data["Dirsea"].iloc[::10] * np.pi / 180, - self._data["Hsea"].iloc[::10], - ) - ax5.scatter(x, y, c=z, s=3, cmap="jet") + sea_dirs = self._data["Dirsea"].iloc[::10] * np.pi / 180 + sea_heights = self._data["Hsea"].iloc[::10] + # Filter out NaN and infinite values + valid_mask = np.isfinite(sea_dirs) & np.isfinite(sea_heights) + sea_dirs_valid = sea_dirs[valid_mask] + sea_heights_valid = sea_heights[valid_mask] + + if len(sea_dirs_valid) > 0: + x, y, z = density_scatter(sea_dirs_valid, sea_heights_valid) + ax5.scatter(x, y, c=z, s=3, cmap="jet") ax5.set_theta_zero_location("N", offset=0) ax5.set_xticklabels(["N", "NE", "E", "SE", "S", "SW", "W", "NW"]) ax5.xaxis.grid(True, color="lavender", linestyle="-") @@ -768,11 +785,16 @@ def plot_calibration_results(self) -> Tuple[Figure, list]: ax5.set_title("SEA $Wave$ $Climate$", pad=35, fontweight="bold") # Plot swell wave climate - x, y, z = density_scatter( - self._data["Dirswell1"].iloc[::10] * np.pi / 180, - self._data["Hswell1"].iloc[::10], - ) - ax6.scatter(x, y, c=z, s=3, cmap="jet") + swell_dirs = self._data["Dirswell1"].iloc[::10] * np.pi / 180 + swell_heights = self._data["Hswell1"].iloc[::10] + # Filter out NaN and infinite values + valid_mask = np.isfinite(swell_dirs) & np.isfinite(swell_heights) + swell_dirs_valid = swell_dirs[valid_mask] + swell_heights_valid = swell_heights[valid_mask] + + if len(swell_dirs_valid) > 0: + x, y, z = density_scatter(swell_dirs_valid, swell_heights_valid) + ax6.scatter(x, y, c=z, s=3, cmap="jet") ax6.set_theta_zero_location("N", offset=0) ax6.set_xticklabels(["N", "NE", "E", "SE", "S", "SW", "W", "NW"]) ax6.xaxis.grid(True, color="lavender", linestyle="-") diff --git a/bluemath_tk/waves/superpoint.py b/bluemath_tk/waves/superpoint.py index 14a0dcf..46cec05 100644 --- a/bluemath_tk/waves/superpoint.py +++ b/bluemath_tk/waves/superpoint.py @@ -7,6 +7,7 @@ def superpoint_calculation( stations_data: xr.DataArray, stations_dimension_name: str, sectors_for_each_station: Dict[str, Tuple[float, float]], + overlap_angle: float = 0.0, ) -> xr.DataArray: """ Join multiple station spectral data for each directional sector. @@ -20,6 +21,8 @@ def superpoint_calculation( sectors_for_each_station : Dict[str, Tuple[float, float]] Dictionary mapping each station ID to a tuple of (min_direction, max_direction) representing the directional sector for that station. + overlap_angle : float + Angle in degrees that defines the overlap between two sectors. Returns ------- @@ -54,13 +57,47 @@ def load_station_data(ds: xr.Dataset) -> xr.DataArray: stations_data.isel({stations_dimension_name: 0}) ) - for station_id, (dir_min, dir_max) in sectors_for_each_station.items(): - # Use where to select specific directions for each point - superpoint_dataarray += stations_data.sel( - {stations_dimension_name: station_id} - ).where( - (stations_data["dir"] >= dir_min) & (stations_data["dir"] < dir_max), - 0.0, + if overlap_angle == 0: + for station_id, (dir_min, dir_max) in sectors_for_each_station.items(): + if dir_min < dir_max: + mask = (stations_data["dir"] >= dir_min) & ( + stations_data["dir"] < dir_max + ) + else: + # Handle wrap-around (e.g., 350° to 10°) + mask = (stations_data["dir"] >= dir_min) | ( + stations_data["dir"] < dir_max + ) + superpoint_dataarray += stations_data.sel( + {stations_dimension_name: station_id} + ).where(mask, 0.0) + + else: + # With overlap - expand sectors and average overlaps + directions = stations_data["dir"] + count_array = xr.zeros_like(superpoint_dataarray) # Counter for overlaps + + for station_id, (dir_min, dir_max) in sectors_for_each_station.items(): + station_data = stations_data.sel({stations_dimension_name: station_id}) + + # Expand sector boundaries by overlap_angle + if (dir_max - dir_min) < 0: + mask = (directions >= dir_min - overlap_angle) | ( + directions <= dir_max + overlap_angle + ) + else: + mask = (directions >= dir_min - overlap_angle) & ( + directions <= dir_max + overlap_angle + ) + + # Add contribution where mask is true + superpoint_dataarray += station_data.where(mask, 0.0) + count_array += xr.where(mask, 1, 0) + + # Average where there are overlaps (count > 1) + overlap_mask = count_array > 1 + superpoint_dataarray = xr.where( + overlap_mask, superpoint_dataarray / count_array, superpoint_dataarray ) return superpoint_dataarray diff --git a/bluemath_tk/wrappers/swan/swan_wrapper.py b/bluemath_tk/wrappers/swan/swan_wrapper.py index fc1115a..a22c410 100644 --- a/bluemath_tk/wrappers/swan/swan_wrapper.py +++ b/bluemath_tk/wrappers/swan/swan_wrapper.py @@ -5,7 +5,9 @@ import numpy as np import pandas as pd import scipy.io as sio +import wavespectra import xarray as xr +from wavespectra.construct import construct_partition from .._base_wrappers import BaseModelWrapper from .._utils_wrappers import write_array_in_file @@ -47,6 +49,22 @@ class SwanModelWrapper(BaseModelWrapper): "value": None, "description": "Directional spread.", }, + "dir_dist": { + "type": str, + "choices": ["CIRCLE", "SECTOR"], + "value": "CIRCLE", + "description": "CIRCLE indicates that the spectral directions cover the full circle. SECTOR indicates that the spectral directions cover a limited sector of the circle.", + }, + "dir1": { + "type": float, + "value": None, + "description": "Only with SECTOR option. The direction of the right-hand boundary of the sector when looking outward from the sector (in degrees).", + }, + "dir2": { + "type": float, + "value": None, + "description": "Only with SECTOR option. The direction of the left-hand boundary of the sector when looking outward from the sector (in degrees).", + }, "mdc": { "type": int, "value": 24, @@ -62,6 +80,16 @@ class SwanModelWrapper(BaseModelWrapper): "value": 0.5, "description": "High value for frequency.", }, + "Freq_array": { + "type": np.ndarray, + "value": None, + "description": "Array of frequencies for the model.", + }, + "Dir_array": { + "type": np.ndarray, + "value": None, + "description": "Array of directions for the model.", + }, } available_launchers = { @@ -109,6 +137,7 @@ def __init__( output_dir: str, templates_name: dict = "all", depth_array: np.ndarray = None, + locations: np.ndarray = None, debug: bool = True, ) -> None: """ @@ -128,28 +157,14 @@ def __init__( ) if depth_array is not None: - self.depth_array = np.round(depth_array, 2) - - def build_case( - self, - case_context: dict, - case_dir: str, - ) -> None: - """ - Build the input files for a case. - - Parameters - ---------- - case_context : dict - The case context. - case_dir : str - The case directory. - """ + self.depth_array = np.round(depth_array, 5) + else: + self.depth_array = None - # Save depth array to file - if self.depth_array is not None: - depth_file = os.path.join(case_dir, "depth.dat") - write_array_in_file(array=self.depth_array, filename=depth_file) + if locations is not None: + self.locations = np.column_stack(locations) + else: + self.locations = None def list_available_output_variables(self) -> List[str]: """ @@ -308,3 +323,124 @@ def join_postprocessed_files( """ return xr.concat(postprocessed_files, dim="case_num") + + +def generate_fixed_parameters( + grid_parameters: dict, + freq_array: np.array, + dir_array: np.array, +) -> dict: + """ + Generate fixed parameters for the SWAN model based on grid parameters and frequency/direction arrays. + Parameters + ---------- + grid_parameters : dict + Dictionary with grid configuration for SWAN input. + freq_array : np.ndarray + Array of frequencies for the SWAN model. + dir_array : np.ndarray + Array of directions for the SWAN model. + Returns + ------- + dict + Dictionary with fixed parameters for the SWAN model. + """ + + dirs = np.sort(np.unique(dir_array)) % 360 + step = np.round(np.median(np.diff(np.sort(dirs))), 4) + + # Compute angular gaps between sorted directions (including wrap-around) + diffs = np.diff(np.concatenate([dirs, [dirs[0] + 360]])) + max_gap_idx = np.argmax(diffs) + + if np.isclose(diffs[max_gap_idx], step, atol=1e-2): + dir_dist = "CIRCLE" + dir1, dir2 = None, None + else: + dir_dist = "SECTOR" + dir1 = float((dirs[(max_gap_idx + 1) % len(dirs)]) % 360) # right-hand boundary + dir2 = float((dirs[max_gap_idx]) % 360) # left-hand boundary + print("Distribución direccional:", dir_dist) + if dir_dist == "SECTOR": + print(f"Direcciones de {dir1}° a {dir2}°") + + return { + "xpc": grid_parameters["xpc"], # origin x + "ypc": grid_parameters["ypc"], # origin y + "alpc": grid_parameters["alpc"], # x-axis direction + "xlenc": grid_parameters["xlenc"], # grid length x + "ylenc": grid_parameters["ylenc"], # grid length y + "mxc": grid_parameters["mxc"], # num mesh x + "myc": grid_parameters["myc"], # num mesh y + "xpinp": grid_parameters["xpinp"], # origin x for input grid + "ypinp": grid_parameters["ypinp"], # origin y for input grid + "alpinp": grid_parameters["alpinp"], # x-axis direction + "mxinp": grid_parameters["mxinp"], # num mesh x for input grid + "myinp": grid_parameters["myinp"], # num mesh y for input grid + "dxinp": grid_parameters["dxinp"], # resolution x for input grid + "dyinp": grid_parameters["dyinp"], # resolution y for input grid + "dir_dist": dir_dist, # direction distribution type + "dir1": dir1, # min direction + "dir2": dir2, # max direction + "freq_discretization": len(np.unique(freq_array)), # frequency discretization + "dir_discretization": int( + 360 / (np.unique(dir_array)[1] - np.unique(dir_array)[0]) + ), # direction discretization + "mdc": int( + 360 / (np.unique(dir_array)[1] - np.unique(dir_array)[0]) + ), # number of depth cases + "flow": float(np.min(np.unique(freq_array))), # low frequency limit + "fhigh": float(np.max(np.unique(freq_array))), # high frequency limit + } + + +class BinWavesWrapper(SwanModelWrapper): + """ + Wrapper example for the BinWaves model. + """ + + def build_case(self, case_dir: str, case_context: dict) -> None: + if self.depth_array is not None: + write_array_in_file(self.depth_array, f"{case_dir}/depth.dat") + if self.locations is not None: + write_array_in_file(self.locations, f"{case_dir}/locations.loc") + + # Construct the input spectrum + input_spectrum = construct_partition( + freq_name="jonswap", + freq_kwargs={ + "freq": np.geomspace( + case_context.get("flow", 0.035), + case_context.get("fhigh", 0.5), + case_context.get("freq_discretization", 29), + ), + # "freq": np.linspace(case_context.get("flow", 0.035), case_context.get("fhigh", 0.5), case_context.get("freq_discretization", 29)), + "fp": 1.0 / case_context.get("tp"), + "hs": case_context.get("hs"), + }, + dir_name="cartwright", + dir_kwargs={ + "dir": np.linspace(0, 360, case_context.get("dir_discretization", 24)), + "dm": case_context.get("dir"), + "dspr": case_context.get("spr"), + }, + ) + argmax_bin = np.argmax(input_spectrum.values) + mono_spec_array = np.zeros(input_spectrum.freq.size * input_spectrum.dir.size) + mono_spec_array[argmax_bin] = input_spectrum.sum(dim=["freq", "dir"]) + mono_spec_array = mono_spec_array.reshape( + input_spectrum.freq.size, input_spectrum.dir.size + ) + mono_input_spectrum = xr.Dataset( + { + "efth": (["freq", "dir"], mono_spec_array), + }, + coords={ + "freq": input_spectrum.freq, + "dir": input_spectrum.dir, + }, + ) + for side in ["N", "S", "E", "W"]: + wavespectra.SpecDataset(mono_input_spectrum).to_swan( + os.path.join(case_dir, f"input_spectra_{side}.bnd") + )