search

Deltares-research / kenmerkendewaarden

Derive indicators from waterlevel measurements
https://deltares-research.github.io/kenmerkendewaarden/
GNU General Public License v3.0
1 stars 0 forks source link

Mean waterlevels are different than in psmsl #110

Open veenstrajelmer opened 2 weeks ago

veenstrajelmer commented 2 weeks ago

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

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)
    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')
    df_psmsl['mean_wl'] /= 1000 # mm to meters
    df_psmsl['mean_wl'] -= 6.987 + 0.114 # RLR to NAP from https://psmsl.org/data/obtaining/rlr.diagrams/22.php

    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')/100
        yearmeanLW = pd.read_csv(file_yearmeanLW, sep='\\s+', skiprows=1, names=['datetime','values'], parse_dates=['datetime'], na_values=-999.9, index_col='datetime')/100
        yearmeanwl = pd.read_csv(file_yearmeanwl, sep='\\s+', skiprows=1, names=['datetime','values'], parse_dates=['datetime'], na_values=-999.9, index_col='datetime')/100
        # 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)
    fig.tight_layout()

Gives: image