Pull request - update extremes app

ruins/apps/extremes.py

@@ -3,6 +3,7 @@
 import pandas as pd
 import numpy as np
 import datetime
+import pickle
 from plotly.subplots import make_subplots
 
 from ruins.core import build_config, debug_view, DataManager, Config
@@ -67,25 +68,24 @@
 """
 )
 
+#load cached events
+with open('cache/events.pkl', 'rb') as file:
+    events = pickle.load(file)
+with open('cache/canals.pkl', 'rb') as file:
+    canal_par = pickle.load(file)
 
 def user_input_defaults():
     # streamlit input stuff:
 
     slr = 400   # sea level rise in mm (0, 400, 800, 1200, 1600)
-
-    #recharge_vis = "absolute"   # "cumulative" or "absolute"
+    prec_increase = 1 #Intensify precipitation by factor
 
     # default event z.B.:
-    # (wählt Jonas noch aus)
+    time = list(events.keys())[1]
 
     #if time == "2012":
     t1 = datetime.date(2011, 12, 28)
     t2 = datetime.date(2012, 1, 12)
-
-
-    ## KGE
-    # kge = st.slider("Canal flow uncertainty [KGE * 100]",71,74, value = 74)
-    kge = 74 # nicht mehr user input
 
     ## canal flow input
     # canal_flow_scale = st.number_input("Factor to canal capacity", min_value=0.5, max_value=3., value= 1.0, step=0.1) 
@@ -94,12 +94,12 @@ def user_input_defaults():
 
     canal_area = 4 # user input: st.radio(
         #"Share of water area on catchment [%].",
-        #(3, 4))
+        #(4, 6))
 
     ## pump before event
     advance_pump = 0. # user input: st.radio(
         #"Forecast Pumping",
-        #(0, 4)
+        #(0, 50)
 
     ## visualisation of used pump capacity
     # pump_vis = st.radio("Pump capacity visualisation", ["absolute", "cumulative"])
@@ -108,52 +108,47 @@ def user_input_defaults():
     # maxdh = st.number_input("Stop pumping if dh at Knock is greater than x dm\n(technical limit = 70dm)", min_value=10, max_value=70, value= 40, step=2) 
     maxdh = 4000 # nicht mehr user input
 
-    return slr, t1, t2, kge, canal_flow_scale, canal_area, advance_pump, maxdh
+    return slr, t1, t2, canal_flow_scale, canal_area, advance_pump, maxdh, prec_increase
 
 
-def timeslice_observed_data(dataManager: DataManager, t1, t2, slr):
-
-    raw = dataManager['levelknock'].read()
-    weather_1h = dataManager['prec'].read()
-
+def timeslice_observed_data(dataManager: DataManager, t1, t2, slr, prec_increase, prec_line):
+    
+    extremes = dataManager['hydro_krummh'].read()
+    # recharge
+    hourly_recharge = extremes[prec_line].to_dataframe().rolling("12h").mean()[t1:t2].squeeze() * prec_increase # changed by Jonas
     # tide data
-    tide = raw['L011_pval'][t1:t2]*1000 + slr
+    tide = (extremes['wl_Knock_Outer'].to_dataframe()[t1:t2] + slr).squeeze()
+    # water level
+    EVEx5_lw_pegel_timesliced = (extremes['wl_LW'].to_dataframe()[t1:t2]/1000).squeeze()
 
-    # hourly recharge data
-    hourly_recharge = weather_1h["Prec"][t1:t2]
-    hourly_recharge = hourly_recharge.rolling("12h").mean() # changed by Jonas
-
-    # EVEx5 observed
-
-    #EVEx5 = pd.read_csv('//home/lucfr/hydrocode/RUINS_hydropaper-newlayout/streamlit/data/levelLW.csv')
-    EVEx5 = dataManager['levelLW'].read()
-    EVEx5.index = pd.to_datetime(EVEx5.iloc[:,0])
-    EVEx5_lw_pegel_timesliced = (EVEx5['LW_Pegel_Aussen_pval']/100+0.095)[t1:t2]
-
-    # pump observed
-    pump_capacity_observed = raw['sumI'][t1:t2] / 12.30
+    pump_capacity_observed = extremes['Knock_pump_obs'].to_dataframe()[t1:t2].squeeze()
 
+    # .squeeze() convertes single pd.DataFrame columns to pd.series objects
+
     return tide, hourly_recharge, EVEx5_lw_pegel_timesliced, pump_capacity_observed
 
 
 def create_initial_x_dataset(tide_data, hourly_recharge):
 
-    wig = tide_data*0
-    # what is this and why? Unused!
-
-    x = pd.DataFrame.from_dict({'recharge' : hourly_recharge,
-                                'h_tide' : tide_data,
-                                'wig' : wig})
+    x = pd.DataFrame()
+    x['h_tide'] = tide_data
+    x['recharge'] = hourly_recharge
+    x['wig'] = 0.  # what is this and why? Unused! # comment Jonas: Experimental option to account for a "Wind Induced Gradient" in the canals, which reduces the effective water flow gradient in the drain_cap model
 
     return x 
 
 
-def create_model_runs_list(all_kge_canal_par_df, kge, canal_flow_scale, canal_area, x_df, advance_pump, maxdh):    
+def create_model_runs_list(canal_flow_scale, canal_area, x_df, advance_pump, maxdh, canal_par_array):
     model_runs = []
 
-    kge_canal_par_df = all_kge_canal_par_df.loc[all_kge_canal_par_df.KGE == kge]
-    canal_par_array = kge_canal_par_df[['parexponent', 'parfactor']].to_numpy()
-
+    ### Canal flow parameters from fitting of runoff to canal gradient data
+    #canal_par_array = [[1.112,4156.],[1.045 , 2820.],[0.9946,2142.]]
+
+    ### Define parameters of the default pumping function / pump chart
+    x = np.array([7.,6.,5.,4.,3.5,3.,2,1,0,5.,4.,3.5,3.,2]) * 1000 # water gradient [mm] (by factor 1000 from m)
+    y = np.array([0,4.2,8.4,12.6,14.5,15.8,17.5,19,20.5,8.4,12.6,14.5,15.8,17.5]) * 3600 / (35000 * 100 * 100) * 1000 * 4  # "*3600 / (35000 * 100 * 100) * 1000 * 4)" converts m^3/s in mm/h
+    pumpcap_fit = np.polynomial.polynomial.Polynomial.fit(x = x, y = y, deg = 2)
+
     for z in canal_par_array:
 
         z[1] /= canal_flow_scale
@@ -163,13 +158,15 @@ def create_model_runs_list(all_kge_canal_par_df, kge, canal_flow_scale, canal_ar
          q_pump,
          x_df['h_min'], 
          x_df['flow_rec'], 
-         pump_cost) = drain_cap.storage_model(x_df,
-                                              z, 
-                                              storage = 0, 
-                                              h_store = -1350, # geändert von Jonas
+         pump_cost,
+         store_vol) = drain_cap.storage_model(forcing_data = x_df,
+                                              canal_par = z, 
+                                              v_store = 0, 
+                                              h_store_target = -1350, # geändert von Jonas
                                               canal_area = canal_area, 
-                                              advance_pump = advance_pump, 
-                                              maxdh = maxdh)
+                                              h_forecast_pump = advance_pump, 
+                                              h_grad_pump_max = maxdh,
+                                              pump_par = pumpcap_fit)
         model_runs.append((x_df['h_store'], pump_cost))
         #pump_capacity_model_runs.append(pump_cost)
 
@@ -189,57 +186,84 @@ def flood_model(dataManager: DataManager, config:Config, **kwargs):
     st.sidebar.header('Control Panel')
 
 
-    slr, t1, t2, kge, canal_flow_scale, canal_area, advance_pump, maxdh = user_input_defaults()
+    slr, t1, t2, canal_flow_scale, canal_area, advance_pump, maxdh, prec_increase = user_input_defaults()
 
     with st.sidebar.expander("Event selection"):
         time = st.radio(
             "Event",
-            ("2012", "2017", "choose custom period")    # reduced to 2 nice events -> "custom period" only in expert mode or for self hosting users?
+            (events.keys())
         )
-        if time == 'choose custom period':
-            t1 = st.date_input("start", datetime.date(2017, 12, 1))
-            dt2 = st.number_input("Number of days", min_value=3, max_value=20, value= 10, step=1) 
-            t2 = t1 + datetime.timedelta(dt2)
+    t1 = events[time]
+    t2 = events[time]+datetime.timedelta(days=14)
 
     with st.sidebar.expander("Sea level rise"):
         slr = st.radio(
             "Set SLR [mm]",
-            (0, 400, 800, 1200, 1600)
+            (0,154,249,379,432,522,730,848,918,1143,1676)
+        )
+
+    with st.sidebar.expander("Precipitation"):
+        prec_increase = st.radio(
+            "Intensify precipitation by factor",
+            (0.8, 0.9, 1, 1.1, 1.2, 1.3)
+        )
+
+        prec_line = st.selectbox(
+            "Choose precipitation realization",
+            ('Prec', 'Prec_dissagg_1', 'Prec_dissagg_2', 'Prec_dissagg_3', 'Prec_dissagg_4',
+           'Prec_dissagg_5', 'Prec_dissagg_6', 'Prec_dissagg_7', 'Prec_dissagg_8',
+           'Prec_dissagg_9', 'Prec_dissagg_10', 'Prec_dissagg_11',
+           'Prec_dissagg_12', 'Prec_dissagg_13', 'Prec_dissagg_14',
+           'Prec_dissagg_15', 'Prec_dissagg_16', 'Prec_dissagg_17',
+           'Prec_dissagg_18', 'Prec_dissagg_19', 'Prec_dissagg_20',
+           'Prec_dissagg_21', 'Prec_dissagg_22', 'Prec_dissagg_23',
+           'Prec_dissagg_24', 'Prec_dissagg_25', 'Prec_dissagg_26',
+           'Prec_dissagg_27', 'Prec_dissagg_28', 'Prec_dissagg_29',
+           'Prec_dissagg_30', 'Prec_dissagg_31', 'Prec_dissagg_32',
+           'Prec_dissagg_33', 'Prec_dissagg_34', 'Prec_dissagg_35',
+           'Prec_dissagg_36', 'Prec_dissagg_37', 'Prec_dissagg_38',
+           'Prec_dissagg_39', 'Prec_dissagg_40', 'Prec_dissagg_41',
+           'Prec_dissagg_42', 'Prec_dissagg_43', 'Prec_dissagg_44',
+           'Prec_dissagg_45', 'Prec_dissagg_46', 'Prec_dissagg_47',
+           'Prec_dissagg_48', 'Prec_dissagg_49', 'Prec_dissagg_50',
+           'Prec_dissagg_51', 'Prec_dissagg_52', 'Prec_dissagg_53',
+           'Prec_dissagg_54', 'Prec_dissagg_55', 'Prec_dissagg_56',
+           'Prec_dissagg_57', 'Prec_dissagg_58', 'Prec_dissagg_59',
+           'Prec_dissagg_60', 'Prec_dissagg_61', 'Prec_dissagg_62',
+           'Prec_dissagg_63', 'Prec_dissagg_64', 'Prec_dissagg_65',
+           'Prec_dissagg_66', 'Prec_dissagg_67', 'Prec_dissagg_68',
+           'Prec_dissagg_69', 'Prec_dissagg_70', 'Prec_dissagg_71',
+           'Prec_dissagg_72', 'Prec_dissagg_73', 'Prec_dissagg_74',
+           'Prec_dissagg_75', 'Prec_dissagg_76', 'Prec_dissagg_77',
+           'Prec_dissagg_78', 'Prec_dissagg_79', 'Prec_dissagg_80')
         )
-
-    if time == "2012":
-        t1 = datetime.date(2011, 12, 28)
-        t2 = datetime.date(2012, 1, 12)
 
-    if time == "2017":
-        t1 = datetime.date(2017, 3, 15)
-        t2 = datetime.date(2017, 3, 25)
-
     with st.sidebar.expander("Management options"):
     # pump before event
     #    advance_pump = st.number_input("Additional spare volume in canals", min_value=-5., max_value=8., value= 0., step=0.1)
         advance_pump = st.radio(
-            "Forecast Pumping",
-            (0, 4)
+            "Lower water level by x mm NHN before event.",
+            (0, 50, 200)
         )
-        Canal_area = st.radio(
+        canal_area = st.radio(
             "Share of water area on catchment [%].",
-            (3, 4)
+            (4, 6)
+        )
+        maxdh = st.radio(
+            "Maximum pump gradient [mm]",
+            (1, 4500, 6000)
         )
 
     # Model runs:
 
     (tide, 
     hourly_recharge, 
     EVEx5_lw_pegel_timesliced, 
-    pump_capacity_observed) = timeslice_observed_data(dataManager, t1, t2, slr)
+    pump_capacity_observed) = timeslice_observed_data(dataManager, t1, t2, slr, prec_increase, prec_line)
 
     x = create_initial_x_dataset(tide, hourly_recharge)
 
-    all_kge_canal_par_df = dataManager['kge_canal_par'].read()
-
-
-    hg_model_runs = create_model_runs_list(all_kge_canal_par_df, kge, canal_flow_scale, canal_area, x, advance_pump, maxdh)
+    hg_model_runs = create_model_runs_list(canal_flow_scale, canal_area, x, advance_pump, maxdh, canal_par)
 
     # plotting:
     col1, col2 = container.columns(2)