# SPDX-License-Identifier: GPL-2.0-only # Astronomical compatibility implementation after SIROCCO/T-UGOm 2013. # Upstream scientific work: Florent Lyard and Damien Allain, LEGOS/CNRS. # Original GPL-2.0 text: ../dist/data/TTB-COPYING.txt. """Reconstruct the 2013 PREVIMER atlas in its original TTB convention. No fitting, datum offset, inferred tide or network access is performed. The source NetCDF identifies tools 2.2.3, Hg 1145:2c4ace44f605, 2013-02-07 17:32. That exact historical revision is preserved in the producer's public Git repo: https://gitlab.in2p3.fr/sirocco/sirocco_tools/sirocco_tools/-/tree/119b33e963f7c508d87e4642c4d725201e573771 Wave definitions: src/tides.def; formulae/astronomy: src/tides.cpp. These numerical definitions intentionally reproduce the ORIGINAL analysis, including M1, KJ2, KQ1 and MK4 conventions later changed by the producer. The modern pyTMD/FES aliases are not equivalent. This is a compatibility model, not a new astronomical theory or a correction to the historical atlas. """ from __future__ import annotations import hashlib import json import sys from pathlib import Path ROOT = Path(__file__).resolve().parents[2] / 'essai-calcul-local' sys.path.insert(0, str(ROOT / '.deps')) import numpy as np import xarray as xr def configure_experiment_root(path): """Select local input NetCDF files; no network I/O or implicit downloads.""" global ROOT ROOT = Path(path).resolve() TTB_REVISION = '119b33e963f7c508d87e4642c4d725201e573771' # Each row: name, T,s,h,p,p1,shift_deg, xi,nu,nu_prime,nu_sec,R, f_formula. # The five fast arguments use the historic Darwin solar angle T=UT*15+180. _ROWS = ''' Mm 0 1 0 -1 0 0 0 0 0 0 0 73 Mf 0 2 0 0 0 0 -2 0 0 0 0 74 2Q1 1 -4 1 2 0 90 2 -1 0 0 0 75 Sig1 1 -4 3 0 0 90 2 -1 0 0 0 75 Q1 1 -3 1 1 0 90 2 -1 0 0 0 75 Ro1 1 -3 3 -1 0 90 2 -1 0 0 0 75 O1 1 -2 1 0 0 90 2 -1 0 0 0 75 MP1 1 -2 3 0 0 -90 0 -1 0 0 0 76 M1 1 -1 1 1 0 -90 1 -1 0 0 0 144 Ki1 1 -1 3 -1 0 -90 0 -1 0 0 0 76 Pi1 1 0 -2 0 1 90 0 0 0 0 0 0 P1 1 0 -1 0 0 90 0 0 0 0 0 0 K1 1 0 1 0 0 -90 0 0 -1 0 0 227 Psi1 1 0 2 0 -1 -90 0 0 0 0 0 0 Phi1 1 0 3 0 0 -90 0 0 0 0 0 0 Tta1 1 1 -1 1 0 -90 0 -1 0 0 0 76 J1 1 1 1 -1 0 -90 0 -1 0 0 0 76 OO1 1 2 1 0 0 -90 -2 -1 0 0 0 77 KQ1 1 3 1 -1 0 -90 -2 -1 0 0 0 15 E2 2 -5 4 1 0 0 2 -2 0 0 0 78 2N2 2 -4 2 2 0 0 2 -2 0 0 0 78 Mu2 2 -4 4 0 0 0 2 -2 0 0 0 78 N2 2 -3 2 1 0 0 2 -2 0 0 0 78 Nu2 2 -3 4 -1 0 0 2 -2 0 0 0 78 M2 2 -2 2 0 0 0 2 -2 0 0 0 78 La2 2 -1 0 1 0 180 2 -2 0 0 0 78 L2 2 -1 2 -1 0 180 2 -2 0 0 -1 215 T2 2 0 -1 0 1 0 0 0 0 0 0 0 S2 2 0 0 0 0 0 0 0 0 0 0 0 R2 2 0 1 0 -1 180 0 0 0 0 0 0 K2 2 0 2 0 0 0 0 0 0 -2 0 235 KJ2 2 1 2 -1 0 0 0 0 0 -2 0 79 MN4 4 -5 4 1 0 0 4 -4 0 0 0 1 M4 4 -4 4 0 0 0 4 -4 0 0 0 1 MS4 4 -2 2 0 0 0 2 -2 0 0 0 2 MK4 4 -2 4 0 0 0 2 -2 0 -2 0 4 M6 6 -6 6 0 0 0 6 -6 0 0 0 4 ''' WAVES = [] for line in _ROWS.strip().splitlines(): parts = line.split() x = [int(v) for v in parts[1:]] WAVES.append({'name': parts[0], 'darwin': x[:6], 'nodal': x[6:11], 'formula': x[11]}) NAMES = [wave['name'] for wave in WAVES] DARWIN = np.array([wave['darwin'][:5] for wave in WAVES]) SHIFT = np.radians([wave['darwin'][5] for wave in WAVES]) NODAL = np.array([wave['nodal'] for wave in WAVES]) RATES_DEG_PER_CENTURY = np.array([13149000.0, 481267.892, 36000.76892, 4069.0322056, 1.719175]) OMEGA = np.radians(DARWIN @ RATES_DEG_PER_CENTURY) / (36525 * 86400) PLACES = { 'billou': {'name': 'Plage du Bilou', 'targetLatitude': 48.704765, 'targetLongitude': -4.050046, 'targetSource': 'https://www.roscoff-tourisme.com/fr/fiche/equipements-de-loisirs/plage-du-billou-santec_TFOLOIBRE029V527UHR/', 'MANW': [117, 44]}, 'theven': {'name': 'Le Théven', 'targetLatitude': 48.710312, 'targetLongitude': -4.035724, 'targetSource': 'https://www.roscoff-tourisme.com/fr/fiche/equipements-de-loisirs/plage-de-theven-santec_TFOLOIBRE000FS0004D/', 'MANW': [121, 51]}, 'laber': {'name': 'Le Laber · Perharidy', 'targetLatitude': 48.72022, 'targetLongitude': -4.0066, 'targetSource': 'https://www.openstreetmap.org/way/163078411', 'MANW': [126, 57], 'targetNote': 'Western side of the Laber inlet, Perharidy; selected by the user. Nominal centre of the mapped beach, not a surveyed prediction point.'}, 'roscoffDiagnostic': {'name': 'Roscoff tide-gauge diagnostic', 'targetLatitude': 48.718426, 'targetLongitude': -3.965679, 'MANW': [121, 71]}, } def astronomical_basis(dates): """Return f and (V+u) radians under the atlas's historical 2013 convention.""" days = (np.asarray(dates, dtype='datetime64[ns]') - np.datetime64('1900-01-01', 'ns')) / np.timedelta64(1, 'D') centuries = np.asarray(days, dtype=float).reshape(-1) / 36525 angles_deg = np.array([180.0, 277.0256206, 280.1895015, 334.3837214, 281.2208568])[None, :] + centuries[:, None] * RATES_DEG_PER_CENTURY angles = np.radians(np.remainder(angles_deg, 360)) node = np.remainder(np.radians(259.1560563 - 1934.1423972 * centuries), 2 * np.pi) inclination = np.arccos(.913694997 - .035692561 * np.cos(node)) a = np.arctan(1.01883 * np.tan(node / 2)) b = np.arctan(.64412 * np.tan(node / 2)) xi = node - a - b xi = np.where(node > np.pi, xi - 2 * np.pi, xi) nu = a - b t2 = np.tan(inclination / 2) ** 2 P = angles[:, 3] - xi x1ra = np.sqrt(1 - 12 * t2 * np.cos(2 * P) + 36 * t2 ** 2) R = np.arctan(np.sin(2 * P) / (1 / (6 * t2) - np.cos(2 * P))) nu_prime = np.arctan(np.sin(2 * inclination) * np.sin(nu) / (np.sin(2 * inclination) * np.cos(nu) + .3347)) nu_sec = .5 * np.arctan(np.sin(inclination) ** 2 * np.sin(2 * nu) / (np.sin(inclination) ** 2 * np.cos(2 * nu) + .0727)) si, ch, sh = np.sin(inclination), np.cos(inclination / 2), np.sin(inclination / 2) factors = { 0: np.ones(len(centuries)), 73: (2 / 3 - si ** 2) / .5021, 74: si ** 2 / .1578, 75: si * ch ** 2 / .3800, 76: np.sin(2 * inclination) / .7214, 77: si * sh ** 2 / .0164, 78: ch ** 4 / .9154, 79: si ** 2 / .1565, 144: (1 - 10 * sh ** 2 + 15 * sh ** 4) * ch / .5873, 227: np.sqrt(.8965 * np.sin(2 * inclination) ** 2 + .6001 * np.sin(2 * inclination) * np.cos(nu) + .1006), 235: np.sqrt(19.0444 * si ** 4 + 2.7702 * si ** 2 * np.cos(2 * nu) + .0981), } factors.update({1: factors[78] ** 2, 2: factors[78], 4: factors[78] ** 3, 15: factors[235] * factors[75], 215: factors[78] * x1ra}) f = np.column_stack([factors[wave['formula']] for wave in WAVES]) nodal_angles = np.column_stack([xi, nu, nu_prime, nu_sec, R]) phase = angles @ DARWIN.T + SHIFT + nodal_angles @ NODAL.T return f, phase def haversine_m(latitude, longitude, target_lat, target_lon): lat, lon, tlat, tlon = np.radians([latitude, longitude, target_lat, target_lon]) a = np.sin((lat - tlat) / 2) ** 2 + np.cos(lat) * np.cos(tlat) * np.sin((lon - tlon) / 2) ** 2 return float(6371000 * 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))) def extract_points(domain='MANW'): """Extract exact unmodified amplitudes/phases at nearest fully wet cells.""" folder = ROOT / 'data' / 'ifremer' / domain manifest = json.loads((folder / 'download-manifest.json').read_text(encoding='utf-8')) paths = [folder / f'{name}-XE-{domain}-atlas.nc' for name in NAMES] with xr.open_dataset(paths[NAMES.index('M2')], engine='scipy') as ds: latitude = ds.latitude.values longitude = ds.longitude.values wet = np.isfinite(ds.XE_a.values) & np.isfinite(ds.XE_G.values) shape = list(ds.XE_a.shape) attrs = dict(ds.attrs) arrays = [] for path in paths: with xr.open_dataset(path, engine='scipy') as ds: a, g = ds.XE_a.values.copy(), ds.XE_G.values.copy() wet &= np.isfinite(a) & np.isfinite(g) arrays.append((a, g)) points = {} for key, place in PLACES.items(): target_lat, target_lon = place['targetLatitude'], place['targetLongitude'] la, lo = np.radians(latitude), np.radians(longitude) tla, tlo = np.radians([target_lat, target_lon]) h = np.sin((la - tla) / 2) ** 2 + np.cos(la) * np.cos(tla) * np.sin((lo - tlo) / 2) ** 2 j, i = np.unravel_index(np.argmin(np.where(wet, h, np.inf)), wet.shape) if domain in place and [int(j), int(i)] != place[domain]: raise ValueError(f'Nearest wet point changed for {key}: {(j, i)}') amp = np.array([a[j, i] for a, g in arrays], dtype=float) phase = np.array([g[j, i] for a, g in arrays], dtype=float) weights = amp * np.exp(-1j * np.radians(phase)) points[key] = {k: v for k, v in place.items() if k != 'MANW'} points[key].update(latitude=float(latitude[j, i]), longitude=float(longitude[j, i]), distanceM=haversine_m(latitude[j, i], longitude[j, i], target_lat, target_lon), cell={'j': int(j), 'i': int(i), 'indexBase': 0}, amplitudeM=amp.tolist(), phaseDegrees=phase.tolist(), weights=[[z.real, z.imag] for z in weights], meanOffsetM=0, datum='mean-relative') return {'domain': domain, 'shape': shape, 'netcdfAttributes': attrs, 'files': manifest['files'], 'ftpDirectory': manifest['directory'], 'waves': WAVES, 'models': points} def reconstruct(point, dates): f, phase = astronomical_basis(dates) weights = np.array([complex(*z) for z in point['weights']]) return np.real(f * np.exp(1j * phase) * weights[None, :]).sum(axis=1) def main(): import argparse parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--domain', choices=['MANW', 'FINIS'], default='MANW') parser.add_argument('--extract', type=Path) parser.add_argument('--experiment-root', type=Path, default=ROOT) args = parser.parse_args() configure_experiment_root(args.experiment_root) source = extract_points(args.domain) if args.extract: args.extract.write_text(json.dumps(source, indent=2, ensure_ascii=False) + '\n', encoding='utf-8') # Independent pre-existing event comparison against the 355 SHOM benchmarks. from predict_local import extrema from validate_local import compare, DATES, benchmarks events = extrema(DATES, reconstruct(source['models']['roscoffDiagnostic'], DATES)) result = {'domain': args.domain, 'convention': 'TTB-2.2.3-exact-2013', 'ttbRevision': TTB_REVISION, 'sourceConstantsFitted': False, 'diagnostic': source['models']['roscoffDiagnostic'], 'sources': {}} for reference in benchmarks['sources']: timing_reference = {**reference, 'events': [{**event, 'coef': None} for event in reference['events']]} _, summary = compare(events, timing_reference, []) result['sources'][reference['id']] = summary target = ROOT / 'results' / f'atlas-{args.domain.lower()}-validation.json' target.write_text(json.dumps(result, indent=2, ensure_ascii=False) + '\n', encoding='utf-8') print(json.dumps({'report': str(target), 'domain': args.domain, 'points': {key: {k: value for k, value in point.items() if k in ['cell', 'distanceM', 'latitude', 'longitude']} for key, point in source['models'].items()}, 'sources': result['sources']}, ensure_ascii=False)) if __name__ == '__main__': main()