Mean waterlevels are different compared to psmsl

[veenstrajelmer]

Should be corrected from RLR to nap. RLR should contain a timeseries with continuous vertical reference.

• nap hoogtes van peilbout worden periodiek ingemeten, staat allemaal in nap-info: https://maps.rijkswaterstaat.nl/geoweb55/index.html?viewer=NAPinfo. Toegang vragen kan, dan krijg je toegang tot extra (groene) puntjes. >> we cannot get this data automated for all stations
• FB also compared psmsl metric to rlr: https://github.com/openearth/sealevel/blob/report/notebooks/validation/metric-vs-rlr.ipynb
• FB also compared waterdata/ddldata to psmsl: https://github.com/openearth/sealevel/blob/report/notebooks/validation/paired%20difference%20analysis%20RWS%20and%20PSMSL.ipynb
• FB/WS use this for correction: https://github.com/Deltares-research/sealevelmonitor/blob/main/data/deltares/main_stations.json

Some code to test:

import os
import pandas as pd
import matplotlib.pyplot as plt
plt.close('all')
import hatyan
import kenmerkendewaarden as kw # pip install git+https://github.com/Deltares-research/kenmerkendewaarden

dir_base = r'p:\11210325-005-kenmerkende-waarden\work'
dir_meas = os.path.join(dir_base,'measurements_wl_18700101_20240101')

station_list = ["HOEKVHLD"]

nap_correction = False
min_coverage = 0.9 # for tidalindicators and slotgemiddelde #TODO: can also be used for havengetallen and gemgetij

for current_station in station_list:
    # timeseries are used for slotgemiddelden, gemgetijkrommen (needs slotgem+havget)
    data_pd_meas_all = kw.read_measurements(dir_output=dir_meas, station=current_station, extremes=False, nap_correction=nap_correction)
    mean_wl_ddl = kw.calc_wltidalindicators(data_pd_meas_all, min_coverage=min_coverage)['wl_mean_peryear']
    mean_wl_ddl.index = mean_wl_ddl.index.to_timestamp()
    
    dir_data = r"p:\archivedprojects\11208031-010-kenmerkende-waarden-k\work\data_vanRWS_20220805\wetransfer_waterstandsgegevens_2022-08-05_1306"
    file_dia_hoek1 = os.path.join(dir_data, r"WATHTE_10min\HOEKVHLD_1.dia")
    file_dia_hoek2 = os.path.join(dir_data, r"WATHTE_oud\HOEK_KW.dia")
    data_wl_tkdia = hatyan.read_dia([file_dia_hoek1,file_dia_hoek2], block_ids="allstation", station="HOEKVHLD", allow_duplicates=True)
    data_wl_tkdia['values'] *= 100
    mean_wl_dia = kw.calc_wltidalindicators(data_wl_tkdia, min_coverage=min_coverage)['wl_mean_peryear']
    mean_wl_dia.index = mean_wl_dia.index.to_timestamp()

    # for HOEKVHLD
    df_psmsl = pd.read_csv("https://psmsl.org/data/obtaining/rlr.annual.data/22.rlrdata", sep=';', names=['year','mean_wl','flag','code'])
    df_psmsl.index = pd.PeriodIndex(df_psmsl['year'], freq='Y').to_timestamp()
    df_psmsl['mean_wl'] /= 10 # mm to centimeters
    # RLR to NAP from https://psmsl.org/data/obtaining/rlr.diagrams/22.php
    nap_correctie = df_psmsl['mean_wl'] * 0 + 4.855*100
    bool_2005 = nap_correctie.index >= "2005-01-01"
    nap_correctie.loc[bool_2005] = 4.827*100
    df_psmsl['mean_wl'] = df_psmsl['mean_wl'] - 11.7*100 + nap_correctie
    
    fig, ax1 = plt.subplots(figsize=(12,7))
    ax1.plot(mean_wl_ddl, 'rx', label='yearmean ddl')
    ax1.plot(mean_wl_dia, 'k2', label='yearmean dia')
    ax1.plot(df_psmsl['mean_wl'], 'c1', label='yearmean psmsl rlr to nap')
    
    station_name_dict = {'HOEKVHLD':'hoek',
                         'HARVT10':'ha10'}
    if current_station in station_name_dict.keys():
        dir_meas_gemHWLWwlAB = r'p:\archivedprojects\11208031-010-kenmerkende-waarden-k\work\data_KW-RMM'
        file_yearmeanHW = os.path.join(dir_meas_gemHWLWwlAB,f'{station_name_dict[current_station]}_hw.txt')
        file_yearmeanLW = os.path.join(dir_meas_gemHWLWwlAB,f'{station_name_dict[current_station]}_lw.txt')
        file_yearmeanwl = os.path.join(dir_meas_gemHWLWwlAB,f'{station_name_dict[current_station]}_Z.txt')
        yearmeanHW = pd.read_csv(file_yearmeanHW, sep='\\s+', skiprows=1, names=['datetime','values'], parse_dates=['datetime'], na_values=-999.9, index_col='datetime')
        yearmeanLW = pd.read_csv(file_yearmeanLW, sep='\\s+', skiprows=1, names=['datetime','values'], parse_dates=['datetime'], na_values=-999.9, index_col='datetime')
        yearmeanwl = pd.read_csv(file_yearmeanwl, sep='\\s+', skiprows=1, names=['datetime','values'], parse_dates=['datetime'], na_values=-999.9, index_col='datetime')
        # ax1.plot(yearmeanHW['values'],'+g', zorder=0)
        # ax1.plot(yearmeanLW['values'],'+g', zorder=0)
        ax1.plot(yearmeanwl['values'],'+g',label='yearmean validation', zorder=0)
    
    ax1.legend(loc=2)
    ax1.grid()
    fig.tight_layout()

Gives: