From 957c07f2acfde801bc547a059ba0915b24d98332 Mon Sep 17 00:00:00 2001
From: veenstrajelmer <60435591+veenstrajelmer@users.noreply.github.com>
Date: Sat, 19 Oct 2024 18:23:29 +0200
Subject: [PATCH] exposed yearly tidal components (#159)

* exposed yearly tidal components

* updated whatsnew
---
 docs/whats-new.md                     |  1 +
 kenmerkendewaarden/gemiddeldgetij.py  | 17 ++++----------
 kenmerkendewaarden/tidalindicators.py | 33 +++++++++++++++++++--------
 tests/test_tidalindicators.py         | 15 ++++++++++++
 4 files changed, 44 insertions(+), 22 deletions(-)

diff --git a/docs/whats-new.md b/docs/whats-new.md
index 422ab8e..5633466 100644
--- a/docs/whats-new.md
+++ b/docs/whats-new.md
@@ -5,6 +5,7 @@
 ### Feat
 - expanded physical_break_dict in [#151](https://github.com/Deltares-research/kenmerkendewaarden/pull/151)
 - linear fit for slotgemiddelden (no nodal) in [#157](https://github.com/Deltares-research/kenmerkendewaarden/pull/157)
+- exposed yearly tidal components with `kw.calc_getijcomponenten()` in [#159](https://github.com/Deltares-research/kenmerkendewaarden/pull/159)
 
 
 ## 0.3.0 (2024-10-01)
diff --git a/kenmerkendewaarden/gemiddeldgetij.py b/kenmerkendewaarden/gemiddeldgetij.py
index 0ba623e..ff19364 100644
--- a/kenmerkendewaarden/gemiddeldgetij.py
+++ b/kenmerkendewaarden/gemiddeldgetij.py
@@ -8,7 +8,8 @@
 import hatyan
 import logging
 import matplotlib.pyplot as plt
-from kenmerkendewaarden.tidalindicators import calc_HWLWtidalrange
+from kenmerkendewaarden.tidalindicators import (calc_HWLWtidalrange,
+                                                calc_getijcomponenten)
 from kenmerkendewaarden.havengetallen import calc_havengetallen
 from kenmerkendewaarden.utils import TimeSeries_TimedeltaFormatter_improved
 from matplotlib.ticker import MaxNLocator, MultipleLocator
@@ -370,18 +371,8 @@ def get_gemgetij_components(data_pd_meas):
     # =============================================================================
     
     # components should not be reduced, since higher harmonics are necessary
-    const_list = hatyan.get_const_list_hatyan("year")  
-    hatyan_settings_ana = dict(
-        nodalfactors=True,
-        fu_alltimes=False,
-        xfac=True,
-        analysis_perperiod="Y",
-        return_allperiods=True,
-    )  # RWS-default settings
-    comp_frommeasurements_avg, comp_frommeasurements_allyears = hatyan.analysis(
-        data_pd_meas, const_list=const_list, **hatyan_settings_ana
-    )
-
+    comp_frommeasurements_avg, _ = calc_getijcomponenten(df_meas=data_pd_meas)
+    
     # #check if all years are available
     # comp_years = comp_frommeasurements_allyears['A'].columns
     # expected_years = tstop_dt.year-tstart_dt.year
diff --git a/kenmerkendewaarden/tidalindicators.py b/kenmerkendewaarden/tidalindicators.py
index 46b6d8e..45f5ab0 100644
--- a/kenmerkendewaarden/tidalindicators.py
+++ b/kenmerkendewaarden/tidalindicators.py
@@ -19,6 +19,7 @@
     "calc_HWLWtidalrange",
     "calc_hat_lat_fromcomponents",
     "calc_hat_lat_frommeasurements",
+    "calc_getijcomponenten",
 ]
 
 logger = logging.getLogger(__name__)
@@ -410,18 +411,32 @@ def calc_hat_lat_frommeasurements(df_meas_19y: pd.DataFrame) -> tuple:
         )
 
     # TODO: fu_alltimes=False makes the process significantly faster (default is True)
-    # TODO: xfac actually varies between stations (default is False), but different xfac has only very limited impact on the resulting hat/lat values
-    _, comp_all = hatyan.analysis(
-        df_meas_19y,
-        const_list="year",
-        analysis_perperiod="Y",
-        return_allperiods=True,
-        fu_alltimes=False,
-    )
+    # TODO: xfac actually varies between stations (default is False),
+    # but different xfac has only very limited impact on the resulting hat/lat values
+    _, comp_all = calc_getijcomponenten(df_meas=df_meas_19y, const_list=None)
 
-    # TODO: a frequency of 1min is better in theory, but 10min is faster and hat/lat values differ only 2mm for HOEKVHLD
+    # TODO: a frequency of 1min is better in theory, but 10min is faster and hat/lat 
+    # values differ only 2mm for HOEKVHLD
     df_pred = hatyan.prediction(comp_all, timestep="10min")
 
     lat = df_pred["values"].min()
     hat = df_pred["values"].max()
     return hat, lat
+
+
+def calc_getijcomponenten(df_meas, const_list=None):
+    if const_list is None:
+        const_list = hatyan.get_const_list_hatyan("year")  
+    # RWS-default settings
+    hatyan_settings_ana = dict(
+        nodalfactors=True,
+        fu_alltimes=False,
+        xfac=True,
+        analysis_perperiod="Y",
+        return_allperiods=True,
+    )
+    # analysis
+    comp_avg, comp_all = hatyan.analysis(
+        df_meas, const_list=const_list, **hatyan_settings_ana
+    )
+    return comp_avg, comp_all
diff --git a/tests/test_tidalindicators.py b/tests/test_tidalindicators.py
index a6e46dc..a5b320e 100644
--- a/tests/test_tidalindicators.py
+++ b/tests/test_tidalindicators.py
@@ -319,6 +319,21 @@ def test_calc_hat_lat_frommeasurements_tooshortperiod(df_meas_2010_2014):
     assert "please provide a timeseries of 19 years instead of 5 years" in str(e.value)
 
 
+@pytest.mark.unittest
+def test_calc_getijcomponenten(df_meas_2010_2014):
+    comp_av, comp_all = kw.calc_getijcomponenten(df_meas_2010_2014)
+    assert comp_av.shape == (95, 2)
+    assert comp_all.shape == (95, 10)
+    m2_av = comp_av.loc["M2"].values
+    m2_av_expected = np.array([ 0.7842112 , 86.07279444])
+    m2_all = comp_all.loc["M2"].values
+    m2_all_expected = np.array([0.79855603,  0.78044153,  0.78365574, 0.77764858,
+                                0.78078641, 85.91979015, 86.36060506, 85.71167985,
+                                86.24874909, 86.12879046])
+    assert np.allclose(m2_av, m2_av_expected)
+    assert np.allclose(m2_all, m2_all_expected)
+
+
 @pytest.mark.unittest
 def test_calc_HWLWtidalrange_aggers_input(df_ext_2010):
     with pytest.raises(ValueError) as e: