diff --git a/cache/events.pkl b/cache/events.pkl new file mode 100644 index 0000000..3105f8f Binary files /dev/null and b/cache/events.pkl differ diff --git a/ruins/apps/extremes.py b/ruins/apps/extremes.py index 67f14a5..f574a0b 100644 --- a/ruins/apps/extremes.py +++ b/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,42 +186,72 @@ 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: @@ -232,14 +259,11 @@ def flood_model(dataManager: DataManager, config:Config, **kwargs): (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) diff --git a/ruins/processing/drain_cap.py b/ruins/processing/drain_cap.py index 1c3f6c5..b32f834 100644 --- a/ruins/processing/drain_cap.py +++ b/ruins/processing/drain_cap.py @@ -2,16 +2,17 @@ import numpy as np import pandas as pd -x = np.array([7.,6.,5.,4.,3.5,3.,2,1,0,5.,4.,3.5,3.,2]) -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]) +### 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) -def drain_cap(h_tide: np.ndarray, h_store: np.ndarray, h_min: int = -2000, pump_par = pumpcap_fit, channel_par: Tuple[float, float] = [1.016 , 2572.], dh: int = 50, wind_safe: int = 0, gradmax: int = 4000): +def drain_cap(h_tide: np.ndarray, h_store: np.ndarray, h_min: int = -2000, pump_par = pumpcap_fit, canal_par: Tuple[float, float] = [1.016 , 2572.], h_increment: int = 50, h_wind_safe: int = 0, h_grad_pump_max: int = 4000): """ - Find sthe maximal flow rate in a system with a pump and a canal. + Find the maximal flow rate in a system with a pump and a canal. The flow through the pump is determined by the gradient from inner to outer water level and a pump function. - The flow through the canal is determined by the gradient from the canal water level to the pump inner water level and a flow function. + The flow through the canal is determined by the gradient from the canal water level to the pump inner water level and a canal flow function. The inner water level is unknown and estimated within this function, but a lower limit can be set. Parameters @@ -24,18 +25,18 @@ def drain_cap(h_tide: np.ndarray, h_store: np.ndarray, h_min: int = -2000, pump_ the lower boundary of the inner water level pump_par : np.ndarray parameters of the pump function - channel_par : Tuple[float, float] + canal_par : Tuple[float, float] parameters of the canal flow function - dh : int + h_increment : int increment of inner water level estimation - wind_safe : int + h_wind_safe : int gradient from canal to inner water level, which is induced by wind and therefore not contributing to flow - gradmax :int + h_grad_pump_max :int maximum gradient from inner to outer water level, at which pumps shall run Returns ------- - a_channel : np.ndarry + q_channel : np.ndarray the actual maximum flow through canal and pump h_min : float the estimated inner water level @@ -44,73 +45,110 @@ def drain_cap(h_tide: np.ndarray, h_store: np.ndarray, h_min: int = -2000, pump_ """ # set h_min to either absolute technical lower limit or maximal pump gradient - h_min = np.maximum(h_min, h_tide - gradmax) + h_min = np.maximum(h_min, h_tide - h_grad_pump_max) pumplim = True # in case drainage ist limited by pumps, the minimum water level at knock increases until the flow is limited by channel_ while pumplim: - if(h_tide <= h_min): - q_pump_m3 = pump_par((1)/1000) # assume 1 mm gradient to pump to get some estimate - eventually sluicing is more effective - q_pump = q_pump_m3 * 3600 / (35000 * 100 * 100) * 1000 * 4 # "*3600 / (35000 * 100 * 100) * 1000 * 4)" converts m^3/s in mm/h + if(h_tide <= h_min): # if outer water level is lower then inner water level we can sluice + q_pump = pump_par(1) # assume 1 mm gradient to pump to get some estimate - eventually sluicing is more effective else: - q_pump_m3 = pump_par((h_tide - h_min)/1000) - q_pump = q_pump_m3 *3600 / (35000 * 100 * 100) * 1000 * 4 # "*3600 / (35000 * 100 * 100) * 1000 * 4)" converts m^3/s in mm/h + q_pump = pump_par(h_tide - h_min) - if(((h_store - h_min) - wind_safe) <= 0): + if(((h_store - h_min) - h_wind_safe) <= 0): q_channel = 0. else: - q_channel = ((((h_store - h_min) - wind_safe)**channel_par[0])/channel_par[1]) + q_channel = ((((h_store - h_min) - h_wind_safe)**canal_par[0])/canal_par[1]) # end the loop, if h_min is set to value that q_channel is smaller than q_pump, otherwise increase h_min if(q_pump >= q_channel or q_channel <= 0) : pumplim = False else: - h_min += dh + h_min += h_increment # set output h_min to either technical lower limit or water level in canals h_min = np.minimum(h_min, h_store) return (q_channel, h_min, q_pump) -def storage_model (x, z, storage = 0, h_store = -1400, canal_area = 4, advance_pump = 0, maxdh = 6000): +def storage_model (forcing_data, canal_par, v_store = 0, h_store_target = -1400, canal_area = 4, h_forecast_pump = 0, h_grad_pump_max = 6000, h_canal_max = -900, pump_par = pumpcap_fit, h_min = -2000): """ - Storage model used for the Krummhörn region + Storage model used for the Krummhoern region + + Parameters + ---------- + forcing_data : pd.DataFrame + time series forcing data of the storage model with the columns: + - 'recharge' the water volume [waterbalace mm] which flows into the storage + - 'h_tide' the outer (tidal) water level at the pumps [mm NHN] + - 'wig' a wind induced gradient [mm height] in the canals, which alters the flow capacity of the canals + canal_par : np.ndarray + parameters of the canal flow capacity + v_store : numeric + initial water volume [waterbalace mm] in the canals above target water level + h_store_target : numeric + target water level [mm NHN] in the canals + canal_area : numeric + share of water area on total landscape area [%] + h_forecast_pump : numeric + water level [mm height] which is substracted from target water level to allow lower canal water level before forecasted event + h_grad_pump_max : numeric + allowed maximal gradient [h_tide - h_min] for pumps, at which pumps will be turned off for their technical protection + pump_par : np.ndarray + parameters of the pump function + h_min : numeric + minimum water level height [mm NHN] at inner side of pumps + + Returns + ------- + h_store_rec : np.ndarray + water height in canals [mm NHN] + q_pump_rec : np.ndarray + possible drainage capacity by pumps [waterbalace mm] + h_min_rec : np.ndarray + evaluated lowest inner water level [mm NHN] + q_rec : np.ndarray + "real" flow [waterbalace mm], either limited by drainage capacity or pump capacity + usage_pump_rec : np.ndarray + Actually used percentage of possible pumping capacity + v_store_rec : np:ndarray + time series of water volume in storage [waterbalace mm] """ - store = [] - h_min = [] - q_pump = [] - pump_cost = [] - flow_rec = [] + v_store_rec = [] + h_min_rec = [] + q_pump_rec = [] + usage_pump_rec = [] + q_rec = [] - for idx, game in x.iterrows(): + for step_id, step in forcing_data.iterrows(): # recharge storage - storage += game['recharge'] - - # get drain_cap, do not pump if dh > than specified limit -# if(game['h_tide'] - (h_store + storage*100/canal_area)>maxdh): -# cap = [0,h_store + storage*100/canal_area, 0.0000000001] -# else: - cap = drain_cap(game['h_tide'], (h_store + storage*100/canal_area), channel_par = z, dh = 1, wind_safe = game['wig'], gradmax = maxdh) + v_store += step['recharge'] + v_store_max_step = v_store + cap = drain_cap(h_tide = step['h_tide'], # parse tidal water level + h_store = np.min(((h_store_target + v_store*100/canal_area), h_canal_max)), # limit maximal water flow at upper canal crest + h_min = h_min, # parse lowest inner water level + pump_par = pump_par, # parse pump parameters + canal_par = canal_par, # parse canal parameters + h_increment = 1, # increment inner water level by 1 mm steps + h_wind_safe = step['wig'], # parse wind induced gradient + h_grad_pump_max = h_grad_pump_max) # parse maximum pump gradient # drain storage - storage -= cap[0] + v_store -= cap[0] # compare new storage value to lower limit of storage - storage = np.maximum(storage, -advance_pump) + v_store = np.maximum(v_store, -h_forecast_pump/100*canal_area) # save time step - store = np.append(store, storage) - h_min = np.append(h_min, cap[1]) - q_pump = np.append(q_pump, cap[2]) - + v_store_rec = np.append(v_store_rec, v_store) + h_min_rec = np.append(h_min_rec, cap[1]) + q_pump_rec = np.append(q_pump_rec, cap[2]) + + v_flow = v_store_max_step - v_store + q_rec = np.append(q_rec, v_flow) # save "power consumption" of pumps - if(storage > -advance_pump): - flow = cap[0] - else: - flow = game['recharge'] - flow_rec = np.append(flow_rec, flow) - pump_cost = np.append(pump_cost, flow/cap[2]) + usage_pump_rec = np.append(usage_pump_rec, v_flow/cap[2]) - h_store_rec = h_store + store*100/canal_area + h_store_rec = h_store_target + v_store_rec*100/canal_area - return (h_store_rec, q_pump, h_min, flow_rec, pump_cost) + return (h_store_rec, q_pump_rec, h_min_rec, q_rec, usage_pump_rec, v_store_rec)