67 add plot gemiddeldgetij (#68)

* added plot_gemiddeldgetij function

* improved timedelta formatting and added optional n-hourly ticks

* add test for plot_gemiddeldgetij

docs/whats-new.md

@@ -9,12 +9,13 @@
 - added neaptide tidal indicators for extremes in [#34](https://github.com/Deltares-research/kenmerkendewaarden/pull/34)
 - used threshold frequency instead of fixed index in `kw.overschrijding.blend_distributions` in [#38](https://github.com/Deltares-research/kenmerkendewaarden/pull/38)
 - dropped timezones consistently in `kw.calc_wltidalindicators()` and `kw.calc_HWLWtidalindicators()` to increase performance [#41](https://github.com/Deltares-research/kenmerkendewaarden/pull/41)
-- simplified methods for gemiddeld getij and reducing public functions to `kw.gemiddeld_getijkromme_av_sp_np()` in [#46](https://github.com/Deltares-research/kenmerkendewaarden/pull/46)
+- simplified methods for gemiddeld getij and reducing public functions to `kw.calc_gemiddeldgetij()` in [#46](https://github.com/Deltares-research/kenmerkendewaarden/pull/46)
 - simplified methods for havengetallen and reducing public functions to `kw.calc_havengetallen()` in [#48](https://github.com/Deltares-research/kenmerkendewaarden/pull/48)
 - simplified methods for slotgemiddelden and reducing public functions to `kw.calc_slotgemiddelden()` in [#62](https://github.com/Deltares-research/kenmerkendewaarden/pull/62)
 - increased test coverage in [#50](https://github.com/Deltares-research/kenmerkendewaarden/pull/50) and [#55](https://github.com/Deltares-research/kenmerkendewaarden/pull/55)
 - 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)
 
 
 ## 0.1.0 (2024-03-11)

examples/KWK_process.py

@@ -178,58 +178,33 @@
         pred_freq = "10s" # frequency decides accuracy of tU/tD and other timings (and is writing freq of BOI timeseries)
 
         # derive getijkrommes: raw, scaled to havengetallen, scaled to havengetallen and 12h25min period
-        prediction_av_raw, prediction_sp_raw, prediction_np_raw = kw.gemiddeld_getijkromme_av_sp_np(
-                                        df_meas=data_pd_meas_10y, df_ext=None,
-                                        freq=pred_freq, nb=0, nf=0, 
-                                        scale_extremes=False, scale_period=False)
-        prediction_av_corr, prediction_sp_corr, prediction_np_corr = kw.gemiddeld_getijkromme_av_sp_np(
-                                        df_meas=data_pd_meas_10y, df_ext=data_pd_HWLW_10y_12,
-                                        freq=pred_freq, nb=2, nf=2, 
-                                        scale_extremes=True, scale_period=False)
-        prediction_av_corr_boi, prediction_sp_corr_boi, prediction_np_corr_boi = kw.gemiddeld_getijkromme_av_sp_np(
-                                        df_meas=data_pd_meas_10y, df_ext=data_pd_HWLW_10y_12,
-                                        freq=pred_freq, nb=0, nf=10, 
-                                        scale_extremes=True, scale_period=True)
+        gemgetij_raw = kw.calc_gemiddeldgetij(df_meas=data_pd_meas_10y, df_ext=None,
+                                              freq=pred_freq, nb=0, nf=0, 
+                                              scale_extremes=False, scale_period=False)
+        gemgetij_corr = kw.calc_gemiddeldgetij(df_meas=data_pd_meas_10y, df_ext=data_pd_HWLW_10y_12,
+                                               freq=pred_freq, nb=1, nf=1, 
+                                               scale_extremes=True, scale_period=False)
+        gemgetij_corr_boi = kw.calc_gemiddeldgetij(df_meas=data_pd_meas_10y, df_ext=data_pd_HWLW_10y_12,
+                                                   freq=pred_freq, nb=0, nf=4, 
+                                                   scale_extremes=True, scale_period=True)
 
         # write boi timeseries to csv files # TODO: maybe convert timedeltaIndex to minutes instead?
-        prediction_av_corr_boi.to_csv(os.path.join(dir_gemgetij,f'gemGetijkromme_BOI_{current_station}_slotgem{year_slotgem}.csv'),float_format='%.3f',date_format='%Y-%m-%d %H:%M:%S')
-        prediction_sp_corr_boi.to_csv(os.path.join(dir_gemgetij,f'springtijkromme_BOI_{current_station}_slotgem{year_slotgem}.csv'),float_format='%.3f',date_format='%Y-%m-%d %H:%M:%S')
-        prediction_np_corr_boi.to_csv(os.path.join(dir_gemgetij,f'doodtijkromme_BOI_{current_station}_slotgem{year_slotgem}.csv'),float_format='%.3f',date_format='%Y-%m-%d %H:%M:%S')
+        for key in gemgetij_corr_boi.keys():
+            file_boi_csv = os.path.join(dir_gemgetij, f'Getijkromme_BOI_{key}_{current_station}_slotgem{year_slotgem}.csv')
+            gemgetij_corr_boi[key].to_csv(file_boi_csv, float_format='%.3f')
 
-
-        cmap = plt.get_cmap("tab10")
-
-        print(f'plot getijkromme trefHW: {current_station}')
-        fig_sum,ax_sum = plt.subplots(figsize=(14,7))
-        ax_sum.set_title(f'getijkromme trefHW {current_station}')
-        prediction_av_raw['values'].plot(ax=ax_sum, linestyle='--', color=cmap(0), linewidth=0.7, label='gem kromme, one')
-        prediction_av_corr['values'].plot(ax=ax_sum, color=cmap(0), label='gem kromme, corr')
-        prediction_sp_raw['values'].plot(ax=ax_sum, linestyle='--', color=cmap(1), linewidth=0.7, label='sp kromme, one')
-        prediction_sp_corr['values'].plot(ax=ax_sum, color=cmap(1), label='sp kromme, corr')
-        prediction_np_raw['values'].plot(ax=ax_sum, linestyle='--', color=cmap(2), linewidth=0.7, label='np kromme, one')
-        prediction_np_corr['values'].plot(ax=ax_sum, color=cmap(2), label='np kromme, corr')
-        ax_sum.set_xticks([x*3600e9 for x in range(-15, 25, 5)]) # nanoseconds units # TODO: make multiple of 12
-        ax_sum.legend(loc=4)
-        ax_sum.grid()
-        ax_sum.set_xlim([x*3600e9 for x in [-15.5,15.5]])
-        ax_sum.set_xlabel('hours since HW (ts are shifted to this reference)')
-        fig_sum.tight_layout()
+        fig_sum, ax_sum = kw.plot_gemiddeldgetij(gemgetij_dict=gemgetij_corr, gemgetij_dict_raw=gemgetij_raw, station=current_station, tick_hours=6)
         fig_sum.savefig(os.path.join(dir_gemgetij,f'gemgetij_trefHW_{current_station}'))
 
         print(f'plot BOI figure and compare to KW2020: {current_station}')
-        fig_boi,ax1_boi = plt.subplots(figsize=(14,7))
-        ax1_boi.set_title(f'getijkromme BOI {current_station}')
-        #plot gemtij/springtij/doodtij
-        prediction_av_corr_boi['values'].plot(ax=ax1_boi,color=cmap(0),label='prediction gemtij')
-        prediction_sp_corr_boi['values'].plot(ax=ax1_boi,color=cmap(1),label='prediction springtij')
-        prediction_np_corr_boi['values'].plot(ax=ax1_boi,color=cmap(2),label='prediction doodtij')
-        ax1_boi.set_xticks([x*3600e9 for x in range(0, 6*24, 12)]) # nanoseconds units
+        fig_boi, ax1_boi = kw.plot_gemiddeldgetij(gemgetij_dict=gemgetij_corr_boi, station=current_station, tick_hours=12)
 
-        #plot validation lines if available
+        # plot validation lines if available
         dir_vali_krommen = r'p:\archivedprojects\11205258-005-kpp2020_rmm-g5\C_Work\00_KenmerkendeWaarden\07_Figuren\figures_ppSCL_2\final20201211'
         file_vali_doodtijkromme = os.path.join(dir_vali_krommen,f'doodtijkromme_{current_station}_havengetallen{year_slotgem}.csv')
         file_vali_gemtijkromme = os.path.join(dir_vali_krommen,f'gemGetijkromme_{current_station}_havengetallen{year_slotgem}.csv')
         file_vali_springtijkromme = os.path.join(dir_vali_krommen,f'springtijkromme_{current_station}_havengetallen{year_slotgem}.csv')        
+        cmap = plt.get_cmap("tab10")
         if os.path.exists(file_vali_gemtijkromme):
             data_vali_gemtij = pd.read_csv(file_vali_gemtijkromme,index_col=0,parse_dates=True)
             ax1_boi.plot(data_vali_gemtij['Water Level [m]'],'--',color=cmap(0),linewidth=0.7,label='validation KW2020 gemtij')
@@ -240,11 +215,6 @@
             data_vali_doodtij = pd.read_csv(file_vali_doodtijkromme,index_col=0,parse_dates=True)
             ax1_boi.plot(data_vali_doodtij['Water Level [m]'],'--',color=cmap(2),linewidth=0.7, label='validation KW2020 doodtij')
 
-        ax1_boi.grid()
-        ax1_boi.legend(loc=4)
-        ax1_boi.set_xlabel('times since first av HW (start of ts)')
-        ax1_boi.set_xlim([x*3600e9 for x in [-2-4, 48-4]]) # TODO: make nicer xrange
-        fig_boi.tight_layout()
         fig_boi.savefig(os.path.join(dir_gemgetij,f'gemspringdoodtijkromme_BOI_{current_station}_slotgem{year_slotgem}.png'))
 
 

kenmerkendewaarden/gemiddeldgetij.py

@@ -7,18 +7,23 @@
 import pandas as pd
 import hatyan
 import logging
+import matplotlib.pyplot as plt
 from kenmerkendewaarden.tidalindicators import calc_HWLWtidalrange
 from kenmerkendewaarden.havengetallen import calc_havengetallen
+from kenmerkendewaarden.utils import TimeSeries_TimedeltaFormatter_improved
+from matplotlib.ticker import MaxNLocator, MultipleLocator
 
-__all__ = ["gemiddeld_getijkromme_av_sp_np",
+
+__all__ = ["calc_gemiddeldgetij",
+           "plot_gemiddeldgetij",
            ]
 
 logger = logging.getLogger(__name__)
 
 
-def gemiddeld_getijkromme_av_sp_np(df_meas: pd.DataFrame, df_ext: pd.DataFrame = None, 
-                             freq: str = "60sec", nb: int = 0, nf: int  = 0, 
-                             scale_extremes: bool = False, scale_period: bool = False, debug: bool = False):
+def calc_gemiddeldgetij(df_meas: pd.DataFrame, df_ext: pd.DataFrame = None, 
+                        freq: str = "60sec", nb: int = 0, nf: int  = 0, 
+                        scale_extremes: bool = False, scale_period: bool = False, debug: bool = False):
     """
     Generate an average tidal signal for average/spring/neap tide by doing a tidal 
     analysis on a timeseries of measurements. The (subsets/adjusted) resulting tidal components 
@@ -48,12 +53,8 @@ def gemiddeld_getijkromme_av_sp_np(df_meas: pd.DataFrame, df_ext: pd.DataFrame =
 
     Returns
     -------
-    prediction_av : pd.DataFrame
-        Dataframe with prediction for average tide.
-    prediction_sp : pd.DataFrame
-        Dataframe with prediction for spring tide.
-    prediction_np : pd.DataFrame
-        Dataframe with prediction for neap tide.
+    gemgetij_dict : dict
+        dictionary with Dataframes with gemiddeld getij for mean, spring and neap tide.
 
     """
     data_pd_meas_10y = df_meas
@@ -180,7 +181,77 @@ def gemiddeld_getijkromme_av_sp_np(df_meas: pd.DataFrame, df_ext: pd.DataFrame =
         prediction_np_corr_one = prediction_np_corr_one.resample(freq).nearest()
     prediction_np = repeat_signal(prediction_np_corr_one, nb=nb, nf=nf)
 
-    return prediction_av, prediction_sp, prediction_np
+    # combine in single dictionary
+    gemgetij_dict = {}
+    gemgetij_dict["mean"] = prediction_av
+    gemgetij_dict["spring"] = prediction_sp
+    gemgetij_dict["neap"] = prediction_np
+
+    return gemgetij_dict
+
+
+def plot_gemiddeldgetij(gemgetij_dict:dict, gemgetij_dict_raw:dict = None, station:str = None, tick_hours:int = None):
+    """
+    Default plotting function for gemiddeldgetij dictionaries.
+
+    Parameters
+    ----------
+    gemgetij_dict : dict
+        dictionary as returned from `kw.calc_gemiddeldgetij()`.
+    gemgetij_raw_dict : dict, optional
+        dictionary as returned from `kw.calc_gemiddeldgetij()` e.g. with uncorrected values. The default is None.
+    station : str, optional
+        station name, used in figure title. The default is None.
+    ticks_12h : bool, optional
+        whether to use xaxis ticks of 12 hours, otherwise automatic but less nice values
+
+    Returns
+    -------
+    fig : matplotlib.figure.Figure
+        Figure handle.
+    ax : matplotlib.axes._axes.Axes
+        Figure axis handle.
+
+    """
+    # TODO: prevent station argument
+    logger.info(f'plot getijkromme trefHW: {station}')
+    fig, ax = plt.subplots(figsize=(14,7))
+    cmap = plt.get_cmap("tab10")
+
+    if gemgetij_dict_raw is not None:
+        prediction_av_raw = gemgetij_dict_raw["mean"]
+        prediction_sp_raw = gemgetij_dict_raw["spring"]
+        prediction_np_raw = gemgetij_dict_raw["neap"]
+        prediction_av_raw['values'].plot(ax=ax, linestyle='--', color=cmap(0), linewidth=0.7, label='gemiddeldgetij mean (raw)')
+        prediction_sp_raw['values'].plot(ax=ax, linestyle='--', color=cmap(1), linewidth=0.7, label='gemiddeldgetij spring (raw)')
+        prediction_np_raw['values'].plot(ax=ax, linestyle='--', color=cmap(2), linewidth=0.7, label='gemiddeldgetij neap (raw)')
+
+    prediction_av_corr = gemgetij_dict["mean"]
+    prediction_sp_corr = gemgetij_dict["spring"]
+    prediction_np_corr = gemgetij_dict["neap"]
+    prediction_av_corr['values'].plot(ax=ax, color=cmap(0), label='gemiddeldgetij mean')
+    prediction_sp_corr['values'].plot(ax=ax, color=cmap(1), label='gemiddeldgetij spring')
+    prediction_np_corr['values'].plot(ax=ax, color=cmap(2), label='gemiddeldgetij neap')
+
+    ax.set_title(f'getijkrommes for {station}')
+    ax.legend(loc=4)
+    ax.grid()
+    ax.set_xlabel('time since high water')
+
+    # fix timedelta ticks
+    ax.xaxis.set_major_formatter(TimeSeries_TimedeltaFormatter_improved())
+    # put ticks at intervals of multiples of 3 and 6, resulting in whole seconds
+    ax.xaxis.set_major_locator(MaxNLocator(steps=[3,6], integer=True))
+    if tick_hours is not None:
+        # put ticks at fixed 12-hour intervals
+        ax.xaxis.set_major_locator(MultipleLocator(base=tick_hours*3600e9))
+    # the above avoids having to manually set tick locations based on hourly intervals (3600e9 nanoseconds)
+    # ax.set_xticks([x*3600e9 for x in range(-15, 25, 5)])
+    # ax.set_xlim([x*3600e9 for x in [-15.5,15.5]])
+
+    fig.tight_layout()
+
+    return fig, ax
 
 
 def get_gemgetij_components(data_pd_meas):

kenmerkendewaarden/utils.py

@@ -1,8 +1,49 @@
 # -*- coding: utf-8 -*-
 
+import numpy as np
+from matplotlib.ticker import Formatter
+
+
 def raise_extremes_with_aggers(df_ext):
     # TODO: alternatively we can convert 12345 to 12 here
     if len(df_ext["HWLWcode"].drop_duplicates()) != 2:
         raise ValueError("df_ext should only contain extremes (HWLWcode 1/2), "
                          "but it also contains aggers (HWLWcode 3/4/5). "
                          "You can convert with `hatyan.calc_HWLW12345to12()`")
+
+# TODO: create pandas issue suggesting this improvement
+class TimeSeries_TimedeltaFormatter_improved(Formatter):
+    """
+    Formats the ticks along an axis controlled by a :class:`TimedeltaIndex`.
+    based on pandas.plotting._matplotlib.converter.TimeSeries_TimedeltaFormatter
+    """
+
+    @staticmethod
+    def format_timedelta_ticks(x, pos, n_decimals: int) -> str:
+        """
+        Convert seconds to 'D days HH:MM:SS.F'
+        """
+        if x < 0:
+            negative = True
+            x = np.abs(x)
+        else:
+            negative = False
+        s, ns = divmod(x, 10**9)
+        m, s = divmod(s, 60)
+        h, m = divmod(m, 60)
+        d, h = divmod(h, 24)
+        decimals = int(ns * 10 ** (n_decimals - 9))
+        s = f"{int(h):02d}:{int(m):02d}:{int(s):02d}"
+        if n_decimals > 0:
+            s += f".{decimals:0{n_decimals}d}"
+        if d != 0:
+            s = f"{int(d):d} days {s}"
+        if negative:
+            s = '-'+s
+        return s
+
+    def __call__(self, x, pos: int = 0) -> str:
+        (vmin, vmax) = tuple(self.axis.get_view_interval())
+        n_decimals = min(int(np.ceil(np.log10(100 * 10**9 / abs(vmax - vmin)))), 9)
+        # print(x, pos, n_decimals)
+        return self.format_timedelta_ticks(x, pos, n_decimals)