diff --git a/docs/whats-new.md b/docs/whats-new.md index 1dec021..84ca3be 100644 --- a/docs/whats-new.md +++ b/docs/whats-new.md @@ -9,6 +9,7 @@ - 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 in [#46](https://github.com/Deltares-research/kenmerkendewaarden/pull/46) ## 0.1.0 (2024-03-11) diff --git a/examples/KWK_process.py b/examples/KWK_process.py index 85585f8..cb0c14e 100644 --- a/examples/KWK_process.py +++ b/examples/KWK_process.py @@ -8,6 +8,12 @@ import hatyan import kenmerkendewaarden as kw # pip install git+https://github.com/Deltares-research/kenmerkendewaarden +# set logging level to INFO to get log messages +# calling basicConfig is essential to set logging level of single module, format is optional +import logging +logging.basicConfig(format='%(message)s') +logging.getLogger("kenmerkendewaarden").setLevel(level="INFO") + # TODO: HW/LW numbers not always increasing (at havengetallen): ['HANSWT','BROUWHVSGT08','PETTZD','DORDT'] tstart_dt = pd.Timestamp(2011,1,1, tz="UTC+01:00") @@ -179,6 +185,7 @@ if compute_havengetallen and data_pd_HWLW_all is not None: print(f'havengetallen for {current_station}') + # TODO: havengetallen are different than p:\archivedprojects\11208031-010-kenmerkende-waarden-k\work\out_havengetallen_2021\havengetallen_2021_HOEKVHLD.csv # TODO: check culm_addtime and HWLWno+4 offsets. culm_addtime could also be 2 days or 2days +1h GMT-MET correction. 20 minutes seems odd since moonculm is about tidal wave from ocean # culm_addtime is a 2d and 2u20min correction, this shifts the x-axis of aardappelgrafiek @@ -214,135 +221,24 @@ if compute_gemgetij and data_pd_meas_all is not None: print(f'gem getijkrommen for {current_station}') - pred_freq_sec = 10 #TODO: frequency decides accuracy of tU/tD and other timings (and is writing freq of BOI timeseries) - - #TODO: add correctie havengetallen HW/LW av/sp/np met slotgemiddelde uit PLSS/modelfit (HW/LW av) + pred_freq = "10s" #TODO: frequency decides accuracy of tU/tD and other timings (and is writing freq of BOI timeseries) file_havget = os.path.join(dir_havget,f'havengetallen_{year_slotgem}_{current_station}.csv') - if not os.path.exists(file_havget): - raise Exception(f'havengetallen file does not exist: {file_havget}') - data_havget = pd.read_csv(file_havget, index_col=0) - for colname in ['HW_delay_median','LW_delay_median','getijperiod_median','duurdaling_median']: - data_havget[colname] = pd.to_timedelta(data_havget[colname]) - - HW_sp, LW_sp = data_havget.loc['spring',['HW_values_median','LW_values_median']] - HW_np, LW_np = data_havget.loc['neap',['HW_values_median','LW_values_median']] - HW_av, LW_av = data_havget.loc['mean',['HW_values_median','LW_values_median']] - - #derive components via TA on measured waterlevels - comp_frommeasurements_avg, comp_av = kw.get_gemgetij_components(data_pd_meas_10y) - - times_pred_1mnth = pd.date_range(start=dt.datetime(tstop_dt.year,1,1,0,0)-dt.timedelta(hours=12), end=dt.datetime(tstop_dt.year,2,1,0,0), freq=f'{pred_freq_sec} s') #start 12 hours in advance, to assure also corrected values on desired tstart - comp_av.attrs['nodalfactors'] = False #nodalfactors=False to guarantee repetative signal - comp_av.attrs['fu_alltimes'] = True # TODO: this is not true, but this setting is the default - prediction_av = hatyan.prediction(comp_av, times=times_pred_1mnth) - prediction_av_ext = hatyan.calc_HWLW(ts=prediction_av, calc_HWLW345=False) - - time_firstHW = prediction_av_ext.loc[prediction_av_ext['HWLWcode']==1].index[0] #time of first HW - ia1 = prediction_av_ext.loc[time_firstHW:].index[0] #time of first HW - ia2 = prediction_av_ext.loc[time_firstHW:].index[2] #time of second HW - prediction_av_one = prediction_av.loc[ia1:ia2] - prediction_av_ext_one = prediction_av_ext.loc[ia1:ia2] - - # ============================================================================= - # Hatyan predictie voor 1 jaar met gemiddelde helling maansbaan (voor afleiden spring-doodtijcyclus) >> predictie zonder nodalfactors instead - # ============================================================================= - """ - uit: gemiddelde getijkrommen 1991.0 - Voor de ruwe krommen voor springtij en doodtij is het getij voorspeld voor een jaar met gemiddelde helling maansbaan - met uitsluitend zuivere combinaties van de componenten M2 en S2: - tabel: Gebruikte componenten voor de spring- en doodtijkromme - SM, 3MS2, mu2, M2, S2, 2SM2, 3MS4, M4, MS4, - 4MS6, M6, 2MS6, M8, 3MS8, M10, 4MS10, M12, 5MS12 - - In het aldus gemodelleerde getij is de vorm van iedere getijslag, gegeven de getijfase, identiek. - Vervolgens is aan de hand van de havengetallen een springtij- en een doodtijkromme geselecteerd. - - #NOTE: background on choice of components - #below is different than provided list, these shallow ones are extra: ['S4','2SM6','M7','4MS4','2(MS)8','3M2S10','4M2S12'] - #shallow relations, derive 'zuivere harmonischen van M2 en S2' (this means averaging over eenmaaldaagse componenten, but why is that chosen?) - #adding above extra components or oneday freqs, gives a modulation and therefore there is no repetative signal anymore. Apperantly all components in this list have an integer number of periods in one springneap cycle? - dummy,shallowrel,dummy = hatyan.foreman.get_foreman_shallowrelations() - bool_M2S2only = shallowrel[1].isin([1,2]) & shallowrel[3].isin(['M2','S2']) & shallowrel[5].isin(['M2','S2',np.nan]) & shallowrel.index.isin(const_list_year) - shallowdeps_M2S2 = shallowrel.loc[bool_M2S2only,:5] - print(shallowdeps_M2S2) - """ - components_sn = ['A0','SM','3MS2','MU2','M2','S2','2SM2','3MS4','M4','MS4','4MS6','M6','2MS6','M8','3MS8','M10','4MS10','M12','5MS12'] - - #make prediction with springneap components with nodalfactors=False (alternative for choosing a year with a neutral nodal factor). Using 1yr instead of 1month does not make a difference in min/max tidal range and shape, also because of nodalfactors=False. (when using more components, there is a slight difference) - comp_frommeasurements_avg_sncomp = comp_frommeasurements_avg.loc[components_sn] - comp_frommeasurements_avg_sncomp.attrs['nodalfactors'] = False #nodalfactors=False to make independent on chosen year - comp_frommeasurements_avg_sncomp.attrs['fu_alltimes'] = True # TODO: this is not true, but this setting is the default - prediction_sn = hatyan.prediction(comp_frommeasurements_avg_sncomp, times=times_pred_1mnth) - - prediction_sn_ext = hatyan.calc_HWLW(ts=prediction_sn, calc_HWLW345=False) - - #selecteer getijslag met minimale tidalrange en maximale tidalrange (werd geselecteerd adhv havengetallen in 1991.0 doc) - prediction_sn_ext = kw.calc_HWLWtidalrange(prediction_sn_ext) - - time_TRmax = prediction_sn_ext.loc[prediction_sn_ext['HWLWcode']==1,'tidalrange'].idxmax() - is1 = prediction_sn_ext.loc[time_TRmax:].index[0] - is2 = prediction_sn_ext.loc[time_TRmax:].index[2] - - time_TRmin = prediction_sn_ext.loc[prediction_sn_ext['HWLWcode']==1,'tidalrange'].idxmin() - in1 = prediction_sn_ext.loc[time_TRmin:].index[0] - in2 = prediction_sn_ext.loc[time_TRmin:].index[2] - - #select one tideperiod for springtide and one for neaptide - prediction_sp_one = prediction_sn.loc[is1:is2] - prediction_sp_ext_one = prediction_sn_ext.loc[is1:is2] - prediction_np_one = prediction_sn.loc[in1:in2] - prediction_np_ext_one = prediction_sn_ext.loc[in1:in2] - - # plot selection of neap/spring - fig, (ax1,ax2) = hatyan.plot_timeseries(ts=prediction_sn,ts_ext=prediction_sn_ext) - ax1.plot(prediction_sp_one['values'],'r') - ax1.plot(prediction_np_one['values'],'r') - ax1.legend(labels=ax1.get_legend_handles_labels()[1]+['kromme spring','kromme neap'],loc=4) - ax1.set_ylabel('waterstand [m]') - ax1.set_title(f'spring- en doodtijkromme {current_station}') - fig.savefig(os.path.join(dir_gemgetij,f'springdoodtijkromme_{current_station}_slotgem{year_slotgem}.png')) - - - #timeseries for gele boekje (av/sp/np have different lengths, time is relative to HW of av and HW of sp/np are shifted there) #TODO: is this product still necessary? - print(f'reshape_signal GEMGETIJ: {current_station}') - prediction_av_one_trefHW = kw.ts_to_trefHW(prediction_av_one) # repeating one is not necessary for av, but easier to do the same for av/sp/np - prediction_av_corr_one = kw.reshape_signal(prediction_av_one, prediction_av_ext_one, HW_goal=HW_av, LW_goal=LW_av, tP_goal=None) - prediction_av_corr_rep5 = kw.repeat_signal(prediction_av_corr_one, nb=2, na=2) - prediction_av_corr_rep5_trefHW = kw.ts_to_trefHW(prediction_av_corr_rep5,HWreftime=ia1) - - print(f'reshape_signal SPRINGTIJ: {current_station}') - prediction_sp_one_trefHW = kw.ts_to_trefHW(prediction_sp_one) - prediction_sp_corr_one = kw.reshape_signal(prediction_sp_one, prediction_sp_ext_one, HW_goal=HW_sp, LW_goal=LW_sp, tP_goal=None) - prediction_sp_corr_rep5 = kw.repeat_signal(prediction_sp_corr_one, nb=2, na=2) - prediction_sp_corr_rep5_trefHW = kw.ts_to_trefHW(prediction_sp_corr_rep5,HWreftime=is1) - - print(f'reshape_signal DOODTIJ: {current_station}') - prediction_np_one_trefHW = kw.ts_to_trefHW(prediction_np_one) - prediction_np_corr_one = kw.reshape_signal(prediction_np_one, prediction_np_ext_one, HW_goal=HW_np, LW_goal=LW_np, tP_goal=None) - prediction_np_corr_rep5 = kw.repeat_signal(prediction_np_corr_one, nb=2, na=2) - prediction_np_corr_rep5_trefHW = kw.ts_to_trefHW(prediction_np_corr_rep5,HWreftime=in1) - - - #12u25m timeseries for BOI computations (no relation between HW and moon, HW has to come at same time for av/sp/np tide, HW timing does differ between stations) - print(f'reshape_signal BOI GEMGETIJ and write to csv: {current_station}') - prediction_av_corrBOI_one = kw.reshape_signal(prediction_av_one, prediction_av_ext_one, HW_goal=HW_av, LW_goal=LW_av, tP_goal=pd.Timedelta(hours=12,minutes=25)) - prediction_av_corrBOI_one_roundtime = prediction_av_corrBOI_one.resample(f'{pred_freq_sec}s').nearest() - prediction_av_corrBOI_one_roundtime.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_av_corrBOI_repn_roundtime = kw.repeat_signal(prediction_av_corrBOI_one_roundtime, nb=0, na=10) - - print(f'reshape_signal BOI SPRINGTIJ and write to csv: {current_station}') - prediction_sp_corrBOI_one = kw.reshape_signal(prediction_sp_one, prediction_sp_ext_one, HW_goal=HW_sp, LW_goal=LW_sp, tP_goal=pd.Timedelta(hours=12,minutes=25)) - prediction_sp_corrBOI_one.index = prediction_sp_corrBOI_one.index - prediction_sp_corrBOI_one.index[0] + prediction_av_corrBOI_one.index[0] #shift times to first HW from gemgetij - prediction_sp_corrBOI_one_roundtime = prediction_sp_corrBOI_one.resample(f'{pred_freq_sec}s').nearest() - prediction_sp_corrBOI_one_roundtime.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_sp_corrBOI_repn_roundtime = kw.repeat_signal(prediction_sp_corrBOI_one_roundtime, nb=0, na=10) - - print(f'reshape_signal BOI DOODTIJ and write to csv: {current_station}') - prediction_np_corrBOI_one = kw.reshape_signal(prediction_np_one, prediction_np_ext_one, HW_goal=HW_np, LW_goal=LW_np, tP_goal=pd.Timedelta(hours=12,minutes=25)) - prediction_np_corrBOI_one.index = prediction_np_corrBOI_one.index - prediction_np_corrBOI_one.index[0] + prediction_av_corrBOI_one.index[0] #shift times to first HW from gemgetij - prediction_np_corrBOI_one_roundtime = prediction_np_corrBOI_one.resample(f'{pred_freq_sec}s').nearest() - prediction_np_corrBOI_one_roundtime.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') - prediction_np_corrBOI_repn_roundtime = kw.repeat_signal(prediction_np_corrBOI_one_roundtime, nb=0, na=10) + + # derive getijkrommes: raw, scaled to havengetallen, scaled to havengetallen and 12h25min period + prediction_av, prediction_sp, prediction_np = kw.gemiddeld_getij_av_sp_np( + df_meas=data_pd_meas_10y, pred_freq=pred_freq, nb=0, nf=0, + scale_extremes=False, scale_period=False) + prediction_av_corr, prediction_sp_corr, prediction_np_corr = kw.gemiddeld_getij_av_sp_np( + df_meas=data_pd_meas_10y, pred_freq=pred_freq, nb=2, nf=2, + scale_extremes=file_havget, scale_period=False) + prediction_av_corr_boi, prediction_sp_corr_boi, prediction_np_corr_boi = kw.gemiddeld_getij_av_sp_np( + df_meas=data_pd_meas_10y, pred_freq=pred_freq, nb=0, nf=10, + scale_extremes=file_havget, 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') cmap = plt.get_cmap("tab10") @@ -350,15 +246,16 @@ 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}') - ax_sum.plot(prediction_av_one_trefHW['values'],'--', color=cmap(0), linewidth=0.7, label='gem kromme, one') - ax_sum.plot(prediction_av_corr_rep5_trefHW['values'], color=cmap(0), label='gem kromme, corr') - ax_sum.plot(prediction_sp_one_trefHW['values'],'--', color=cmap(1), linewidth=0.7, label='sp kromme, one') - ax_sum.plot(prediction_sp_corr_rep5_trefHW['values'], color=cmap(1), label='sp kromme, corr') - ax_sum.plot(prediction_np_one_trefHW['values'],'--', color=cmap(2), linewidth=0.7, label='np kromme, one') - ax_sum.plot(prediction_np_corr_rep5_trefHW['values'], color=cmap(2), label='np kromme, corr') + prediction_av['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['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['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(-15.5,15.5) + 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.savefig(os.path.join(dir_gemgetij,f'gemgetij_trefHW_{current_station}')) @@ -366,11 +263,11 @@ 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 - ax1_boi.plot(prediction_av_corrBOI_repn_roundtime['values'],color=cmap(0),label='prediction gemtij') - ax1_boi.plot(prediction_sp_corrBOI_repn_roundtime['values'],color=cmap(1),label='prediction springtij') - ax1_boi.plot(prediction_np_corrBOI_repn_roundtime['values'],color=cmap(2),label='prediction 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 #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' @@ -390,7 +287,7 @@ ax1_boi.grid() ax1_boi.legend(loc=4) ax1_boi.set_xlabel('times since first av HW (start of ts)') - ax1_boi.set_xlim(tstop_dt-dt.timedelta(hours=2),tstop_dt+dt.timedelta(hours=48)) + 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')) diff --git a/kenmerkendewaarden/gemiddeldgetij.py b/kenmerkendewaarden/gemiddeldgetij.py index d6f87ee..18448c3 100644 --- a/kenmerkendewaarden/gemiddeldgetij.py +++ b/kenmerkendewaarden/gemiddeldgetij.py @@ -3,16 +3,165 @@ Computation of gemiddelde getijkromme """ +import os import numpy as np import pandas as pd import hatyan +import logging +from kenmerkendewaarden.tidalindicators import calc_HWLWtidalrange -__all__ = ["get_gemgetij_components", - "ts_to_trefHW", - "reshape_signal", - "repeat_signal", +__all__ = ["gemiddeld_getij_av_sp_np", ] +logger = logging.getLogger(__name__) + + +def gemiddeld_getij_av_sp_np(df_meas, pred_freq="60sec", nb=0, nf=0, scale_extremes=False, scale_period=False, debug=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 + are then used to make a raw prediction for average/spring/neap tide. + These raw predictions can optionally be scaled in height (with havengetallen) + and in time (to a fixed period of 12h25min). An n-number of backwards and forward repeats + are added before the timeseries are returned, resulting in nb+nf+1 tidal periods. + + df_meas: timeseries of 10 years of waterlevel measurements + pred_freq: frequency of the prediction, a value of 60 seconds or lower is adivisable for decent results. + nb: amount of periods to repeat backward + nf: amount of periods to repeat forward + scale_extremes: scale extremes with havengetallen. Should be a boolean, but is now a filepath to the havengetallen csv files + scale_period: scale to 12h25min (for boi) + """ + data_pd_meas_10y = df_meas + tstop_dt = df_meas.index.max() + + current_station = data_pd_meas_10y.attrs["station"] + + # TODO: deprecate debug argument+plot (maybe use max HW instead of max tidalrange?) + # TODO: we now call this function three times and deriving the unscaled krommes takes quite some time. Put in different function and cache it. + # TODO: add correctie havengetallen HW/LW av/sp/np met slotgemiddelde uit PLSS/modelfit (HW/LW av) + + if scale_period: + tP_goal = pd.Timedelta(hours=12,minutes=25) + else: + tP_goal = None + + # TODO: derive havengetallen on the fly instead, no hardcoded files at least (and deprecate file_havget argument) + if scale_extremes: # if not None, so e.g. file_havget + file_havget = scale_extremes # TODO: this is a temporary solution to pass file_havget + if not os.path.exists(file_havget): + raise Exception(f'havengetallen file does not exist: {file_havget}') + data_havget = pd.read_csv(file_havget, index_col=0) + for colname in ['HW_delay_median','LW_delay_median','getijperiod_median','duurdaling_median']: + data_havget[colname] = pd.to_timedelta(data_havget[colname]) + + HW_sp, LW_sp = data_havget.loc['spring',['HW_values_median','LW_values_median']] + HW_np, LW_np = data_havget.loc['neap',['HW_values_median','LW_values_median']] + HW_av, LW_av = data_havget.loc['mean',['HW_values_median','LW_values_median']] + else: + HW_av = LW_av = None + HW_sp = LW_sp = None + HW_np = LW_np = None + + #derive components via TA on measured waterlevels + comp_frommeasurements_avg, comp_av = get_gemgetij_components(data_pd_meas_10y) + + times_pred_1mnth = pd.date_range(start=pd.Timestamp(tstop_dt.year,1,1,0,0)-pd.Timedelta(hours=12), end=pd.Timestamp(tstop_dt.year,2,1,0,0), freq=pred_freq) #start 12 hours in advance, to assure also corrected values on desired tstart + comp_av.attrs['nodalfactors'] = False #nodalfactors=False to guarantee repetative signal + comp_av.attrs['fu_alltimes'] = True # TODO: this is not true, but this setting is the default + prediction_av = hatyan.prediction(comp_av, times=times_pred_1mnth) + prediction_av_ext = hatyan.calc_HWLW(ts=prediction_av, calc_HWLW345=False) + + time_firstHW = prediction_av_ext.loc[prediction_av_ext['HWLWcode']==1].index[0] #time of first HW + ia1 = prediction_av_ext.loc[time_firstHW:].index[0] #time of first HW + ia2 = prediction_av_ext.loc[time_firstHW:].index[2] #time of second HW + prediction_av_one = prediction_av.loc[ia1:ia2] + prediction_av_ext_one = prediction_av_ext.loc[ia1:ia2] + + # ============================================================================= + # Hatyan predictie voor 1 jaar met gemiddelde helling maansbaan (voor afleiden spring-doodtijcyclus) >> predictie zonder nodalfactors instead + # ============================================================================= + """ + uit: gemiddelde getijkrommen 1991.0 + Voor de ruwe krommen voor springtij en doodtij is het getij voorspeld voor een jaar met gemiddelde helling maansbaan + met uitsluitend zuivere combinaties van de componenten M2 en S2: + tabel: Gebruikte componenten voor de spring- en doodtijkromme + SM, 3MS2, mu2, M2, S2, 2SM2, 3MS4, M4, MS4, + 4MS6, M6, 2MS6, M8, 3MS8, M10, 4MS10, M12, 5MS12 + + In het aldus gemodelleerde getij is de vorm van iedere getijslag, gegeven de getijfase, identiek. + Vervolgens is aan de hand van de havengetallen een springtij- en een doodtijkromme geselecteerd. + + #NOTE: background on choice of components + #below is different than provided list, these shallow ones are extra: ['S4','2SM6','M7','4MS4','2(MS)8','3M2S10','4M2S12'] + #shallow relations, derive 'zuivere harmonischen van M2 en S2' (this means averaging over eenmaaldaagse componenten, but why is that chosen?) + #adding above extra components or oneday freqs, gives a modulation and therefore there is no repetative signal anymore. Apperantly all components in this list have an integer number of periods in one springneap cycle? + dummy,shallowrel,dummy = hatyan.foreman.get_foreman_shallowrelations() + bool_M2S2only = shallowrel[1].isin([1,2]) & shallowrel[3].isin(['M2','S2']) & shallowrel[5].isin(['M2','S2',np.nan]) & shallowrel.index.isin(const_list_year) + shallowdeps_M2S2 = shallowrel.loc[bool_M2S2only,:5] + print(shallowdeps_M2S2) + """ + components_sn = ['A0','SM','3MS2','MU2','M2','S2','2SM2','3MS4','M4','MS4','4MS6','M6','2MS6','M8','3MS8','M10','4MS10','M12','5MS12'] + + #make prediction with springneap components with nodalfactors=False (alternative for choosing a year with a neutral nodal factor). Using 1yr instead of 1month does not make a difference in min/max tidal range and shape, also because of nodalfactors=False. (when using more components, there is a slight difference) + comp_frommeasurements_avg_sncomp = comp_frommeasurements_avg.loc[components_sn] + comp_frommeasurements_avg_sncomp.attrs['nodalfactors'] = False #nodalfactors=False to make independent on chosen year + comp_frommeasurements_avg_sncomp.attrs['fu_alltimes'] = True # TODO: this is not true, but this setting is the default + prediction_sn = hatyan.prediction(comp_frommeasurements_avg_sncomp, times=times_pred_1mnth) + + prediction_sn_ext = hatyan.calc_HWLW(ts=prediction_sn, calc_HWLW345=False) + + #selecteer getijslag met minimale tidalrange en maximale tidalrange (werd geselecteerd adhv havengetallen in 1991.0 doc) + prediction_sn_ext = calc_HWLWtidalrange(prediction_sn_ext) + + time_TRmax = prediction_sn_ext.loc[prediction_sn_ext['HWLWcode']==1,'tidalrange'].idxmax() + is1 = prediction_sn_ext.loc[time_TRmax:].index[0] + is2 = prediction_sn_ext.loc[time_TRmax:].index[2] + + time_TRmin = prediction_sn_ext.loc[prediction_sn_ext['HWLWcode']==1,'tidalrange'].idxmin() + in1 = prediction_sn_ext.loc[time_TRmin:].index[0] + in2 = prediction_sn_ext.loc[time_TRmin:].index[2] + + #select one tideperiod for springtide and one for neaptide + prediction_sp_one = prediction_sn.loc[is1:is2] + prediction_sp_ext_one = prediction_sn_ext.loc[is1:is2] + prediction_np_one = prediction_sn.loc[in1:in2] + prediction_np_ext_one = prediction_sn_ext.loc[in1:in2] + + # plot selection of neap/spring + if debug: + fig, (ax1,ax2) = hatyan.plot_timeseries(ts=prediction_sn,ts_ext=prediction_sn_ext) + ax1.plot(prediction_sp_one['values'],'r') + ax1.plot(prediction_np_one['values'],'r') + ax1.legend(labels=ax1.get_legend_handles_labels()[1]+['kromme spring','kromme neap'],loc=4) + ax1.set_ylabel('waterstand [m]') + ax1.set_title(f'spring- en doodtijkromme {current_station}') + # fig.savefig(os.path.join(dir_gemgetij,f'springdoodtijkromme_{current_station}_slotgem{year_slotgem}.png')) + + #timeseries for gele boekje (av/sp/np have different lengths, time is relative to HW of av and HW of sp/np are shifted there) #TODO: is this product still necessary? + logger.info(f'reshape_signal GEMGETIJ: {current_station}') + prediction_av_corr_one = reshape_signal(prediction_av_one, prediction_av_ext_one, HW_goal=HW_av, LW_goal=LW_av, tP_goal=tP_goal) + prediction_av_corr_one.index = prediction_av_corr_one.index - prediction_av_corr_one.index[0] # make relative to first timestamp (=HW) + if scale_period: # resampling required because of scaling + prediction_av_corr_one = prediction_av_corr_one.resample(pred_freq).nearest() + prediction_av = repeat_signal(prediction_av_corr_one, nb=nb, nf=nf) + + logger.info(f'reshape_signal SPRINGTIJ: {current_station}') + prediction_sp_corr_one = reshape_signal(prediction_sp_one, prediction_sp_ext_one, HW_goal=HW_sp, LW_goal=LW_sp, tP_goal=tP_goal) + prediction_sp_corr_one.index = prediction_sp_corr_one.index - prediction_sp_corr_one.index[0] # make relative to first timestamp (=HW) + if scale_period: # resampling required because of scaling + prediction_sp_corr_one = prediction_sp_corr_one.resample(pred_freq).nearest() + prediction_sp = repeat_signal(prediction_sp_corr_one, nb=nb, nf=nf) + + logger.info(f'reshape_signal DOODTIJ: {current_station}') + prediction_np_corr_one = reshape_signal(prediction_np_one, prediction_np_ext_one, HW_goal=HW_np, LW_goal=LW_np, tP_goal=tP_goal) + prediction_np_corr_one.index = prediction_np_corr_one.index - prediction_np_corr_one.index[0] # make relative to first timestamp (=HW) + if scale_period: # resampling required because of scaling + prediction_np_corr_one = prediction_np_corr_one.resample(pred_freq).nearest() + prediction_np = repeat_signal(prediction_np_corr_one, nb=nb, nf=nf) + + return prediction_av, prediction_sp, prediction_np + def get_gemgetij_components(data_pd_meas): # ============================================================================= @@ -70,8 +219,8 @@ def get_gemgetij_components(data_pd_meas): comp_av.loc['A0'] = comp_frommeasurements_avg.loc['A0'] - print('verhouding tussen originele en kwadratensom componenten') - print(comp_av/comp_frommeasurements_avg.loc[components_av]) # values are different than 1991.0 document and differs per station while the document states "Zoals te verwachten is de verhouding per component tussen deze wortel en de oorspronkelijke amplitude voor alle plaatsen gelijk" + logger.debug('verhouding tussen originele en kwadratensom componenten:\n' + f'{comp_av/comp_frommeasurements_avg.loc[components_av]}') # values are different than 1991.0 document and differs per station while the document states "Zoals te verwachten is de verhouding per component tussen deze wortel en de oorspronkelijke amplitude voor alle plaatsen gelijk" return comp_frommeasurements_avg, comp_av @@ -83,6 +232,12 @@ def reshape_signal(ts, ts_ext, HW_goal, LW_goal, tP_goal=None): time_down was scaled with havengetallen before, but not anymore to avoid issues with aggers """ + # early escape # TODO: consider not calling function in this case + if HW_goal is None and LW_goal is None: + return ts + + # TODO: consider removing the need for ts_ext, it should be possible with min/max, although the HW of the raw timeseries are not exactly equal + TR_goal = HW_goal-LW_goal #selecteer alle hoogwaters en opvolgende laagwaters @@ -121,26 +276,13 @@ def reshape_signal(ts, ts_ext, HW_goal, LW_goal, tP_goal=None): return ts_corr -def ts_to_trefHW(ts,HWreftime=None): - """ - converts to hours relative to HWreftime, to plot av/sp/np tidal signals in one plot - """ - ts.index.name = 'times' #just to be sure - ts_trefHW = ts.reset_index() - if HWreftime is None: - ts_trefHW.index = (ts_trefHW['times']-ts_trefHW['times'].iloc[0]).dt.total_seconds()/3600 - else: - ts_trefHW.index = (ts_trefHW['times']-HWreftime).dt.total_seconds()/3600 - return ts_trefHW - - -def repeat_signal(ts_one_HWtoHW, nb, na): +def repeat_signal(ts_one_HWtoHW, nb, nf): """ repeat tidal signal, necessary for sp/np, since they are computed as single tidal signal first """ tidalperiod = ts_one_HWtoHW.index[-1] - ts_one_HWtoHW.index[0] ts_rep = pd.DataFrame() - for iAdd in np.arange(-nb,na+1): + for iAdd in np.arange(-nb,nf+1): ts_add = pd.DataFrame({'values':ts_one_HWtoHW['values'].values}, index=ts_one_HWtoHW.index + iAdd*tidalperiod) ts_rep = pd.concat([ts_rep,ts_add])