72 add hatlat computation from measurements (#74)

* added calc_hat_lat_frommeasurements including test

* extended plot_tidalindicators to also show constant values

* removed hat/lat example script and merged with kwk_process.py

* updated whatsnew

docs/whats-new.md

@@ -16,6 +16,7 @@
 - created `kw.plot_stations()` in [#57](https://github.com/Deltares-research/kenmerkendewaarden/pull/57)
 - clipping of timeseries on physical breaks with `kw.data_retrieve.clip_timeseries_physical_break()` (private) in [#61](https://github.com/Deltares-research/kenmerkendewaarden/pull/61) and [#64](https://github.com/Deltares-research/kenmerkendewaarden/pull/64)
 - added dedicated plotting functions in [#64](https://github.com/Deltares-research/kenmerkendewaarden/pull/64), [#66](https://github.com/Deltares-research/kenmerkendewaarden/pull/66) and [#68](https://github.com/Deltares-research/kenmerkendewaarden/pull/68)
+- added computation of hat/lat from measurements with `kw.calc_hat_lat_frommeasurements()` in [#74](https://github.com/Deltares-research/kenmerkendewaarden/pull/74)
 
 ### Fix
 - implemented workaround for pandas 2.2.0 with different rounding behaviour in [#69](https://github.com/Deltares-research/kenmerkendewaarden/pull/69)

examples/KWK_process.py

@@ -95,9 +95,17 @@
         # compute and plot tidal indicators
         dict_wltidalindicators = kw.calc_wltidalindicators(data_pd_meas_all)
         dict_HWLWtidalindicators = kw.calc_HWLWtidalindicators(data_pd_HWLW_all_12)
-        fig,ax = kw.plot_tidalindicators(dict_wltidalindicators, dict_HWLWtidalindicators)
-        fig.savefig(os.path.join(dir_indicators,f'tidal_indicators_{current_station}'))
 
+        # add hat/lat
+        df_meas_19y = data_pd_meas_all.loc["2001":"2019"]
+        hat, lat = kw.calc_hat_lat_frommeasurements(df_meas_19y)
+        dict_HWLWtidalindicators["hat"] = hat
+        dict_HWLWtidalindicators["lat"] = lat
+
+        # plot
+        fig, ax = kw.plot_tidalindicators(dict_wltidalindicators, dict_HWLWtidalindicators)
+        fig.savefig(os.path.join(dir_indicators,f'tidal_indicators_{current_station}'))
+
 
 
     #### SLOTGEMIDDELDEN

kenmerkendewaarden/tidalindicators.py

@@ -17,6 +17,7 @@
            "plot_tidalindicators",
            "calc_HWLWtidalrange",
            "calc_hat_lat_fromcomponents",
+           "calc_hat_lat_frommeasurements",
            ]
 
 logger = logging.getLogger(__name__)
@@ -148,6 +149,27 @@ def calc_wltidalindicators(data_wl_pd, min_count=None):
     return dict_wltidalindicators
 
 
+def plot_pd_series(indicators_dict, ax):
+    for key in indicators_dict.keys():
+        value = indicators_dict[key]
+        if not isinstance(value, pd.Series):
+            continue
+        if key.endswith("peryear"):
+            linestyle = "-"
+        elif key.endswith("permonth"):
+            linestyle = "--"
+        value.plot(ax=ax, label=key, linestyle=linestyle)
+        xmin = value.index.min()
+        xmax = value.index.max()
+
+    # separate loop for floats to make sure the xlim is already correct
+    for key in indicators_dict.keys():
+        value = indicators_dict[key]
+        if not isinstance(value, float):
+            continue
+        ax.hlines(value, xmin, xmax, linestyle="--", color="k", label=key, zorder=1)
+
+
 def plot_tidalindicators(indicators_wl:dict = None, indicators_ext = None):
     """
     plot tidalindicators
@@ -169,19 +191,7 @@ def plot_tidalindicators(indicators_wl:dict = None, indicators_ext = None):
     """
 
     fig, ax = plt.subplots(figsize=(12,6))
-
-    def plot_pd_series(indicators_dict, ax):
-
-        for key in indicators_dict.keys():
-            value = indicators_dict[key]
-            if not isinstance(value, pd.Series):
-                continue
-            if key.endswith("peryear"):
-                linestyle = "-"
-            elif key.endswith("permonth"):
-                linestyle = "--"
-            value.plot(ax=ax, label=key, linestyle=linestyle)
-
+
     if indicators_wl is not None:
         # TODO: maybe add an escape for if the station attr is not present
         station = indicators_wl['wl_mean_peryear'].attrs["station"]
@@ -216,9 +226,12 @@ def calc_HWLWtidalrange(ts_ext):
 def calc_hat_lat_fromcomponents(comp: pd.DataFrame) -> tuple:
     """
     Derive highest and lowest astronomical tide (HAT/LAT) from a component set.
-    The component set is used to make a tidal prediction for an arbitrary period of 19 years with a 1 minute interval. The max/min values of the predictions of all years are the HAT/LAT values.
-    The HAT/LAT is very dependent on the A0 of the component set. Therefore, the HAT/LAT values are relevant for the same year as the slotgemiddelde that is used to replace A0 in the component set. For instance, if the slotgemiddelde is valid for 2021.0, HAT and LAT are also relevant for that year.
-    The HAT/LAT values are also very dependent on the hatyan_settings used, in general it is important to use the same settings as used to derive the tidal components.
+    The component set is used to make a tidal prediction for an arbitrary period of 19 years with a 10 minute interval. 
+    The max/min values of the predictions of all years are the HAT/LAT values.
+    The HAT/LAT is very dependent on the A0 of the component set. Therefore, the HAT/LAT values are 
+    relevant for the same year as the slotgemiddelde that is used to replace A0 in the component set. 
+    For instance, if the slotgemiddelde is valid for 2021.0, HAT and LAT are also relevant for that year.
+    It is important to use the same tidal prediction settings as used to derive the tidal components.
     
     Parameters
     ----------
@@ -228,15 +241,15 @@ def calc_hat_lat_fromcomponents(comp: pd.DataFrame) -> tuple:
     Returns
     -------
     tuple
-        DESCRIPTION.
+        hat and lat values.
 
     """
 
     min_vallist_allyears = pd.Series(dtype=float)
     max_vallist_allyears = pd.Series(dtype=float)
     logger.info("generating prediction for 19 years")
     for year in range(2020,2039): # 19 arbitrary consequtive years to capture entire nodal cycle
-        times_pred_all = pd.date_range(start=dt.datetime(year,1,1), end=dt.datetime(year+1,1,1), freq='1min')
+        times_pred_all = pd.date_range(start=dt.datetime(year,1,1), end=dt.datetime(year+1,1,1), freq="10min")
         ts_prediction = hatyan.prediction(comp=comp, times=times_pred_all)
 
         min_vallist_allyears.loc[year] = ts_prediction['values'].min()
@@ -246,3 +259,41 @@ def calc_hat_lat_fromcomponents(comp: pd.DataFrame) -> tuple:
     hat = max_vallist_allyears.max()
     lat = min_vallist_allyears.min()
     return hat, lat
+
+
+def calc_hat_lat_frommeasurements(df_meas_19y: pd.DataFrame) -> tuple:
+    """
+    Derive highest and lowest astronomical tide (HAT/LAT) from a measurement timeseries of 19 years.
+    Tidal components are derived for each year of the measurement timeseries.
+    The resulting component sets are used to make a tidal prediction each year of the measurement timeseries with a 10 minute interval. 
+    The max/min values of the predictions of all years are the HAT/LAT values.
+    The HAT/LAT is very dependent on the A0 of the component sets. Therefore, the HAT/LAT values are 
+    relevant for the same period as the measurement timeseries. 
+
+    Parameters
+    ----------
+    df_meas_19y : pd.DataFrame
+        Measurements timeseries of 19 years.
+
+    Returns
+    -------
+    tuple
+        hat and lat values.
+
+    """
+
+    num_years = len(df_meas_19y.index.year.unique())
+    if num_years != 19:
+        raise ValueError(f"please provide a timeseries of 19 years instead of {num_years} years")
+
+    # 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: 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
+

tests/test_tidalindicators.py

@@ -114,8 +114,23 @@ def test_calc_hat_lat_fromcomponents(df_components_2010):
     df_components_2010_sel = df_components_2010.loc[["M2","M4","S2"]]
     # generate prediction and derive hat/lat
     hat, lat = kw.calc_hat_lat_fromcomponents(df_components_2010_sel)
-    assert np.isclose(hat, 1.2262061956302408)
-    assert np.isclose(lat, -0.8369650960044822)
+    assert np.isclose(hat, 1.2259179749801052)
+    assert np.isclose(lat, -0.8368954797393148)
+
+
+@pytest.mark.unittest
+def test_calc_hat_lat_frommeasurements(df_meas):
+    df_meas_19y = df_meas.loc["2001":"2019"]
+    hat, lat = kw.calc_hat_lat_frommeasurements(df_meas_19y)
+    assert np.isclose(hat, 1.6856114961274238)
+    assert np.isclose(lat, -1.0395726747948162)
+
+
+@pytest.mark.unittest
+def test_calc_hat_lat_frommeasurements_tooshortperiod(df_meas_2010_2014):
+    with pytest.raises(ValueError) as e:
+        hat, lat = kw.calc_hat_lat_frommeasurements(df_meas_2010_2014)
+    assert "please provide a timeseries of 19 years instead of 5 years" in str(e.value)
 
 
 @pytest.mark.unittest