From 344c213dcf5f1274753a238b0fadb26a6980a075 Mon Sep 17 00:00:00 2001 From: Javier Tausia Hoyal Date: Thu, 19 Jun 2025 08:58:37 +0200 Subject: [PATCH 1/8] [JTH] add kapi functions to clean binwaves --- bluemath_tk/topo_bathy/swan_grid.py | 183 ++++++++++++++++++---- bluemath_tk/waves/binwaves.py | 68 ++++++-- bluemath_tk/wrappers/swan/swan_wrapper.py | 178 ++++++++++++++++++--- 3 files changed, 365 insertions(+), 64 deletions(-) diff --git a/bluemath_tk/topo_bathy/swan_grid.py b/bluemath_tk/topo_bathy/swan_grid.py index 1aa1857..19b0644 100644 --- a/bluemath_tk/topo_bathy/swan_grid.py +++ b/bluemath_tk/topo_bathy/swan_grid.py @@ -2,7 +2,34 @@ import xarray as xr -def generate_grid_parameters(bathy_data: xr.DataArray) -> dict: +def generate_grid_parameters( + bathy_data: xr.DataArray, + buffer_distance: float = None, +) -> dict: + """ + Generate grid parameters for the SWAN model based on bathymetry. + + Parameters + ---------- + bathy_data : xr.DataArray + Bathymetry data with dimensions 'lat' and 'lon'. + + + Returns + ------- + dict + Dictionary with grid configuration for SWAN input. + + Raises + ------ + ValueError + If coord_type is not 'geographic' or 'cartesian'. + + Contact + ------- + @bellidog on GitHub + """ + """ Generate the grid parameters for the SWAN model. @@ -11,38 +38,136 @@ def generate_grid_parameters(bathy_data: xr.DataArray) -> dict: bathy_data : xr.DataArray Bathymetry data. Must have the following dimensions: - - lon: longitude - - lat: latitude + - lon/x: longitude or x coordinate + - lat/y: latitude or y coordinate Returns ------- dict Grid parameters for the SWAN model. - - Contact @bellidog on GitHub for more information. """ - 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"]) + # coord_type = 'geographic' + 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"] + ) + # coord_type = 'cartesian' + + # if coord_type not in ["geographic", "cartesian"]: + # raise ValueError("coord_type must be 'geographic' or 'cartesian'.") + + # use_int = coord_type == "cartesian" + # cast = int if use_int else float + + # Get the main parameters from user input + print("Please enter the following parameters:") + alpc = float(input("Enter rotation angle in degrees (alpc): ")) + + # 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: + xpc = float(input("Enter x origin (xpc): ")) + ypc = float(input("Enter y origin (ypc): ")) + xlenc = float(input("Enter grid length in x (xlenc): ")) + ylenc = float(input("Enter grid length in y (ylenc): ")) + + 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 + + # TODO: This is a temporary buffer distance, it should be adjusted to the actual grid size + # Get bounds with buffer + 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 + ) + + fixed_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 fixed_parameters, cropped + + else: + # Compute parameters from full bathymetry + return { + "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": alpc, # 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 + } 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/wrappers/swan/swan_wrapper.py b/bluemath_tk/wrappers/swan/swan_wrapper.py index fc1115a..9c454d4 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: + # Build the input files for the BinWaves model. + write_array_in_file(self.depth_array, f"{case_dir}/depth.dat") + # Write location file + 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") + ) From d638775507a75ecca9f1e67559cc25785375ecf1 Mon Sep 17 00:00:00 2001 From: Javier Tausia Hoyal Date: Thu, 19 Jun 2025 10:34:57 +0200 Subject: [PATCH 2/8] [JTH] little edit in wrapper --- bluemath_tk/wrappers/swan/swan_wrapper.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/bluemath_tk/wrappers/swan/swan_wrapper.py b/bluemath_tk/wrappers/swan/swan_wrapper.py index 9c454d4..a22c410 100644 --- a/bluemath_tk/wrappers/swan/swan_wrapper.py +++ b/bluemath_tk/wrappers/swan/swan_wrapper.py @@ -400,10 +400,10 @@ class BinWavesWrapper(SwanModelWrapper): """ def build_case(self, case_dir: str, case_context: dict) -> None: - # Build the input files for the BinWaves model. - write_array_in_file(self.depth_array, f"{case_dir}/depth.dat") - # Write location file - write_array_in_file(self.locations, f"{case_dir}/locations.loc") + 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( From f86e1bb46ab292fcab9ac54a4723dbaede918baa Mon Sep 17 00:00:00 2001 From: JenMontano Date: Thu, 31 Jul 2025 12:17:46 +0200 Subject: [PATCH 3/8] superPoint function including overlaping and fixing error around zero --- bluemath_tk/waves/superpoint.py | 43 ++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 9 deletions(-) diff --git a/bluemath_tk/waves/superpoint.py b/bluemath_tk/waves/superpoint.py index 14a0dcf..0c0d2f0 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 ------- @@ -53,14 +56,36 @@ def load_station_data(ds: xr.Dataset) -> xr.DataArray: superpoint_dataarray = xr.zeros_like( 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 From f4f7b545a4df3cc9f50bc5ed0f51885a6cc5bca5 Mon Sep 17 00:00:00 2001 From: JenMontano Date: Thu, 31 Jul 2025 17:33:48 +0200 Subject: [PATCH 4/8] [JMM] swan_grid changed and calibration CalVal fixed for NaNs --- bluemath_tk/topo_bathy/swan_grid.py | 66 +++++++++++++---------------- bluemath_tk/waves/calibration.py | 45 +++++++++++++------- 2 files changed, 59 insertions(+), 52 deletions(-) diff --git a/bluemath_tk/topo_bathy/swan_grid.py b/bluemath_tk/topo_bathy/swan_grid.py index 19b0644..5ac1bb3 100644 --- a/bluemath_tk/topo_bathy/swan_grid.py +++ b/bluemath_tk/topo_bathy/swan_grid.py @@ -1,10 +1,14 @@ import numpy as np import xarray as xr - def generate_grid_parameters( bathy_data: xr.DataArray, - buffer_distance: float = None, + alpc: float = 0, + xpc: float = None, + ypc: float = None, + xlenc: float = None, + ylenc: float = None, + buffer_distance: float = None, ) -> dict: """ Generate grid parameters for the SWAN model based on bathymetry. @@ -12,41 +16,35 @@ def generate_grid_parameters( Parameters ---------- bathy_data : xr.DataArray - Bathymetry data with dimensions 'lat' and 'lon'. - + Bathymetry data. + Must have the following dimensions: + - 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 Dictionary with grid configuration for SWAN input. - Raises - ------ - ValueError - If coord_type is not 'geographic' or 'cartesian'. - + Contact ------- - @bellidog on GitHub + @ """ - """ Generate the grid parameters for the SWAN model. - Parameters - ---------- - bathy_data : xr.DataArray - Bathymetry data. - Must have the following dimensions: - - lon/x: longitude or x coordinate - - lat/y: latitude or y coordinate - Returns ------- dict Grid parameters for the SWAN model. """ - # Determine coordinate system based on coordinate names coord_names = list(bathy_data.coords) @@ -54,6 +52,7 @@ def generate_grid_parameters( 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"]) + is_geographic = True # coord_type = 'geographic' else: x_coord = next( @@ -62,6 +61,7 @@ def generate_grid_parameters( y_coord = next( name for name in coord_names if name in ["y", "Y", "cy", "northing"] ) + is_geographic = False # coord_type = 'cartesian' # if coord_type not in ["geographic", "cartesian"]: @@ -70,9 +70,6 @@ def generate_grid_parameters( # use_int = coord_type == "cartesian" # cast = int if use_int else float - # Get the main parameters from user input - print("Please enter the following parameters:") - alpc = float(input("Enter rotation angle in degrees (alpc): ")) # Get resolution from cropped data grid_resolution_x = abs( @@ -83,13 +80,7 @@ def generate_grid_parameters( ) if alpc != 0: - xpc = float(input("Enter x origin (xpc): ")) - ypc = float(input("Enter y origin (ypc): ")) - xlenc = float(input("Enter grid length in x (xlenc): ")) - ylenc = float(input("Enter grid length in y (ylenc): ")) - angle_rad = np.radians(alpc) - # Create rotation matrix R = np.array( [ @@ -109,9 +100,8 @@ def generate_grid_parameters( # Translate to corner position x = rotated[:, 0] + xpc y = rotated[:, 1] + ypc - - # TODO: This is a temporary buffer distance, it should be adjusted to the actual grid size - # Get bounds with buffer + 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 @@ -127,7 +117,7 @@ def generate_grid_parameters( } # Note: slice from max to min for descending coordinates ) - fixed_parameters = { + grid_parameters = { "xpc": xpc, "ypc": ypc, "alpc": alpc, @@ -143,11 +133,11 @@ def generate_grid_parameters( "dxinp": grid_resolution_x, # resolution from cropped data "dyinp": grid_resolution_y, # resolution from cropped data } - return fixed_parameters, cropped + return grid_parameters, cropped, corners else: # Compute parameters from full bathymetry - return { + 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 @@ -161,7 +151,7 @@ def generate_grid_parameters( "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": alpc, # x-axis direction + "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( @@ -170,4 +160,6 @@ def generate_grid_parameters( "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/calibration.py b/bluemath_tk/waves/calibration.py index acd61b8..2878b3e 100644 --- a/bluemath_tk/waves/calibration.py +++ b/bluemath_tk/waves/calibration.py @@ -348,8 +348,13 @@ 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 @@ -557,7 +562,7 @@ def correct( correction_coeffs[n_part, :] = np.array( [ self.calibration_params["sea_correction"][ - int(peak_direction / self.direction_bin_size) + 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 @@ -568,7 +573,7 @@ def correct( correction_coeffs[n_part, :] = np.array( [ self.calibration_params["swell_correction"][ - int(peak_direction / self.direction_bin_size) + 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 @@ -594,7 +599,7 @@ def correct( * np.array( [ self.calibration_params["sea_correction"][ - int(peak_direction / self.direction_bin_size) + int(peak_direction / self.direction_bin_size) if peak_direction < 360 else 0 ] for peak_direction in corrected_data["Dirsea"] ] @@ -608,7 +613,7 @@ def correct( * np.array( [ self.calibration_params["swell_correction"][ - int(peak_direction / self.direction_bin_size) + 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 +758,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 +778,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="-") From 74d230b66d034c56a325a4a4ad01e0b2be59f968 Mon Sep 17 00:00:00 2001 From: JenMontano Date: Mon, 4 Aug 2025 09:23:44 +0200 Subject: [PATCH 5/8] bug NOAA downloader fixed --- bluemath_tk/downloaders/noaa/noaa_downloader.py | 2 +- ini_lat | 0 2 files changed, 1 insertion(+), 1 deletion(-) create mode 100644 ini_lat diff --git a/bluemath_tk/downloaders/noaa/noaa_downloader.py b/bluemath_tk/downloaders/noaa/noaa_downloader.py index 7883dfe..5483a0b 100644 --- a/bluemath_tk/downloaders/noaa/noaa_downloader.py +++ b/bluemath_tk/downloaders/noaa/noaa_downloader.py @@ -740,7 +740,7 @@ def _read_directional_file(self, file_path: Path) -> Optional[pd.DataFrame]: header_lines = [] while True: line = f.readline().strip() - if not line.startswith("#") or not line.startswith("YYYY"): + if not line.startswith("#") and not line.startswith("YYYY"): break header_lines.append(line) diff --git a/ini_lat b/ini_lat new file mode 100644 index 0000000..e69de29 From 37473661b074a5240c86108572d6601d6afe9268 Mon Sep 17 00:00:00 2001 From: JenMontano Date: Mon, 4 Aug 2025 13:05:47 +0200 Subject: [PATCH 6/8] small changes swan_grid --- bluemath_tk/topo_bathy/swan_grid.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/bluemath_tk/topo_bathy/swan_grid.py b/bluemath_tk/topo_bathy/swan_grid.py index 5ac1bb3..6df25f9 100644 --- a/bluemath_tk/topo_bathy/swan_grid.py +++ b/bluemath_tk/topo_bathy/swan_grid.py @@ -61,14 +61,6 @@ def generate_grid_parameters( y_coord = next( name for name in coord_names if name in ["y", "Y", "cy", "northing"] ) - is_geographic = False - # coord_type = 'cartesian' - - # if coord_type not in ["geographic", "cartesian"]: - # raise ValueError("coord_type must be 'geographic' or 'cartesian'.") - - # use_int = coord_type == "cartesian" - # cast = int if use_int else float # Get resolution from cropped data @@ -136,6 +128,8 @@ def generate_grid_parameters( return grid_parameters, cropped, corners else: + + # Compute parameters from full bathymetry grid_parameters = { "xpc": float(np.nanmin(bathy_data[x_coord])), # origin x From 9dcbe43b628194c17b2e72e48530c8e5c6482326 Mon Sep 17 00:00:00 2001 From: JenMontano Date: Tue, 5 Aug 2025 12:51:36 +0200 Subject: [PATCH 7/8] [JMM] remove not used file --- ini_lat | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 ini_lat diff --git a/ini_lat b/ini_lat deleted file mode 100644 index e69de29..0000000 From 4a5ced0141f563abc35c18edfe78167982466646 Mon Sep 17 00:00:00 2001 From: Javier Tausia Date: Tue, 5 Aug 2025 13:29:41 +0200 Subject: [PATCH 8/8] [JTH] little edits to format code --- bluemath_tk/topo_bathy/swan_grid.py | 27 ++++--------------- bluemath_tk/waves/calibration.py | 21 ++++++++++----- bluemath_tk/waves/superpoint.py | 40 +++++++++++++++++++---------- 3 files changed, 45 insertions(+), 43 deletions(-) diff --git a/bluemath_tk/topo_bathy/swan_grid.py b/bluemath_tk/topo_bathy/swan_grid.py index 6df25f9..219234e 100644 --- a/bluemath_tk/topo_bathy/swan_grid.py +++ b/bluemath_tk/topo_bathy/swan_grid.py @@ -1,6 +1,7 @@ import numpy as np import xarray as xr + def generate_grid_parameters( bathy_data: xr.DataArray, alpc: float = 0, @@ -8,7 +9,7 @@ def generate_grid_parameters( ypc: float = None, xlenc: float = None, ylenc: float = None, - buffer_distance: float = None, + buffer_distance: float = None, ) -> dict: """ Generate grid parameters for the SWAN model based on bathymetry. @@ -31,20 +32,8 @@ def generate_grid_parameters( ------- dict Dictionary with grid configuration for SWAN input. - - - Contact - ------- - @ - """ """ - Generate the grid parameters for the SWAN model. - Returns - ------- - dict - Grid parameters for the SWAN model. - """ # Determine coordinate system based on coordinate names coord_names = list(bathy_data.coords) @@ -52,8 +41,6 @@ def generate_grid_parameters( 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"]) - is_geographic = True - # coord_type = 'geographic' else: x_coord = next( name for name in coord_names if name in ["x", "X", "cx", "easting"] @@ -62,7 +49,6 @@ def generate_grid_parameters( 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 @@ -93,7 +79,7 @@ def generate_grid_parameters( 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 @@ -128,8 +114,6 @@ def generate_grid_parameters( return grid_parameters, cropped, corners else: - - # Compute parameters from full bathymetry grid_parameters = { "xpc": float(np.nanmin(bathy_data[x_coord])), # origin x @@ -145,7 +129,7 @@ def generate_grid_parameters( "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 + "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( @@ -154,6 +138,5 @@ def generate_grid_parameters( "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/calibration.py b/bluemath_tk/waves/calibration.py index 2878b3e..ec884e8 100644 --- a/bluemath_tk/waves/calibration.py +++ b/bluemath_tk/waves/calibration.py @@ -354,7 +354,6 @@ def _create_vec_direc(self, waves: np.ndarray, direcs: np.ndarray) -> np.ndarray else: bin_idx = int(direcs[i] / self.direction_bin_size) data[i, bin_idx] = waves[i] - return data @@ -562,7 +561,9 @@ def correct( correction_coeffs[n_part, :] = np.array( [ self.calibration_params["sea_correction"][ - int(peak_direction / self.direction_bin_size) if peak_direction < 360 else 0 #TODO: Check if this with Javi + 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 @@ -573,7 +574,9 @@ def correct( correction_coeffs[n_part, :] = np.array( [ self.calibration_params["swell_correction"][ - int(peak_direction / self.direction_bin_size) if peak_direction < 360 else 0 #TODO: Check if this with Javi + 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 @@ -599,7 +602,9 @@ def correct( * np.array( [ self.calibration_params["sea_correction"][ - int(peak_direction / self.direction_bin_size) if peak_direction < 360 else 0 + int(peak_direction / self.direction_bin_size) + if peak_direction < 360 + else 0 ] for peak_direction in corrected_data["Dirsea"] ] @@ -613,7 +618,9 @@ def correct( * np.array( [ self.calibration_params["swell_correction"][ - int(peak_direction / self.direction_bin_size) if peak_direction < 360 else 0 + int(peak_direction / self.direction_bin_size) + if peak_direction < 360 + else 0 ] for peak_direction in corrected_data[f"Dirswell{n_part}"] ] @@ -764,7 +771,7 @@ def plot_calibration_results(self) -> Tuple[Figure, list]: 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") @@ -784,7 +791,7 @@ def plot_calibration_results(self) -> Tuple[Figure, list]: 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") diff --git a/bluemath_tk/waves/superpoint.py b/bluemath_tk/waves/superpoint.py index 0c0d2f0..46cec05 100644 --- a/bluemath_tk/waves/superpoint.py +++ b/bluemath_tk/waves/superpoint.py @@ -56,36 +56,48 @@ def load_station_data(ds: xr.Dataset) -> xr.DataArray: superpoint_dataarray = xr.zeros_like( stations_data.isel({stations_dimension_name: 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) + 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) - + 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) - + 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) + superpoint_dataarray = xr.where( + overlap_mask, superpoint_dataarray / count_array, superpoint_dataarray + ) return superpoint_dataarray