diff --git a/bluemath_tk/config/paths.py b/bluemath_tk/config/paths.py index 35f38a4..57bb108 100644 --- a/bluemath_tk/config/paths.py +++ b/bluemath_tk/config/paths.py @@ -32,6 +32,10 @@ GEOOCEAN_THREDDS_DATA, "GEOOCEAN/SHyTCWaves_bulk/mda_mask_indices_lowres.nc", ), + "SHYTCWAVES_SPECTRA": op.join( + GEOOCEAN_CLUSTER_DATA, + "GEOOCEAN/SHyTCWaves/", + ), } @@ -52,16 +56,31 @@ def update_paths(new_paths: dict) -> None: PATHS.update(new_paths) -def get_paths() -> dict: +def get_paths(verbose: bool = True) -> dict: """ Get the paths dictionary. + Parameters + ---------- + verbose : bool, optional + Whether to print warnings about Thredds paths. Default is True. + Returns ------- dict Dictionary containing the paths. """ + if verbose: + for file_name, file_path in PATHS.items(): + if "thredds" in file_path: + print(f"WARNING: {file_name} is a Thredds path.") + print( + "You can update any path or add new paths with the update_paths function," + " from bluemath_tk.config.paths." + ) + print("Example: update_paths({'SHYTCWAVES_COEFS': '/new/path/to/data'})") + return PATHS diff --git a/bluemath_tk/core/geo.py b/bluemath_tk/core/geo.py index 656874a..c645ada 100644 --- a/bluemath_tk/core/geo.py +++ b/bluemath_tk/core/geo.py @@ -15,6 +15,9 @@ RAD2DEG = 180.0 / pi +# TODO: Check which functions are implemented in Pyproj library! + + def convert_to_radians(*args: Union[float, np.ndarray]) -> tuple: """ Convert degree inputs to radians. diff --git a/bluemath_tk/downloaders/noaa/noaa_downloader.py b/bluemath_tk/downloaders/noaa/noaa_downloader.py index 225df7f..7883dfe 100644 --- a/bluemath_tk/downloaders/noaa/noaa_downloader.py +++ b/bluemath_tk/downloaders/noaa/noaa_downloader.py @@ -151,10 +151,11 @@ def download_data( ---------- data_type : str The data type to download. - - 'bulk_parameters' - - 'wave_spectra' - - 'directional_spectra' - - 'wind_forecast' + - 'bulk_parameters' + - 'wave_spectra' + - 'directional_spectra' + - 'wind_forecast' + load_df : bool, optional Whether to load and return the DataFrame after downloading. Default is False. @@ -739,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("#"): + if not line.startswith("#") or not line.startswith("YYYY"): break header_lines.append(line) diff --git a/bluemath_tk/tcs/shytcwaves.py b/bluemath_tk/tcs/shytcwaves.py index 66a3521..420fb7c 100644 --- a/bluemath_tk/tcs/shytcwaves.py +++ b/bluemath_tk/tcs/shytcwaves.py @@ -7,11 +7,12 @@ import xarray as xr from tqdm import tqdm -from ..config.paths import PATHS from ..core.constants import EARTH_RADIUS from ..core.geo import geodesic_azimuth, geodesic_distance_azimuth, shoot from ..datamining.mda import find_nearest_indices +from ..waves.superpoint import superpoint_calculation from .tracks import ( + PATHS, get_vmean, historic_track_interpolation, historic_track_preprocessing, @@ -31,7 +32,7 @@ def storm2stopmotion(df_storm: pd.DataFrame) -> pd.DataFrame: """ - Generate stopmotion units from storm track segments. + Generate stopmotion segments from storm track. Parameters ---------- @@ -48,7 +49,7 @@ def storm2stopmotion(df_storm: pd.DataFrame) -> pd.DataFrame: Returns ------- pd.DataFrame - Stopmotion units with columns: + Stopmotion segments with columns: - vseg : Mean translational speed (kt) - dvseg : Speed variation (kt) - pseg : Segment pressure (mbar) @@ -71,7 +72,7 @@ def storm2stopmotion(df_storm: pd.DataFrame) -> pd.DataFrame: Notes ----- - The function generates stopmotion units methodology from 6-hour storm track segments: + The function generates stopmotion segments methodology from 6-hour storm track segments: A. Warmup segment (24h): - 4 segments to define start/end coordinates @@ -87,10 +88,7 @@ def storm2stopmotion(df_storm: pd.DataFrame) -> pd.DataFrame: is multiplied by -1). """ - # no need to remove NaTs,NaNs -> historic_track_preprocessing - # no need to remove wind/rmw -> historic_track_interpolation filled gaps - - # remove NaN + # remove NaNs df_ = df_storm.dropna() df_["time"] = df_.index.values @@ -243,12 +241,12 @@ def stopmotion_interpolation( t_prop: int = 42, ) -> Tuple[List[pd.DataFrame], List[pd.DataFrame]]: """ - Generate SWAN cases in cartesian coordinates from stopmotion parameterized units. + Generate SWAN cases in cartesian coordinates from stopmotion parameterized segments. Parameters ---------- df_seg : pd.DataFrame - Stopmotion parameterized units containing: + Stopmotion segments containing: - vseg : Mean translational speed (kt) - pseg : Segment pressure (mbar) - wseg : Maximum winds (kt) @@ -259,6 +257,7 @@ def stopmotion_interpolation( - dwseg : Wind variation (kt) - drseg : RMW variation (nmile) - daseg : Azimuth variation (degrees) + st : pd.DataFrame, optional Real storm track data. If None, MDA segments are used (unrelated to historic tracks). t_warm : int, optional @@ -272,6 +271,7 @@ def stopmotion_interpolation( ------- Tuple[List[pd.DataFrame], List[pd.DataFrame]] Two lists containing: + 1. List of storm track DataFrames with columns: - x, y : Cartesian coordinates (m) - lon, lat : Geographic coordinates (degrees) @@ -282,6 +282,7 @@ def stopmotion_interpolation( - vmax : Maximum winds (kt) - rmw : RMW (nmile) - latitude : Latitude with hemisphere sign + 2. List of empty wave event DataFrames with columns: - hs, t02, dir, spr : Wave parameters - U10, V10 : Wind components @@ -682,6 +683,7 @@ def stopmotion_st_bmu( - case : Storm segments - point : Control points - time : Time steps + Contains variables: - hs : Significant wave height - tp : Peak period @@ -836,6 +838,192 @@ def stopmotion_st_bmu( return xds_out +def stopmotion_st_spectra( + df_analogue: pd.DataFrame, + df_seg: pd.DataFrame, + st: pd.DataFrame, + cp_lon_ls: List[float], + cp_lat_ls: List[float], + cp_names: List[str] = [], + max_dist: float = 60, + list_out: bool = False, + tqdm_out: bool = False, + text_out: bool = True, + mode: str = "", +) -> Tuple[xr.Dataset, xr.Dataset]: + """ + Function to access the library analogue cases for a given storm track, + calculate distance and angle from the target segment origin to the control + point (relative coordinate system), and extract the directional wave + spectra at the closest node (for every analogue segment) + + Parameters + ---------- + df_analogue : pd.DataFrame + Analogue prerun segments from library. + df_seg : pd.DataFrame + Storm 6h-segments parameters. + st : pd.DataFrame + Storm track interpolated every 6h. + cp_lon_ls : List[float] + Control point geographical coordinates. + cp_lat_ls : List[float] + Control point geographical coordinates. + cp_names : List[str], optional + Control point names. Default is []. + max_dist : float, optional + Maximum distance [km] to extract closest node. Default is 60. + list_out : bool, optional + Whether to list output. Default is False. + tqdm_out : bool, optional + Whether to use tqdm. Default is False. + text_out : bool, optional + Whether to print text. Default is True. + mode : str, optional + Mode. Default is "". + + Returns + ------- + Tuple[xr.Dataset, xr.Dataset] + - xds_spec : Wave directional spectra (dim 'case') + - xds_bmu : BMU indices + + Notes + ----- + The function: + 1. Removes NaN values from df_seg and df_analogue + 2. Assigns time to df_seg + 3. Gets hemisphere sign + 4. Gets bmu (wavespectra reconstructed) + 5. Opens seg_sim dataset + """ + + # remove NaN + df_seg = df_seg[~df_seg.isna().any(axis=1)] + df_analogue = df_analogue[~df_analogue.isna().any(axis=1)] + + # assign time + df_seg["time"] = df_seg.index.values + + # get hemisphere + sign = np.sign(df_seg.lat.values[0]) + + # get bmu (wavespectra reconstructed) + # it provides 'time,bmu,ix_near,pos_nonan' (point,time) + xds_bmu = stopmotion_st_bmu( + df_analogue=df_analogue, + df_seg=df_seg, + st=st, + cp_lon_ls=cp_lon_ls, + cp_lat_ls=cp_lat_ls, + max_dist=max_dist, + list_out=list_out, + tqdm_out=tqdm_out, + text_out=text_out, + mode=mode, + ) + + # spectral energy + seg_sim = xr.open_dataset( + op.join(PATHS["SHYTCWAVES_SPECTRA"], "0000/spec_outpts_main.nc") + ) + efth_arr = np.full( + ( + seg_sim.frequency.size, + seg_sim.direction.size, + xds_bmu.point.size, + xds_bmu.time.size, + ), + np.nan, + ) # 38,72,cp,t + + if tqdm_out: + array = tqdm(range(xds_bmu.case.size)) + else: + array = range(xds_bmu.case.size) + for iseg in array: + # get storm segment 'i' + df_icase = df_seg.iloc[iseg] + df_ianalogue = df_analogue.iloc[iseg] + iseg_analogue = df_ianalogue.name # analogue id + + # --------------------------------------------------------------------- + # get storm coordinates at "seg_time" + seg_time = np.datetime64(df_icase.time) # timestamp to datetime64 + st_time = st.index.values[st.index.values == seg_time][0] + + # --------------------------------------------------------------------- + # get analogue segment from library + filename = "spec_outpts_main.nc" + p_analogue = op.join( + PATHS["SHYTCWAVES_SPECTRA"], f"{iseg_analogue:04d}", filename + ) + # load file + seg_sim = xr.open_dataset(p_analogue) # freq,dir,2088,48 + + # time array + hour_intervals = seg_sim.time.size + time = [st_time + np.timedelta64(1, "h") * i for i in range(hour_intervals)] + time_array = np.array(time) + + # get intersect time iseg vs xds_bmu + _, ix_time_st, ix_time_shy = np.intersect1d( + time_array, xds_bmu.time.values, return_indices=True + ) + + # find all closest grid points + shy_inear = xds_bmu.ix_near.values[iseg, :].astype("int64") # case,point + + # find bmu indices for iseg + in_pt, in_t = np.where(xds_bmu.bmu.values == iseg) + + # get indices of casei + in_t_ = in_t - ix_time_shy[0] + + # reorder spectral directions + base = 5 + if mode == "_lowres": + base = 10 # depends on the library dirs delta + efth_case = seg_sim.isel(point=shy_inear) # .isel(point=in_pt, time=in_t_) + if sign < 0: + efth_case["direction"] = 360 - seg_sim.direction.values + new_dirs = np.round( + efth_case.direction.values + base * round(df_icase.aseg / base) + 90, 1 + ) + else: + new_dirs = np.round( + efth_case.direction.values + base * round(df_icase.aseg / base) - 90, 1 + ) + new_dirs = np.mod(new_dirs, 360) + new_dirs[new_dirs > 270] -= 360 + efth_case["direction"] = new_dirs + efth_case = efth_case.sel(direction=seg_sim.direction.values) + + # insert spectral values for bmu=iseg + efth_arr[:, :, in_pt, in_t] = efth_case.efth.values[:, :, in_pt, in_t_] + + if text_out: + print("Inserting envelope spectra...", datetime.datetime.now()) + + # store dataset + xds_spec = xr.Dataset( + { + "efth": (("freq", "dir", "point", "time"), efth_arr), + "lon": (("point"), np.array(cp_lon_ls)), + "lat": (("point"), np.array(cp_lat_ls)), + "station": (("point"), np.array(cp_names)), + }, + coords={ + "freq": seg_sim.frequency.values, + "dir": seg_sim.direction.values, + "point": xds_bmu.point.values, + "time": xds_bmu.time.values, + }, + ) + + return xds_spec, xds_bmu + + def get_cp_radii_angle( st_lat: float, st_lon: float, @@ -961,133 +1149,6 @@ def get_mask_radii_angle(icase: int, mode: str = "") -> Tuple[np.ndarray, np.nda return rad, ang # [km], [º] -def get_mask_indices(mode: str = "") -> xr.Dataset: - """ - Create file with polar coordinates for SWAN case grids. - - Parameters - ---------- - mode : str, optional - Option to select SHyTCWaves library resolution ('', '_lowres'). Default is "". - - Returns - ------- - xr.Dataset - Dataset containing radii and angle arrays for each grid size: - - radii_sma, angle_sma : Small grid coordinates - - radii_med, angle_med : Medium grid coordinates - - radii_lar, angle_lar : Large grid coordinates - - Notes - ----- - The function: - 1. Loads library indices file with grid size information - 2. Gets one representative case for each grid size - 3. Calculates polar coordinates (radii/angle) for each grid - 4. For low resolution mode, uses half the number of radial points - """ - - # load library indices file - xds_mda = xr.open_dataset(PATHS["SHYTCWAVES_MDA_INDICES"]) - - # get one case for each grid size - case_sma = xds_mda.indices_small.values[0] - case_med = xds_mda.indices_medium.values[0] - case_lar = xds_mda.indices_large.values[0] - - rr, aa = [], [] - # get analogue case - for ic, case_id in enumerate([case_sma, case_med, case_lar]): - # get polar coords (radii/distance, angle) - rc, thetac = find_analogue_grid_coords(case_id) - - if mode == "_lowres": # new library of cases with half rings - rc = rc[:, ::2] - thetac = thetac[:, ::2] - - # to point dimension - rr.append(np.reshape(rc, -1)) # [m] - aa.append(np.reshape(thetac, -1)) # [º] - - # store all grid sizes mask indices - xd = xr.Dataset( - { - "radii_sma": (("point_sma",), rr[0]), - "angle_sma": (("point_sma",), aa[0]), - "radii_med": (("point_med",), rr[1]), - "angle_med": (("point_med",), aa[1]), - "radii_lar": (("point_lar",), rr[2]), - "angle_lar": (("point_lar",), aa[2]), - } - ) - - return xd - - -def find_analogue_grid_coords(icase: int) -> Tuple[np.ndarray, np.ndarray]: - """ - Find grid size and calculate polar coordinates for analogue segment case. - - Parameters - ---------- - icase : int - Analogue case ID. - - Returns - ------- - Tuple[np.ndarray, np.ndarray] - Two arrays containing: - - rc : Radial coordinates (m) - - thetac : Angular coordinates (degrees) - - Notes - ----- - The function: - 1. Loads MDA indices file to determine grid size - 2. Identifies if case uses small, medium, or large grid - 3. Sets appropriate number of rings based on grid size - 4. Creates domain centered at (x,y)=(0,0) - 5. Generates polar coordinates with custom spacing - """ - - # load MDA indices (grid sizes) - xds_ind_mda = xr.open_dataset(PATHS["SHYTCWAVES_MDA_INDICES"]) - - # get grid code - pos_small = np.where(icase == xds_ind_mda.indices_small)[0] - pos_medium = np.where(icase == xds_ind_mda.indices_medium)[0] - pos_large = np.where(icase == xds_ind_mda.indices_large)[0] - - # number of rings - if len(pos_small) == 1: - num, rings = 0, 29 - elif len(pos_medium) == 1: - num, rings = 1, 33 - elif len(pos_large) == 1: - num, rings = 2, 33 - - # get sizes - xsize_1 = xds_ind_mda.x_size_left.values[num] * 10**3 - xsize_2 = xds_ind_mda.x_size_right.values[num] * 10**3 - ysize_1 = xds_ind_mda.y_size_left.values[num] * 10**3 - ysize_2 = xds_ind_mda.y_size_right.values[num] * 10**3 - - # (cartesian convention) - # computational grid extent - res = 15 * 10**3 # km to m # TODO: resolution variable - - # create domain centered at (x,y)=(0,0) - lon = np.arange(-xsize_1, xsize_2, res) - lat = np.arange(-ysize_1, ysize_2, res) - - # custom output coordinates - _, _, _, _, rc, thetac = generate_polar_coords( - lon, lat, res, rings=rings, radii_ini=5000, angle_inc=5 - ) - - return rc, thetac - - def find_nearest( cp_rad_ls: List[float], cp_ang_ls: List[float], @@ -1198,19 +1259,22 @@ def historic2shytcwaves_cluster( Parameters ---------- path_save : str - Path to store results. + Base path to store results, without the file name. tc_name : str Storm name. storm : xr.Dataset Storm track dataset with standard IBTrACS variables: + Required: - - longitude, latitude : Storm coordinates - - pressure : Central pressure (mbar) - - maxwinds : Maximum sustained winds (kt) + - longitude, latitude : Storm coordinates + - pressure : Central pressure (mbar) + - maxwinds : Maximum sustained winds (kt) + Optional: - - rmw : Radius of maximum winds (nmile) - - dist2land : Distance to land - - basin : Basin identifier + - rmw : Radius of maximum winds (nmile) + - dist2land : Distance to land + - basin : Basin identifier + center : str IBTrACS center code. lon : np.ndarray @@ -1219,6 +1283,14 @@ def historic2shytcwaves_cluster( Latitude coordinates for swath calculation. dict_site : dict, optional Site data for superpoint building. Default is None. + Must contain: + - lonpts : Longitude coordinates + - latpts : Latitude coordinates + - namepts : Site names + - site : Site identifier + - sectors : Sectors + - deg_superposition : Superposition degree + calibration : bool, optional Whether to apply SHyTCWaves calibration. Default is True. mode : str, optional @@ -1243,25 +1315,25 @@ def historic2shytcwaves_cluster( 6. Saves results to specified directory """ - # A: stopmotion segmentation, 6h interval + # stopmotion segmentation, 6h interval df = historic_track_preprocessing( - storm, + xds=storm, center=center, - database_on=database_on, forecast_on=False, + database_on=database_on, st_param=st_param, ) - dt_int_seg = 6 * 60 # [minutes] constant segments + dt_int_minutes = 6 * 60 # [minutes] constant segments # optional: shytcwaves calibration of track parameters if calibration: - coef = get_coef_calibration() # lineal fitting + coef = get_coef_calibration() # linear fitting df["pressure"] = df["pressure"].values * (1 + coef[0]) + coef[1] df["maxwinds"] = np.nan st, _ = historic_track_interpolation( df, - dt_int_seg, + dt_int_minutes, interpolation=False, mode="mean", fit=True, @@ -1269,7 +1341,7 @@ def historic2shytcwaves_cluster( ) else: st, _ = historic_track_interpolation( - df, dt_int_seg, interpolation=False, mode="mean" + df, dt_int_minutes, interpolation=False, mode="mean" ) # skip when only NaN or 0 @@ -1283,18 +1355,17 @@ def historic2shytcwaves_cluster( st_trim = track_triming(st, lat[0], lon[0], lat[-1], lon[-1]) # store tracks for shytcwaves - st.to_pickle(op.join(path_save, "{0}_track.pkl".format(tc_name))) - st_trim.to_pickle(op.join(path_save, "{0}_track_trim.pkl".format(tc_name))) + st.to_pickle(op.join(path_save, f"{tc_name}_track.pkl")) + st_trim.to_pickle(op.join(path_save, f"{tc_name}_track_trim.pkl")) # parameterized segts (24h warmup + 6htarget) df_seg = storm2stopmotion(st_trim) if df_seg.shape[0] > 2: - print("st:", st.shape[0], "df_seg:", df_seg.shape[0]) + print(f"st: {st.shape[0]}, df_seg: {df_seg.shape[0]}") st_list, we_list = stopmotion_interpolation(df_seg, st=st_trim) - ####################################################################### - # B: analogue segments from library + # analogue segments from library df_mda = xr.open_dataset(PATHS["SHYTCWAVES_MDA"]).to_dataframe() ix_weights = [1] * 10 # equal weights ix = find_analogue(df_mda, df_seg, ix_weights) @@ -1302,19 +1373,11 @@ def historic2shytcwaves_cluster( df_analogue = df_mda.iloc[ix] df_analogue = analogue_endpoints(df_seg, df_analogue) - st_list_analogue, we_list_analogue = stopmotion_interpolation( - df_analogue, st=st_trim - ) - - ####################################################################### - # C: extract bulk envelope (to plot swaths) - + # extract bulk envelope (to plot swaths) if extract_bulk: mesh_lo, mesh_la = np.meshgrid(lon, lat) print( - "Number of segments: {0}, number of swath nodes: {1}".format( - len(st_list), mesh_lo.size - ) + f"Number of segments: {len(st_list)}, number of swath nodes: {mesh_lo.size}" ) if len(st_list) < max_segments: @@ -1329,40 +1392,40 @@ def historic2shytcwaves_cluster( ) # store xds_shy_bulk.to_netcdf( - op.join(path_save, "{0}_xds_shy_bulk.nc".format(tc_name)) + op.join(path_save, f"{tc_name}_xds_shy_bulk.nc") ) - ####################################################################### - # D: extract spectra envelope - if type(dict_site) == dict: + # extract spectra envelope + if isinstance(dict_site, dict): xds_shy_spec, _ = stopmotion_st_spectra( df_analogue, df_seg, st_trim, - dict_site["lonpts"], - dict_site["latpts"], + cp_lon_ls=dict_site["lonpts"], + cp_lat_ls=dict_site["latpts"], cp_names=dict_site["namepts"], - max_dist=60, - list_out=False, mode=mode, ) # store xds_shy_spec.to_netcdf( op.join( path_save, - "{0}_xds_shy_spec_{1}.nc".format(tc_name, dict_site["site"]), + f"{tc_name}_xds_shy_spec_{dict_site['site']}.nc", + ) + ) + + # build superpoint + xds_shy_sp = superpoint_calculation( + xds_shy_spec.efth, + "point", + dict_site["sectors"], + ) + # store + xds_shy_sp.to_netcdf( + op.join( + path_save, + f"{tc_name}_xds_shy_sp_{dict_site['site']}.nc", ) ) - ################################################################### - # # E: build superpoint - # stations = list(np.arange(0, xds_shy_spec.point.size)) - # xds_shy_sp = SuperPoint_Superposition(xds_shy_spec, stations, - # dict_site['sectors'], - # dict_site['deg_superposition']) - # # store - # xds_shy_sp.to_netcdf(op.join(p_save, - # '{0}_xds_shy_sp_{2}.nc'.format( - # tc_name, dict_site['site']))) - - print("Files stored.\n") + print("Files stored.") diff --git a/bluemath_tk/tcs/tracks.py b/bluemath_tk/tcs/tracks.py index 843bedc..df71065 100644 --- a/bluemath_tk/tcs/tracks.py +++ b/bluemath_tk/tcs/tracks.py @@ -7,11 +7,13 @@ import xarray as xr from numpy import polyfit -from ..config.paths import PATHS +from ..config.paths import get_paths from ..core.constants import EARTH_RADIUS from ..core.geo import geodesic_distance, geodesic_distance_azimuth, shoot -# Configuration dictionaries and constants +PATHS = get_paths(verbose=True) + +# Configuration dictionaries and constants for IBTrACS data centers_config_params: Dict[str, Dict[str, Union[str, List[str], float]]] = { "USA": { "id": "usa", @@ -99,7 +101,7 @@ def get_center_information(center: str = "WMO") -> Dict[str, Union[str, None]]: """ - Get the center information from the configuration parameters. + Get the center information from the configuration parameters for IBTrACS data. Parameters ---------- @@ -183,11 +185,6 @@ def check_and_plot_track_data(track_data: xr.Dataset) -> plt.Figure: - Longitude values are converted to [0-360º] convention - Wind speeds are converted to 1-minute average using center-specific factors - Variables that are entirely NaN are omitted from plots - - Examples - -------- - >>> fig = check_and_plot_track_data(track_data) - >>> plt.show() """ fig, axes = plt.subplots(1, 4, figsize=(20, 4)) @@ -259,6 +256,7 @@ def filter_track_by_basin(tracks_data: xr.Dataset, id_basin: str) -> xr.Dataset: >>> print(sp_tracks.storm.size) """ + # TODO: check whether just [0, :] is needed, as it is the only one used return tracks_data.isel( storm=np.where(tracks_data.basin.values[0, :].astype(str) == id_basin) ) @@ -300,7 +298,7 @@ def ibtracs_fit_pmin_wmax(ibtracs_data: xr.Dataset = None, N: int = 3) -> xr.Dat return xr.open_dataset(PATHS["SHYTCWAVES_COEFS"]) except Exception as e: - print(f"File could not be found: {PATHS['SHYTCWAVES_COEFS']}. Error: {e}") + print(f"File could not be opened: {PATHS['SHYTCWAVES_COEFS']}. Error: {e}") coef_fit = np.nan * np.zeros((len(all_centers), len(all_basins), N + 1)) pres_data = np.nan * np.zeros((len(all_centers), len(all_basins), 200000)) @@ -438,8 +436,7 @@ def wind2rmw(wmax: np.ndarray, vmean: np.ndarray, latitude: np.ndarray) -> np.nd Examples -------- >>> rmw = wind2rmw(np.array([100]), np.array([10]), np.array([25])) - >>> print(f"{rmw[0]:.1f}") - 30.5 + array([26.51]) """ pifac = np.arccos(-1) / 180 # pi/180 @@ -509,8 +506,7 @@ def get_vmean( Examples -------- >>> gamma, v, vx, vy = get_vmean(0, 0, 0, 1, 6) - >>> print(f"{v:.1f}") # Translation speed in km/h - 18.5 + (270.0, 18.55, 18.55, 5.68e-15) """ arcl_h, gamma_h = geodesic_distance_azimuth(lat2, lon2, lat1, lon1) # great circle @@ -518,10 +514,10 @@ def get_vmean( r = arcl_h * np.pi / 180.0 * EARTH_RADIUS # distance between coordinates [km] vmean = r / deltat # translation speed [km/h] - vu = vmean * np.sin((gamma_h + 180) * np.pi / 180) + vx = vmean * np.sin((gamma_h + 180) * np.pi / 180) vy = vmean * np.cos((gamma_h + 180) * np.pi / 180) - return gamma_h, vmean, vu, vy # [º], [km/h] + return gamma_h, vmean, vx, vy # [º], [km/h] def wind2pres( @@ -529,17 +525,18 @@ def wind2pres( ) -> np.ndarray: """ Convert maximum wind speeds to minimum central pressure using fitted coefficients. + As many other functions in this module, this works for IBTrACS data for the moment. Parameters ---------- xds_coef : xr.Dataset - Dataset containing polynomial fitting coefficients for Pmin-Wmax relationship + Dataset containing polynomial fitting coefficients for Pmin-Wmax relationship. st_wind : np.ndarray - Storm maximum winds in knots (1-min average) + Storm maximum winds in knots (1-min average). st_center : str - Storm center/agency identifier (e.g., 'WMO', 'TOKYO') + Storm center/agency identifier (e.g., 'WMO', 'TOKYO'). st_basin : str - Storm basin identifier (e.g., 'NA', 'WP') + Storm basin identifier (e.g., 'NA', 'WP'). Returns ------- @@ -556,11 +553,11 @@ def wind2pres( Examples -------- >>> pres = wind2pres(coef_dataset, np.array([100]), 'WMO', 'NA') - >>> print(f"{pres[0]:.1f}") - 955.2 + array([955.2]) """ pmin_pred = [] + for iwind in st_wind: # select Pmin-Wmax coefficients coeff = xds_coef.sel(center=st_center, basin=st_basin).coef.values[:] @@ -998,10 +995,16 @@ def historic_track_interpolation( # select Pmin-Wmax polynomial fitting coefficients (IBTrACS center,basin) xds_coef = ibtracs_fit_pmin_wmax() - coefs = xds_coef.sel( - center=st_center.encode("utf-8"), - basin=st_basin.astype("bytes"), - ).coef.values[:] + try: + coefs = xds_coef.sel( + center=st_center, + basin=st_basin, + ).coef.values[:] + except Exception: + coefs = xds_coef.sel( + center=st_center.encode("utf-8"), + basin=st_basin.astype("bytes"), + ).coef.values[:] p1, p2, p3, p4 = coefs[:, 0], coefs[:, 1], coefs[:, 2], coefs[:, 3] wind_estimate = ( diff --git a/bluemath_tk/waves/superpoint.py b/bluemath_tk/waves/superpoint.py index e69de29..14a0dcf 100644 --- a/bluemath_tk/waves/superpoint.py +++ b/bluemath_tk/waves/superpoint.py @@ -0,0 +1,66 @@ +from typing import Dict, Tuple + +import xarray as xr + + +def superpoint_calculation( + stations_data: xr.DataArray, + stations_dimension_name: str, + sectors_for_each_station: Dict[str, Tuple[float, float]], +) -> xr.DataArray: + """ + Join multiple station spectral data for each directional sector. + + Parameters + ---------- + stations_data : xr.DataArray + DataArray containing spectral data for multiple stations. + stations_dimension_name : str + Name of the dimension representing different stations in the DataArray. + 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. + + Returns + ------- + xr.DataArray + A new DataArray where each point is the sum of spectral data from all stations + for the specified directional sector. + + Notes + ----- + If your stations_data is saved in different files, you can load them all and then + concatenate them using xr.open_mfdataset function. Example below: + + ```python + files = [ + "path/to/station1.nc", + "path/to/station2.nc", + "path/to/station3.nc" + ] + + def load_station_data(ds: xr.Dataset) -> xr.DataArray: + return ds.efth.expand_dims("station") + + stations_data = xr.open_mfdataset( + files, + concat_dim="station", + preprocess=load_station_data, + ) + ``` + """ + + 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, + ) + + return superpoint_dataarray