From 1f889c5c25f863746163681c16237ba7d7025d2c Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Thu, 11 May 2023 11:56:03 +0200 Subject: [PATCH 01/22] First attempt at dynamic update --- model/HBM.py | 124 ++++++++++++++++++++++++++------------------------- 1 file changed, 63 insertions(+), 61 deletions(-) diff --git a/model/HBM.py b/model/HBM.py index 942583c..37c1250 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -19,6 +19,7 @@ from StandardArraysLUE import StandardArraysLUE from RetrieveData import RetrieveData from CalculateFlux import CalculateFlux +from utilityFunctionsHBM import utilityFunctions # Other tools import tools.MakeGIF @@ -72,17 +73,41 @@ def __init__(self, configuration): self.gw_base = float(configuration.modelSettings['groundWaterBase']) * self.standard_LUE.one() # self.ldd = lfr.d8_flow_direction(self.dem) self.ldd = lfr.from_gdal(self.input_dir + configuration.dataSettings['ldd'], partition_shape) - # Set constants + + + #### SETTING CONSTANTS #### + ## GENERAL ## self.resolution = float(configuration.modelSettings['resolution']) self.cell_area = self.resolution * self.resolution self.slope = lfr.slope(self.dem, self.resolution) + self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(self.slope, 0.05, 0.00001) + self.channel_width = 1 + self.coefficient = self.mannings / (self.slope_sqrd * self.width) self.imperm_below_dem = float(configuration.modelSettings['impermeableLayerBelowDEM']) self.imperm_lay_height = self.dem - self.imperm_below_dem self.water_below_dem = float(configuration.modelSettings['waterBelowDEM']) - self.max_gw_s = self.imperm_below_dem * self.cell_area # Full storage of porosity - self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point + self.max_gw_s = self.imperm_below_dem * self.cell_area # Full storage of porosity + self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point + # Refactorings value from mm/hour to m/h times the cell area. + self.refactor = (self.cell_area / 1000) / 3600 + + # Channel length and area + self.channel_length = self.resolution * self.standard_LUE.one() + self.channel_area = self.width * self.channel_length + self.channel_rat = self.channel_area / self.cell_area + self.infil_to_gw_s = self.channel_area / self.porosity # self.notBoundaryCells = generate.boundaryCell() # Currently not working + # Kinematic Surface Water Routing Constants + self.alpha = 1.5 + self.beta = 0.6 + self.timestep = 1.0 * float(configuration.modelSettings['timestep']) + self.c = 5/3 + + # Static, really small value because inflow = 0 is not accepted + self.inflow = self.standard_LUE.one()*1E-20 + + # Load initial groundWaterStorage, if no raster is supplied, use the waterBelowDEM in combination with DEM to create a initialGroundWaterStorage layer. try: self.ini_gw_s = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniGroundWaterStorage'], partition_shape) @@ -109,8 +134,33 @@ def __init__(self, configuration): print("\n") - def update_and_route(self): - pass + def update_and_route(self, gw_s, gw_flux, height, sw_flux): + # The groundwater is adjusted by the fluxes + # channel_infiltation = lfr.where(height > pot_channel_infiltation, pot_channel_infiltation, height) + gw_s = gw_s + gw_flux #+ channel_infiltation*infil_to_gw_s + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) + + # Discharge is affected by the surfacewater fluxes, and seepage is added + height = height + ((sw_flux + seepage)/self.channel_area) #- channel_infiltation + + discharge = lfr.pow(height, self.c) / self.coefficient + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < 1E-20, 1E-20, discharge) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, self.inflow,\ + self.alpha, self.beta, self.timestep,\ + self.channel_length) + + height = lfr.pow(self.constants.coefficient*discharge, 0.6) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + gw_s = gw_s - (seepage / self.porosity) + + return height, discharge, gw_s @lfr.runtime_scope def dynamic_model(self, configuration, report): @@ -119,38 +169,13 @@ def dynamic_model(self, configuration, report): eD = list(map(int, configuration.modelSettings['endDate'].split(", "))) start_date = datetime.datetime(sD[0], sD[1], sD[2], sD[3], sD[4], sD[5]) end_date = datetime.datetime(eD[0], eD[1], eD[2], eD[3], eD[4], eD[5]) - dT = int((end_date - start_date).seconds / dt) + dT = int((end_date - start_date).total_seconds() / dt) # Loading initial conditions height = self.ini_water_h gw_s = self.ini_gw_s int_s = self.ini_int_s - gw_height = self.imperm_lay_height + gw_s/self.cell_area - - # Values for discharge to height calculation - slope_sqrd = lfr.sqrt(self.slope) - slope_sqrd = lfr.where(slope_sqrd < 0.001, 0.001, slope_sqrd) - slope_sqrd = lfr.where(slope_sqrd > 0.05, 0.05, slope_sqrd) - width = 1 - coefficient = self.mannings / (slope_sqrd * width) - - # Channel length and area - channel_length = self.resolution * self.standard_LUE.one() - channel_area = width * channel_length - channel_rat = channel_area / self.cell_area - infil_to_gw_s = channel_area / self.porosity - - # Refactorings value from mm/hour to m/h times the cell area. - refactor = (self.cell_area / 1000) / 3600 - - # Kinematic Surface Water Routing Constants - alpha = 1.5 - beta = 0.6 - timestep = 1.0 * float(configuration.modelSettings['timestep']) - c = 5/3 - - # Static, really small value because inflow = 0 is not accepted - inflow = self.standard_LUE.one()*1E-20 + gw_height = self.imperm_lay_height + gw_s/self.cell_area # Open file to write maximum discharge values to for post simulation validation. with open(self.output_dir + "/maximumDischarge.csv", "w", newline="") as f: @@ -159,16 +184,16 @@ def dynamic_model(self, configuration, report): # Start model for dT large periods for i in range(dT): # Time in minutes is the small iteration multiplied with the timestep (both in seconds) divided by 60 seconds. - date = start_date + datetime.timedelta(seconds = i * (dt*timestep)) + date = start_date + datetime.timedelta(seconds = i * (dt*self.timestep)) # Load flux and storage values precipitation = self.retrieve_data.csv_timeseries_to_flux(configuration.generalSettings['inputDir'] + configuration.dataSettings['precipitationData'], - refactor, date) # m/s + self.refactor, date) # m/s ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(configuration.generalSettings['inputDir'] + configuration.dataSettings['evapotranspirationData'], - refactor, date) # m/s + self.refactor, date) # m/s int_s, precipitation, evapotranspiration_surface = self.calculate_flux.interception(int_s, self.max_int_s, @@ -188,13 +213,13 @@ def dynamic_model(self, configuration, report): evapotranspiration_surface) # The infiltration happens only in the region that is used by the channel and therefore this factor should be accounted for - pot_channel_infiltation = pot_channel_infiltation * channel_rat # is in m/s + pot_channel_infiltation = pot_channel_infiltation * self.channel_rat # is in m/s # Groundwater LDD, gradient and flow flux gw_ldd = lfr.d8_flow_direction(gw_height) del_h_gw = gw_height - lfr.downstream(gw_ldd, gw_height) gw_grad = (del_h_gw) / self.resolution - gw_flow = self.Ks * gw_grad * timestep * (gw_height - self.imperm_lay_height) * self.resolution # Groundwater velocity in m2/s + gw_flow = self.Ks * gw_grad * self.timestep * (gw_height - self.imperm_lay_height) * self.resolution # Groundwater velocity in m2/s # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. gw_flow = lfr.where(gw_flow * dt > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/dt, gw_flow) @@ -205,30 +230,7 @@ def dynamic_model(self, configuration, report): sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters for j in range(dt): - # The groundwater is adjusted by the fluxes - # channel_infiltation = lfr.where(height > pot_channel_infiltation, pot_channel_infiltation, height) - gw_s = gw_s + gw_flux #+ channel_infiltation*infil_to_gw_s - - # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. - seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) - - # Discharge is affected by the surfacewater fluxes, and seepage is added - height = height + ((sw_flux + seepage)/channel_area) #- channel_infiltation - - discharge = lfr.pow(height, c) / coefficient - - # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. - discharge = lfr.where(discharge < 1E-20, 1E-20, discharge) - - # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. - discharge = lfr.kinematic_wave(self.ldd, discharge, inflow,\ - alpha, beta, timestep,\ - channel_length,) - - height = lfr.pow(coefficient*discharge, 0.6) - - # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage - gw_s = gw_s - (seepage / self.porosity) + height, discharge, gw_s = self.update_and_route(gw_s, gw_flux, height, sw_flux) # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() From f308b67094ea9ccde984370453aba83306aa7b1b Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Thu, 11 May 2023 12:02:12 +0200 Subject: [PATCH 02/22] fix --- config/config.ini | 2 +- model/HBM.py | 6 +++--- model/reporting.py | 25 ++++++++++++++++++++----- 3 files changed, 24 insertions(+), 9 deletions(-) diff --git a/config/config.ini b/config/config.ini index 4f97352..d6a27d4 100644 --- a/config/config.ini +++ b/config/config.ini @@ -27,7 +27,7 @@ demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 [modelSettings] # Date = y, m, d, h, m, s startDate = 2023, 4, 6, 12, 30, 0 -endDate = 2023, 4, 7, 13, 30, 0 +endDate = 2023, 4, 6, 22, 30, 0 iterationsBeforeReport = 60 timestep = 1 diff --git a/model/HBM.py b/model/HBM.py index 37c1250..dbbd2b0 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -82,7 +82,7 @@ def __init__(self, configuration): self.slope = lfr.slope(self.dem, self.resolution) self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(self.slope, 0.05, 0.00001) self.channel_width = 1 - self.coefficient = self.mannings / (self.slope_sqrd * self.width) + self.coefficient = self.mannings / (self.slope_sqrd * self.channel_width) self.imperm_below_dem = float(configuration.modelSettings['impermeableLayerBelowDEM']) self.imperm_lay_height = self.dem - self.imperm_below_dem self.water_below_dem = float(configuration.modelSettings['waterBelowDEM']) @@ -93,7 +93,7 @@ def __init__(self, configuration): # Channel length and area self.channel_length = self.resolution * self.standard_LUE.one() - self.channel_area = self.width * self.channel_length + self.channel_area = self.channel_width * self.channel_length self.channel_rat = self.channel_area / self.cell_area self.infil_to_gw_s = self.channel_area / self.porosity # self.notBoundaryCells = generate.boundaryCell() # Currently not working @@ -155,7 +155,7 @@ def update_and_route(self, gw_s, gw_flux, height, sw_flux): self.alpha, self.beta, self.timestep,\ self.channel_length) - height = lfr.pow(self.constants.coefficient*discharge, 0.6) + height = lfr.pow(self.coefficient*discharge, 0.6) # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage gw_s = gw_s - (seepage / self.porosity) diff --git a/model/reporting.py b/model/reporting.py index 8542009..2b45b9b 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -9,10 +9,12 @@ import numpy as np from osgeo import gdal import pandas as pd +from StandardArraysLUE import StandardArraysLUE # Reporting for the HydrologicBaseModel class Report: def __init__(self, configuration): + self.standard_LUE = StandardArraysLUE(configuration) self.timestep = configuration.modelSettings['timestep'] self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings["scenario"] self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] @@ -31,7 +33,7 @@ def balance_report(self, configuration): end_date = self.string_to_datetime(configuration.modelSettings['endDate'], seperator= ", ") start_date_txt = start_date.strftime("%d/%m/%Y %H:%M") end_date_txt = end_date.strftime("%d/%m/%Y %H:%M") - print(end_date_txt) + if configuration.generalSettings['includePrecipitation'] == "True": p = pd.read_csv(configuration.generalSettings["inputDir"] + configuration.dataSettings["precipitationData"], names=["date", "p"]) @@ -41,7 +43,7 @@ def balance_report(self, configuration): else: mean_precipitation = 0 start_idx = 0 - end_idx = 1 + end_idx = (end_date - start_date).total_seconds() / 300 if configuration.generalSettings['includeEvapotranspiration'] == "True": e = pd.read_csv(configuration.generalSettings["inputDir"] + configuration.dataSettings["evapotranspirationData"], names=["date", "e"]) @@ -53,9 +55,22 @@ def balance_report(self, configuration): end_date = end_date - datetime.timedelta(minutes=1) - ini_sur_stor = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings["iniWaterHeight"]) - ini_gro_stor = self.input_dir + configuration.dataSettings['iniGroundWaterStorage'] - ini_int_stor = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings['iniInterceptionStorage']) + # Load initial files, if they cannot be loaded, assume zero (same as model does) + try: + ini_sur_stor = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings["iniWaterHeight"]) + except: + self.standard_LUE.zero() + + try: + ini_gro_stor = self.input_dir + configuration.dataSettings['iniGroundWaterStorage'] + except: + self.standard_LUE.zero() + + try: + ini_int_stor = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings['iniInterceptionStorage']) + except: + self.standard_LUE.zero() + end_sur_stor = self.tiff_to_np_sum(self.output_dir + "/{}_height_{}.tiff".format(int(configuration.modelSettings["timestep"]), end_date.strftime("%Y-%m-%d-%H%M"))) end_gro_stor = self.output_dir + "/{}_gw_s_{}.tiff".format(int(configuration.modelSettings["timestep"]), end_date.strftime("%Y-%m-%d-%H%M")) end_int_stor = self.tiff_to_np_sum(self.output_dir + "/{}_int_s_{}.tiff".format(int(configuration.modelSettings["timestep"]), end_date.strftime("%Y-%m-%d-%H%M"))) From b1938558bc3937a197a98562acb978ea1a4841f7 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Thu, 11 May 2023 12:21:26 +0200 Subject: [PATCH 03/22] Fix --- config/config.ini | 3 ++- model/HBM.py | 15 ++++++++------- model/utilityFunctionsHBM.py | 13 ++++++++++++- 3 files changed, 22 insertions(+), 9 deletions(-) diff --git a/config/config.ini b/config/config.ini index d6a27d4..888c4a1 100644 --- a/config/config.ini +++ b/config/config.ini @@ -27,9 +27,10 @@ demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 [modelSettings] # Date = y, m, d, h, m, s startDate = 2023, 4, 6, 12, 30, 0 -endDate = 2023, 4, 6, 22, 30, 0 +endDate = 2023, 4, 6, 13, 30, 0 iterationsBeforeReport = 60 timestep = 1 +standard_minimal_value = 1E-20 waterBelowDEM = 0.0 impermeableLayerBelowDEM= 2.00 diff --git a/model/HBM.py b/model/HBM.py index dbbd2b0..57ab1d7 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -79,6 +79,7 @@ def __init__(self, configuration): ## GENERAL ## self.resolution = float(configuration.modelSettings['resolution']) self.cell_area = self.resolution * self.resolution + self.stdmin = float(configuration.modelSettings['standard_minimal_value']) self.slope = lfr.slope(self.dem, self.resolution) self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(self.slope, 0.05, 0.00001) self.channel_width = 1 @@ -105,7 +106,7 @@ def __init__(self, configuration): self.c = 5/3 # Static, really small value because inflow = 0 is not accepted - self.inflow = self.standard_LUE.one()*1E-20 + self.inflow = self.standard_LUE.one()*self.stdmin # Load initial groundWaterStorage, if no raster is supplied, use the waterBelowDEM in combination with DEM to create a initialGroundWaterStorage layer. @@ -148,7 +149,7 @@ def update_and_route(self, gw_s, gw_flux, height, sw_flux): discharge = lfr.pow(height, self.c) / self.coefficient # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. - discharge = lfr.where(discharge < 1E-20, 1E-20, discharge) + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. discharge = lfr.kinematic_wave(self.ldd, discharge, self.inflow,\ @@ -165,10 +166,8 @@ def update_and_route(self, gw_s, gw_flux, height, sw_flux): @lfr.runtime_scope def dynamic_model(self, configuration, report): dt = int(configuration.modelSettings['iterationsBeforeReport']) - sD = list(map(int, configuration.modelSettings['startDate'].split(", "))) - eD = list(map(int, configuration.modelSettings['endDate'].split(", "))) - start_date = datetime.datetime(sD[0], sD[1], sD[2], sD[3], sD[4], sD[5]) - end_date = datetime.datetime(eD[0], eD[1], eD[2], eD[3], eD[4], eD[5]) + start_date = utilityFunctions.string_to_datetime(configuration.modelSettings['startDate'], ", ") + end_date = utilityFunctions.string_to_datetime(configuration.modelSettings['endDate'], ", ") dT = int((end_date - start_date).total_seconds() / dt) # Loading initial conditions @@ -223,13 +222,15 @@ def dynamic_model(self, configuration, report): # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. gw_flow = lfr.where(gw_flow * dt > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/dt, gw_flow) - gw_flow = lfr.where(gw_s < self.min_gw_s, 1E-20, gw_flow) + gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) # Add all vertical processes for the surfacewater and all processes groundwater gw_flux = ((direct_infiltration - evapotranspiration_soil)/self.porosity) + lfr.upstream(gw_ldd, gw_flow) - gw_flow # Is now in cubic meters sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters for j in range(dt): + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. height, discharge, gw_s = self.update_and_route(gw_s, gw_flux, height, sw_flux) # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) diff --git a/model/utilityFunctionsHBM.py b/model/utilityFunctionsHBM.py index 2b300f1..fbef096 100644 --- a/model/utilityFunctionsHBM.py +++ b/model/utilityFunctionsHBM.py @@ -6,6 +6,7 @@ """ import lue.framework as lfr +import datetime class utilityFunctions: def calculate_sqrd_slope(slope, max: float, min: float): @@ -13,4 +14,14 @@ def calculate_sqrd_slope(slope, max: float, min: float): slope_sqrd = lfr.where(slope_sqrd < min, min, slope_sqrd) slope_sqrd = lfr.where(slope_sqrd > max, max, slope_sqrd) - return slope_sqrd \ No newline at end of file + return slope_sqrd + + def string_to_datetime(date_string: str, seperator: str): + date_int_list = list(map(int, date_string.split(seperator))) + datetime_date = datetime.datetime(date_int_list[0], + date_int_list[1], + date_int_list[2], + date_int_list[3], + date_int_list[4], + date_int_list[5]) + return datetime_date \ No newline at end of file From 353045bb1d5faaf11a3225cd05fc44f2de780f0f Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Mon, 15 May 2023 11:00:26 +0200 Subject: [PATCH 04/22] Added routing methods --- config/config.ini | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/config/config.ini b/config/config.ini index 888c4a1..727f8da 100644 --- a/config/config.ini +++ b/config/config.ini @@ -27,20 +27,23 @@ demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 [modelSettings] # Date = y, m, d, h, m, s startDate = 2023, 4, 6, 12, 30, 0 -endDate = 2023, 4, 6, 13, 30, 0 +endDate = 2023, 4, 6, 13, 0, 0 iterationsBeforeReport = 60 timestep = 1 standard_minimal_value = 1E-20 +routing = both waterBelowDEM = 0.0 impermeableLayerBelowDEM= 2.00 groundWaterBase = 23.25 porosity = 0.35 +channel_width = 1 arrayExtent = 1000 partitionExtent = 1000 resolution = 5 -validCellsPercentage = 35.62 +validCellsPercentageD = 24.7 +validCellsPercentageG = 25.04 [dataSettings] iniGroundWaterStorage = From 71633b156375d749b062385f9e6390a18d8dc961 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Mon, 15 May 2023 11:01:35 +0200 Subject: [PATCH 05/22] Added reporting submodule and routing methods --- model/reporting.py | 64 ++++++++++++++++++++++++++-------------------- 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/model/reporting.py b/model/reporting.py index 2b45b9b..da88306 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -18,6 +18,12 @@ def __init__(self, configuration): self.timestep = configuration.modelSettings['timestep'] self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings["scenario"] self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] + + def initial(self, variables): + for variable, data in variables.items(): + lfr.to_gdal(data, self.output_dir + "/{}_{}.tiff".format(self.timestep, + variable + )) def dynamic(self, date, variables): dateTime = date.strftime("%Y-%m-%d-%H%M") @@ -26,7 +32,6 @@ def dynamic(self, date, variables): variable, dateTime )) - return 0 def balance_report(self, configuration): start_date = self.string_to_datetime(configuration.modelSettings['startDate'], seperator= ", ") @@ -57,40 +62,45 @@ def balance_report(self, configuration): # Load initial files, if they cannot be loaded, assume zero (same as model does) try: - ini_sur_stor = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings["iniWaterHeight"]) + ini_sur_height = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings["iniWaterHeight"]) except: - self.standard_LUE.zero() + ini_sur_height = self.tiff_to_np_sum(self.output_dir + "/{}_ini_sur_h.tiff".format(int(configuration.modelSettings["timestep"]))) + ini_sur_height = self.tiff_to_np_sum(self.output_dir + "/{}_height_{}.tiff".format(int(configuration.modelSettings["timestep"]), start_date.strftime("%Y-%m-%d-%H%M"))) try: ini_gro_stor = self.input_dir + configuration.dataSettings['iniGroundWaterStorage'] + ini_gro_data = gdal.Open(ini_gro_stor) + img = ini_gro_data.GetRasterBand(1) except: - self.standard_LUE.zero() + ini_gro_stor = (self.output_dir + "/{}_ini_gw_s.tiff".format(int(configuration.modelSettings["timestep"]))) try: ini_int_stor = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings['iniInterceptionStorage']) except: - self.standard_LUE.zero() + ini_int_stor = self.tiff_to_np_sum(self.output_dir + "/{}_ini_int_s.tiff".format(int(configuration.modelSettings["timestep"]))) + print(ini_int_stor) - end_sur_stor = self.tiff_to_np_sum(self.output_dir + "/{}_height_{}.tiff".format(int(configuration.modelSettings["timestep"]), end_date.strftime("%Y-%m-%d-%H%M"))) + end_sur_height = self.tiff_to_np_sum(self.output_dir + "/{}_height_{}.tiff".format(int(configuration.modelSettings["timestep"]), end_date.strftime("%Y-%m-%d-%H%M"))) end_gro_stor = self.output_dir + "/{}_gw_s_{}.tiff".format(int(configuration.modelSettings["timestep"]), end_date.strftime("%Y-%m-%d-%H%M")) end_int_stor = self.tiff_to_np_sum(self.output_dir + "/{}_int_s_{}.tiff".format(int(configuration.modelSettings["timestep"]), end_date.strftime("%Y-%m-%d-%H%M"))) resolution = (int(configuration.modelSettings["resolution"])) + channel_area= resolution * int(configuration.modelSettings['channel_width']) cell_area = resolution ** 2 - del_sur_stor = (end_sur_stor - ini_sur_stor) * resolution + del_sur_stor = (end_sur_height - ini_sur_height) * channel_area del_gro_stor = self.tiff_to_np_sum_difference(end_gro_stor, ini_gro_stor) * float(configuration.modelSettings["porosity"]) del_int_stor = (end_int_stor - ini_int_stor) net_balance = del_sur_stor + del_int_stor + del_gro_stor - precipitation = (((end_idx - start_idx) / 12) * mean_precipitation) / 1000 * cell_area * (int(configuration.modelSettings["arrayExtent"]) ** 2) * (float(configuration.modelSettings["validCellsPercentage"]))/100 - evapotranspiration = (((end_idx - start_idx) / 12) * mean_evapotranspiration) / 1000 * cell_area * (int(configuration.modelSettings["arrayExtent"]) ** 2) * (float(configuration.modelSettings["validCellsPercentage"]))/100 + precipitation = (((end_idx - start_idx) / 12) * mean_precipitation) / 1000 * cell_area * (int(configuration.modelSettings["arrayExtent"]) ** 2) * (float(configuration.modelSettings["validCellsPercentageG"]))/100 + evapotranspiration = (((end_idx - start_idx) / 12) * mean_evapotranspiration) / 1000 * cell_area * (int(configuration.modelSettings["arrayExtent"]) ** 2) * (float(configuration.modelSettings["validCellsPercentageG"]))/100 atmospheric_balance = precipitation - evapotranspiration print("total precipitation: ", precipitation, "m3") print("total evapotranspiration: ", evapotranspiration, "m3 \n") - print("iniTotalSurfaceHeight: ", ini_sur_stor * resolution, "m3") - print("endTotalSurfaceHeight: ", end_sur_stor * resolution, "m3") + print("iniTotalSurfaceHeight: ", ini_sur_height * channel_area, "m3") + print("endTotalSurfaceHeight: ", end_sur_height * channel_area, "m3") print("delta Surface Storage: ", del_sur_stor, "m3 \n") print("delta GroundWater Storage: ", del_gro_stor, "m3 \n") @@ -113,13 +123,10 @@ def balance_report(self, configuration): return 0 def tiff_to_np_sum(self, file): - try: - data = gdal.Open(file) - img = data.GetRasterBand(1) - raster = img.ReadAsArray() - np_array_sum = np.nansum(raster) - except: - np_array_sum = 0 + data = gdal.Open(file) + img = data.GetRasterBand(1) + raster = img.ReadAsArray() + np_array_sum = np.nansum(raster) return np_array_sum def string_to_datetime(self, date_string: str, seperator: str): @@ -134,14 +141,15 @@ def string_to_datetime(self, date_string: str, seperator: str): def tiff_to_np_sum_difference(self, file1, file2): - try: - data1 = gdal.Open(file1) - img1 = data1.GetRasterBand(1) - raster1 = img1.ReadAsArray() - data2 = gdal.Open(file2) - img2 = data2.GetRasterBand(1) - raster2 = img2.ReadAsArray() - np_array_sum = np.nansum(raster1 - raster2) - except: - np_array_sum = 0 + data1 = gdal.Open(file1) + img1 = data1.GetRasterBand(1) + raster1 = img1.ReadAsArray() + + data2 = gdal.Open(file2) + img2 = data2.GetRasterBand(1) + raster2 = img2.ReadAsArray() + + np_array_sum = np.nansum(raster1 - raster2) + + return np_array_sum \ No newline at end of file From 9365d3edc4d47c433d1be7804e142f068fd92c6d Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Mon, 15 May 2023 11:02:30 +0200 Subject: [PATCH 06/22] Added report submodule and routing methods --- model/HBM.py | 131 +++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 106 insertions(+), 25 deletions(-) diff --git a/model/HBM.py b/model/HBM.py index 57ab1d7..cc0c974 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -41,12 +41,13 @@ ) class mainModel: - def __init__(self, configuration): + def __init__(self, configuration, report): print("Initializing the program...") # Initialize submodules self.standard_LUE = StandardArraysLUE(configuration) self.retrieve_data = RetrieveData(configuration) self.calculate_flux = CalculateFlux(configuration) + self.report = report # Set directories self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] @@ -82,8 +83,8 @@ def __init__(self, configuration): self.stdmin = float(configuration.modelSettings['standard_minimal_value']) self.slope = lfr.slope(self.dem, self.resolution) self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(self.slope, 0.05, 0.00001) - self.channel_width = 1 - self.coefficient = self.mannings / (self.slope_sqrd * self.channel_width) + self.channel_width = int(configuration.modelSettings['channel_width']) + self.coefficient = self.mannings / (0.008 * self.channel_width) self.imperm_below_dem = float(configuration.modelSettings['impermeableLayerBelowDEM']) self.imperm_lay_height = self.dem - self.imperm_below_dem self.water_below_dem = float(configuration.modelSettings['waterBelowDEM']) @@ -108,7 +109,6 @@ def __init__(self, configuration): # Static, really small value because inflow = 0 is not accepted self.inflow = self.standard_LUE.one()*self.stdmin - # Load initial groundWaterStorage, if no raster is supplied, use the waterBelowDEM in combination with DEM to create a initialGroundWaterStorage layer. try: self.ini_gw_s = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniGroundWaterStorage'], partition_shape) @@ -121,23 +121,28 @@ def __init__(self, configuration): # Load initial discharge, if no raster is supplied, set to zero. try: - self.ini_water_h = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniWaterHeight'], partition_shape) + self.ini_sur_h = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniWaterHeight'], partition_shape) except: print("Did not find a initial discharge file, looked at: {}".format(configuration.dataSettings['iniWaterHeight'])) - self.ini_water_h = self.standard_LUE.zero() + self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) # Initial InterceptionStorage and groundWaterStorage try: self.ini_int_s = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniInterceptionStorage'], partition_shape) except: print("Did not find a initial interception storage file, looked at: {}".format(configuration.dataSettings['iniInterceptionStorage'])) - self.ini_int_s = self.standard_LUE.zero() + self.ini_int_s = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + # REPORTING INITIAL CONDITIONS + variables = {"ini_gw_s": self.ini_gw_s, "ini_int_s": self.ini_int_s, "ini_sur_h": self.ini_sur_h, + } + self.report.initial(variables) + print("\n") - def update_and_route(self, gw_s, gw_flux, height, sw_flux): + def update_and_route_both(self, gw_s, gw_flux, height, sw_flux): # The groundwater is adjusted by the fluxes - # channel_infiltation = lfr.where(height > pot_channel_infiltation, pot_channel_infiltation, height) + #channel_infiltation = lfr.where(height > pot_channel_infiltation, pot_channel_infiltation, height) gw_s = gw_s + gw_flux #+ channel_infiltation*infil_to_gw_s # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. @@ -163,18 +168,51 @@ def update_and_route(self, gw_s, gw_flux, height, sw_flux): return height, discharge, gw_s + def update_and_route_surface(self, height, sw_flux): + # Discharge is affected by the surfacewater fluxes, and seepage is added + height = height + ((sw_flux)/self.channel_area) #- channel_infiltation + + discharge = lfr.pow(height, self.c) / self.coefficient + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, self.inflow,\ + self.alpha, self.beta, self.timestep,\ + self.channel_length) + + height = lfr.pow(self.coefficient*discharge, 0.6) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + return height, discharge + + def update_and_route_subsurface(self, gw_s, gw_flux): + # The groundwater is adjusted by the fluxes + gw_s = gw_s + gw_flux + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) + + gw_s = gw_s - (seepage / self.porosity) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + + return gw_s, seepage + @lfr.runtime_scope - def dynamic_model(self, configuration, report): + def dynamic_model(self, configuration): dt = int(configuration.modelSettings['iterationsBeforeReport']) start_date = utilityFunctions.string_to_datetime(configuration.modelSettings['startDate'], ", ") end_date = utilityFunctions.string_to_datetime(configuration.modelSettings['endDate'], ", ") dT = int((end_date - start_date).total_seconds() / dt) # Loading initial conditions - height = self.ini_water_h + height = self.ini_sur_h gw_s = self.ini_gw_s int_s = self.ini_int_s gw_height = self.imperm_lay_height + gw_s/self.cell_area + discharge = lfr.pow(height, self.c) / self.coefficient # Open file to write maximum discharge values to for post simulation validation. with open(self.output_dir + "/maximumDischarge.csv", "w", newline="") as f: @@ -223,29 +261,72 @@ def dynamic_model(self, configuration, report): # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. gw_flow = lfr.where(gw_flow * dt > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/dt, gw_flow) gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) + gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), 0.000001) + out_flow = gw_flow + in_flow = lfr.upstream(gw_ldd, gw_flow) # Add all vertical processes for the surfacewater and all processes groundwater - gw_flux = ((direct_infiltration - evapotranspiration_soil)/self.porosity) + lfr.upstream(gw_ldd, gw_flow) - gw_flow # Is now in cubic meters + gw_flux = in_flow - out_flow sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters - for j in range(dt): - # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly - # Then route the water lateraly and repeat this for the amount of iterations required. - height, discharge, gw_s = self.update_and_route(gw_s, gw_flux, height, sw_flux) - - # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) - outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() - print("outflow: ", outflow) + + ############## ROUTING AND UPDATE PER TIMESTEP ############## + + if configuration.modelSettings['routing'] == 'both': + for j in range(dt): + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + height, discharge, gw_s = self.update_and_route_both(gw_s, gw_flux, height, sw_flux) + + # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) + outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() + print("outflow: ", outflow) + + # Write value to csv for later validation + writer.writerow([i*60 + j, outflow]) + + elif configuration.modelSettings['routing'] == 'surface': + for j in range(dt): + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + height, discharge = self.update_and_route_surface(height, sw_flux) + + # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) + outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() + total_volume = lfr.minimum(lfr.zonal_sum(height, self.dem>0)).get() + print("outflow: ", outflow, "total: ", total_volume) + + # Write value to csv for later validation + writer.writerow([i*60 + j, outflow]) - # Write value to csv for later validation - writer.writerow([i*60 + j, outflow]) + elif configuration.modelSettings['routing'] == 'subsurface': + for j in range(dt): + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + gw_s, seepage = self.update_and_route_subsurface(gw_s, gw_flux) + + outflow = lfr.maximum(lfr.zonal_sum(seepage, self.dem>0)).get() + total_volume = lfr.minimum(lfr.zonal_sum(gw_s, self.dem>0)).get() + print("outflow: ", outflow, "total_volume: ", total_volume) + + # Write value to csv for later validation + writer.writerow([i*60 + j, outflow]) + + else: + raise("Unknown 'routing' mode, change configuration setting to one of the following: both, surface, subsurface") + break + + # Adjust the GW Table for the LDD creation of the next timestep. gw_height = self.imperm_lay_height + gw_s/self.cell_area + + + # Save / Report data print(f"Done: {i+1}/{dT}") - variables = {"discharge": discharge, "int_s": int_s, "height": height, "gw_s": gw_s, + variables = {"discharge": discharge, "int_s": int_s, "height": height, "gw_s": gw_s, "gw_flow": gw_flow, "gw_flux": gw_flux } report.dynamic(date, variables) return 0 @@ -277,8 +358,8 @@ def dynamic_model(self, configuration, report): # Run the main model configuration = Configuration("F:/Projecten intern (2023)/Stage Steven Hosper/Model/v1/config/config.ini") report = Report(configuration) - main = mainModel(configuration) - main.dynamic_model(configuration, report) + main = mainModel(configuration, report) + main.dynamic_model(configuration) report.balance_report(configuration) # Process the results into a gif From 55f7cdeb112cbf0408e87601c7b1c3afeed56d96 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Fri, 19 May 2023 10:23:47 +0200 Subject: [PATCH 07/22] subsurface test --- config/config.ini | 6 ++--- model/HBM.py | 64 ++++++++++++++++++++++++++++++++++++---------- model/reporting.py | 2 +- 3 files changed, 54 insertions(+), 18 deletions(-) diff --git a/config/config.ini b/config/config.ini index 727f8da..7ae2b1f 100644 --- a/config/config.ini +++ b/config/config.ini @@ -27,11 +27,11 @@ demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 [modelSettings] # Date = y, m, d, h, m, s startDate = 2023, 4, 6, 12, 30, 0 -endDate = 2023, 4, 6, 13, 0, 0 +endDate = 2023, 4, 8, 12, 30, 0 iterationsBeforeReport = 60 -timestep = 1 +timestep = 15 standard_minimal_value = 1E-20 -routing = both +routing = subsurface waterBelowDEM = 0.0 impermeableLayerBelowDEM= 2.00 diff --git a/model/HBM.py b/model/HBM.py index cc0c974..ed3b2b5 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -81,7 +81,7 @@ def __init__(self, configuration, report): self.resolution = float(configuration.modelSettings['resolution']) self.cell_area = self.resolution * self.resolution self.stdmin = float(configuration.modelSettings['standard_minimal_value']) - self.slope = lfr.slope(self.dem, self.resolution) + self.slope = self.standard_LUE.one() * 0.008 self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(self.slope, 0.05, 0.00001) self.channel_width = int(configuration.modelSettings['channel_width']) self.coefficient = self.mannings / (0.008 * self.channel_width) @@ -92,6 +92,7 @@ def __init__(self, configuration, report): self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point # Refactorings value from mm/hour to m/h times the cell area. self.refactor = (self.cell_area / 1000) / 3600 + self.d = (5**(2/3))/(4**(2/3)) # Channel length and area self.channel_length = self.resolution * self.standard_LUE.one() @@ -124,7 +125,7 @@ def __init__(self, configuration, report): self.ini_sur_h = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniWaterHeight'], partition_shape) except: print("Did not find a initial discharge file, looked at: {}".format(configuration.dataSettings['iniWaterHeight'])) - self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.one() * 0.0001, self.standard_LUE.zero()) # Initial InterceptionStorage and groundWaterStorage try: @@ -151,7 +152,9 @@ def update_and_route_both(self, gw_s, gw_flux, height, sw_flux): # Discharge is affected by the surfacewater fluxes, and seepage is added height = height + ((sw_flux + seepage)/self.channel_area) #- channel_infiltation - discharge = lfr.pow(height, self.c) / self.coefficient + #discharge = lfr.pow(height, self.c) / self.coefficient + + discharge = self.height_to_dis(height) # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) @@ -161,7 +164,9 @@ def update_and_route_both(self, gw_s, gw_flux, height, sw_flux): self.alpha, self.beta, self.timestep,\ self.channel_length) - height = lfr.pow(self.coefficient*discharge, 0.6) + height = self.dis_to_height(discharge) + + #height = lfr.pow(self.coefficient*discharge, 0.6) # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage gw_s = gw_s - (seepage / self.porosity) @@ -172,7 +177,9 @@ def update_and_route_surface(self, height, sw_flux): # Discharge is affected by the surfacewater fluxes, and seepage is added height = height + ((sw_flux)/self.channel_area) #- channel_infiltation - discharge = lfr.pow(height, self.c) / self.coefficient + #discharge = lfr.pow(height, self.c) / self.coefficient + + discharge = self.height_to_dis(height) # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) @@ -182,7 +189,9 @@ def update_and_route_surface(self, height, sw_flux): self.alpha, self.beta, self.timestep,\ self.channel_length) - height = lfr.pow(self.coefficient*discharge, 0.6) + #height = lfr.pow(self.coefficient*discharge, 0.6) + + height = self.dis_to_height(discharge) # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage return height, discharge @@ -200,12 +209,44 @@ def update_and_route_subsurface(self, gw_s, gw_flux): return gw_s, seepage + def dis_to_height(self, discharge): + """Use the open channel flow formula to convert discharge to height + + The Open Channel Flow calculation for rectangular channels using the known width, slope and mannings. + + Args: + discharge (lue partitioned array): discharge array that is a result of the kinematic wave formule. + + Returns: + height (lue partitioned array): height array that can be used to add vertical and lateral in- and outflows. + + """ + height = (1/self.d)*self.mannings * discharge / self.slope_sqrd + + return height + + def height_to_dis(self, height): + """Use the open channel flow formula to convert discharge to height + + The Open Channel Flow calculation for rectangular channels using the known width, slope and mannings. + + Args: + height (lue partitioned array): height array after all in- and outfluxes have been added. + + Returns: + discharge (lue partitioned array): discharge array that can be used as an input for the kinematic wave formula. + + """ + discharge = self.d * (height * self.slope_sqrd)/self.mannings + + return discharge + @lfr.runtime_scope def dynamic_model(self, configuration): dt = int(configuration.modelSettings['iterationsBeforeReport']) start_date = utilityFunctions.string_to_datetime(configuration.modelSettings['startDate'], ", ") end_date = utilityFunctions.string_to_datetime(configuration.modelSettings['endDate'], ", ") - dT = int((end_date - start_date).total_seconds() / dt) + dT = int((end_date - start_date).total_seconds() / (dt*self.timestep)) # Loading initial conditions height = self.ini_sur_h @@ -262,11 +303,9 @@ def dynamic_model(self, configuration): gw_flow = lfr.where(gw_flow * dt > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/dt, gw_flow) gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), 0.000001) - - out_flow = gw_flow - in_flow = lfr.upstream(gw_ldd, gw_flow) + # Add all vertical processes for the surfacewater and all processes groundwater - gw_flux = in_flow - out_flow + gw_flux = lfr.upstream(gw_ldd, gw_flow) - gw_flow sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters @@ -307,7 +346,6 @@ def dynamic_model(self, configuration): outflow = lfr.maximum(lfr.zonal_sum(seepage, self.dem>0)).get() total_volume = lfr.minimum(lfr.zonal_sum(gw_s, self.dem>0)).get() - print("outflow: ", outflow, "total_volume: ", total_volume) # Write value to csv for later validation writer.writerow([i*60 + j, outflow]) @@ -322,8 +360,6 @@ def dynamic_model(self, configuration): gw_height = self.imperm_lay_height + gw_s/self.cell_area - - # Save / Report data print(f"Done: {i+1}/{dT}") variables = {"discharge": discharge, "int_s": int_s, "height": height, "gw_s": gw_s, "gw_flow": gw_flow, "gw_flux": gw_flux diff --git a/model/reporting.py b/model/reporting.py index da88306..b5f59dc 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -65,7 +65,7 @@ def balance_report(self, configuration): ini_sur_height = self.tiff_to_np_sum(self.input_dir + configuration.dataSettings["iniWaterHeight"]) except: ini_sur_height = self.tiff_to_np_sum(self.output_dir + "/{}_ini_sur_h.tiff".format(int(configuration.modelSettings["timestep"]))) - ini_sur_height = self.tiff_to_np_sum(self.output_dir + "/{}_height_{}.tiff".format(int(configuration.modelSettings["timestep"]), start_date.strftime("%Y-%m-%d-%H%M"))) + #ini_sur_height = self.tiff_to_np_sum(self.output_dir + "/{}_height_{}.tiff".format(int(configuration.modelSettings["timestep"]), start_date.strftime("%Y-%m-%d-%H%M"))) try: ini_gro_stor = self.input_dir + configuration.dataSettings['iniGroundWaterStorage'] From c8eb0a47fd8c282c695e39ec6cf3c82056e2a2fa Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Fri, 2 Jun 2023 11:12:07 +0200 Subject: [PATCH 08/22] Updates and fixes --- config/config.ini | 22 +++++---- model/CalculateFlux.py | 109 ++++++++++++++++++++++++++++------------ model/HBM.py | 110 ++++++++++++++++++++++++----------------- model/reporting.py | 5 +- 4 files changed, 158 insertions(+), 88 deletions(-) diff --git a/config/config.ini b/config/config.ini index 7ae2b1f..f4d92a4 100644 --- a/config/config.ini +++ b/config/config.ini @@ -10,8 +10,8 @@ scenario = De Hupsel5 makeGIF = False # Include Processes -includePrecipitation = False -includeEvapotranspiration = False +includePrecipitation = True +includeEvapotranspiration = True includeInfiltration = True includeInterception = True includePercolation = False @@ -28,10 +28,11 @@ demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 # Date = y, m, d, h, m, s startDate = 2023, 4, 6, 12, 30, 0 endDate = 2023, 4, 8, 12, 30, 0 -iterationsBeforeReport = 60 -timestep = 15 +iterationsBeforeData = 300 +iterationsBeforeReport = 900 +timestep = 1 standard_minimal_value = 1E-20 -routing = subsurface +routing = both waterBelowDEM = 0.0 impermeableLayerBelowDEM= 2.00 @@ -46,17 +47,18 @@ validCellsPercentageD = 24.7 validCellsPercentageG = 25.04 [dataSettings] -iniGroundWaterStorage = -iniWaterHeight = +iniGroundWaterStorage = /1_gw_s_2022-04-07-1215.tiff +iniWaterHeight = /1_height_2022-04-07-1215.tiff +iniDischarge = /1_discharge_2022-04-07-1215.tiff iniInterceptionStorage = -dem = /v2/dem_f64_shaped.tiff +dem = /dem_f64_pcr_shaped.tiff ldd = /ldd_f64_pcr_shaped.tiff soilMap = /bodem.tiff landUseMap = /landgebruik.tiff soilData = /soil_conversion.csv landUseData = /landuse_conversion_v2.csv -precipitationData = -evapotranspirationData = +precipitationData = /precipitation.csv +evapotranspirationData = /evapotranspiration.csv [reportSettings] variables = discharge, seepage, groundWaterHeight, Qgw, Sgw, swFlux, gwFlux, infiltration, evapotranspirationSoil diff --git a/model/CalculateFlux.py b/model/CalculateFlux.py index 5cbd834..a641c6e 100644 --- a/model/CalculateFlux.py +++ b/model/CalculateFlux.py @@ -16,14 +16,15 @@ def __init__(self, configuration): 2) Set the processes that should be included. """ self.std_arr_lue = StandardArraysLUE(configuration) - - self.iterations = int(configuration.modelSettings['iterationsBeforeReport']) + self.iterationsData = int(configuration.modelSettings['iterationsBeforeData']) + self.timestep = int(configuration.modelSettings['timestep']) self.include_precipitation = configuration.generalSettings['includePrecipitation'] self.include_evapotranspiration = configuration.generalSettings['includeEvapotranspiration'] self.include_infiltration = configuration.generalSettings['includeInfiltration'] self.include_interception = configuration.generalSettings['includeInterception'] self.include_percolation = configuration.generalSettings['includePercolation'] self.ca = int(configuration.modelSettings['resolution'])**2 + self.T = 0 def interception(self, int_stor, int_stor_max, pre, rev, th_f): @@ -35,16 +36,16 @@ def interception(self, int_stor, int_stor_max, pre, rev, th_f): at the surface, because the evaporative potential is already met by the vegetation. Args: - int_stor (lpa*): Current water stored from interception in the vegetation - int_stor_max (lpa*): Maximum amount of interception stored - pre (lpa*): Precipitation rate - rev (lpa*): Reference evapotranspiration rate - thF (lpa*): Throughfall fraction + int_stor (lpa*): Current water stored from interception in the vegetation (m3) + int_stor_max (lpa*): Maximum amount of interception stored (m3) + pre (lpa*): Precipitation rate (m3/s) + rev (lpa*): Reference evapotranspiration rate (m3/s) + th_f (lpa*): Throughfall fraction (-) Returns: - int_stor (lpa*): Updated water stored from interception in the vegetation + int_stor (lpa*): Updated water stored from interception in the vegetation precipitation (lpa*): Precipitation rate that falls through the vegetation towards the soil - evapotranspiration_surface (lpa*): Evaporative potential left after interception evaporation + evapotranspiration_surface (lpa*): Evaporative potential left after interception evaporation lpa*: lue partitioned array """ @@ -52,10 +53,10 @@ def interception(self, int_stor, int_stor_max, pre, rev, th_f): interception = (self.std_arr_lue.one() - th_f) * pre # Determine if there is enough water in the canopy to meet all evaporative potential - enough_water_int = (interception + int_stor/self.iterations) > rev + enough_water_int = (interception + int_stor/self.iterationsData) > (rev * (self.std_arr_lue.one() - th_f)) # If there is enough water, update the interception storage, otherwise set at zero - int_stor = lfr.where(enough_water_int, int_stor + (interception - rev) * self.iterations, 0) + int_stor = lfr.where(enough_water_int, int_stor + (interception - (rev * (self.std_arr_lue.one() - th_f))) * self.iterationsData, 0) # If the store is exceeded, add to variable called interception storage surplus int_storSurplus = lfr.where(int_stor > int_stor_max, int_stor - int_stor_max, 0) @@ -64,21 +65,21 @@ def interception(self, int_stor, int_stor_max, pre, rev, th_f): int_stor = int_stor - int_storSurplus # Add the interception storage surplus back to the precipitation - precipitation = pre - (interception - int_storSurplus/self.iterations) + precipitation = pre - (interception - int_storSurplus/self.iterationsData) - evapotranspirationSurface = lfr.where(enough_water_int, 0, rev - (interception + int_stor/self.iterations)) + evapotranspirationSurface = lfr.where(enough_water_int, (rev * th_f), rev - (interception + int_stor/self.iterationsData) * th_f) return int_stor, precipitation, evapotranspirationSurface def evapotranspiration(self, pre, ev_s): """Determine the actual surface and soil evapotranspiration Args: - pre (lpa*): Precipitation rate - ev_s (lpa*): Evaporative potential left for the surface + pre (lpa*): Precipitation rate (m3/s) + ev_s (lpa*): Evaporative potential left for the surface (m3/s) Returns: - evapotranspiration_surface (lpa*): Actual surface evapotranspiration rate - evapotranspiration_soil (lpa*): Actual soil evaporation rate + evapotranspiration_surface (lpa*): Actual surface evapotranspiration rate (m3/s) + evapotranspiration_soil (lpa*): Actual soil evaporation rate (m3/s) lpa*: lue partitioned array """ @@ -97,13 +98,13 @@ def infiltration(self, sgw, max_sgw, Ks, prm, por, pre, ev_s): is turned into potential channel infiltration. Args: - sgw (lpa*): Groundwater storage - max_sgw (lpa*): Maximum groundwater storage - Ks (lpa*): Hydraulic conductivity - prm (lpa*): Permeability - por (lpa*): Porosity - pre (lpa*): Precipitation rate reaching the soil - ev_s (lpa*): Actual surface evapotranspiration + sgw (lpa*): Groundwater storage in m3 + max_sgw (lpa*): Maximum groundwater storage (m3) + Ks (lpa*): Hydraulic conductivity (m/s) + prm (lpa*): Permeability coefficient (-) + por (lpa*): Porosity (-) + pre (lpa*): Precipitation rate reaching the soil (m3/s) + ev_s (lpa*): Actual surface evapotranspiration (m3/s) Returns: direct_infiltration (lpa*): precipitation that directly infiltrates the soil @@ -112,13 +113,59 @@ def infiltration(self, sgw, max_sgw, Ks, prm, por, pre, ev_s): lpa*: lue partitioned array """ - pot_infiltration = Ks * prm # meters that can infiltrate the soil - pot_infiltration = lfr.where(((max_sgw - sgw)/self.ca)*por < pot_infiltration, \ - ((max_sgw - sgw)/self.ca)*por, pot_infiltration) * self.ca # The amount that can infiltrate because of capacity times the area + pot_infiltration = Ks * prm * self.ca # m3 that can infiltrate the soil + pot_infiltration = lfr.where((max_sgw - sgw)*por < pot_infiltration, \ + (max_sgw - sgw)*por, pot_infiltration) # The amount that can infiltrate because of capacity times the area pot_infiltration = lfr.where(pot_infiltration < 0, 0, pot_infiltration) enough_water_inf = pre - ev_s > pot_infiltration # If there is more water on the surface available than can infiltrate direct_infiltration = lfr.where(enough_water_inf, pot_infiltration, pre - ev_s) # Either the pot_infiltration will fully infiltrate - pot_infil_channel = lfr.where(enough_water_inf, 0, - pot_infiltration - direct_infiltration) / self.ca - return direct_infiltration, pot_infil_channel # or the available water at the surface will. - \ No newline at end of file + return direct_infiltration # or the available water at the surface will. + + def adjusted_Horton(self, sgw, max_sgw, Ks, Ki, k, prm, por, pre, ev_s): + """Created a version of the Horton infiltration that deals with dry periods. + It returns to high infiltration rates over time with a relatively simple function. + One of the drawbacks is that by using 'time' the intensity does not matter. + So a very light rain could affect the high infiltration rate, which would not be used at all. + It should however still significantly improve the model from the function used prior. + + Args: + sgw (lpa*): Groundwater storage in m3 + max_sgw (lpa*): Maximum groundwater storage (m3) + Ks (lpa*): Saturated hydraulic conductivity (m/s) + Ki (lpa*): Intial hydraulic conductivity (m/s) + k (lpa*): Decay factor (-) + prm (lpa*): Permeability coefficient (-) + por (lpa*): Porosity (-) + pre (lpa*): Precipitation rate reaching the soil (m3/s) + ev_s(lpa*): Actual surface evapotranspiration (m3/s) + + Returns: + infiltration (lpa*): the infiltration rate based on horton (m3/s) + """ + + wet = pre > ev_s + + if wet: + # Horton infiltration + pot_infiltration_hort = Ks + (Ki - Ks)*2.71828^(-k*self.T) + pot_infiltration = pot_infiltration_hort * prm * self.ca # m3 that can infiltrate the soil + pot_infiltration = lfr.where((max_sgw -sgw)*por / (self.timestep * self.iterationsData) > pot_infiltration, + pot_infiltration, + (max_sgw -sgw)*por / (self.timestep * self.iterationsData)) + + + # Add time to rainfall event + self.T = self.T + (self.iterationsData * self.timestep / 3600) + + # Check if the amount of precipitation available for infiltration is larger than the potential infiltration + enough_rain = pre - ev_s > pot_infiltration + + # Calculate actual infiltration + infiltration = lfr.where(enough_rain, pot_infiltration, pre - ev_s) + + else: + infiltration = self.std_arr_lue.zero() + self.T = self.T - (self.iterationsData * self.timestep / 3600) + self.T = lfr.where(self.T < 0, 0, self.T) + + return infiltration \ No newline at end of file diff --git a/model/HBM.py b/model/HBM.py index ed3b2b5..cf6bfc5 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -61,6 +61,8 @@ def __init__(self, configuration, report): self.dem = lfr.where(self.dem < 0.1, 35, self.dem) land_use = lfr.from_gdal(self.input_dir + configuration.dataSettings['landUseMap'], partition_shape) # Land-use, example: road soil_type = lfr.from_gdal(self.input_dir + configuration.dataSettings['soilMap'], partition_shape) # example: sand or clay + self.resolution = float(configuration.modelSettings['resolution']) + self.cell_area = self.resolution * self.resolution # Retrieve the soil properties self.Ks, self.porosity, self.wilting_point = self.retrieve_data.soil_csv( @@ -69,7 +71,8 @@ def __init__(self, configuration, report): # Retrieve the land-use properties self.mannings, self.permeability, self.max_int_s, self.throughfall_frac = \ self.retrieve_data.land_characteristics_csv( - configuration.generalSettings['inputDir'] + configuration.dataSettings['landUseData'], land_use) # land-use characteristics + configuration.generalSettings['inputDir'] + configuration.dataSettings['landUseData'], land_use, + self.cell_area) # land-use characteristics self.gw_base = float(configuration.modelSettings['groundWaterBase']) * self.standard_LUE.one() # self.ldd = lfr.d8_flow_direction(self.dem) @@ -78,8 +81,6 @@ def __init__(self, configuration, report): #### SETTING CONSTANTS #### ## GENERAL ## - self.resolution = float(configuration.modelSettings['resolution']) - self.cell_area = self.resolution * self.resolution self.stdmin = float(configuration.modelSettings['standard_minimal_value']) self.slope = self.standard_LUE.one() * 0.008 self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(self.slope, 0.05, 0.00001) @@ -90,7 +91,7 @@ def __init__(self, configuration, report): self.water_below_dem = float(configuration.modelSettings['waterBelowDEM']) self.max_gw_s = self.imperm_below_dem * self.cell_area # Full storage of porosity self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point - # Refactorings value from mm/hour to m/h times the cell area. + # Refactorings value from mm/hour to m/s times the cell area. self.refactor = (self.cell_area / 1000) / 3600 self.d = (5**(2/3))/(4**(2/3)) @@ -125,7 +126,13 @@ def __init__(self, configuration, report): self.ini_sur_h = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniWaterHeight'], partition_shape) except: print("Did not find a initial discharge file, looked at: {}".format(configuration.dataSettings['iniWaterHeight'])) - self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.one() * 0.0001, self.standard_LUE.zero()) + self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + try: + self.ini_dis = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniDischarge'], partition_shape) + except: + print("Did not find a initial discharge file, looked at: {}".format(configuration.dataSettings['iniDischarge'])) + self.ini_dis = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) # Initial InterceptionStorage and groundWaterStorage try: @@ -141,45 +148,48 @@ def __init__(self, configuration, report): print("\n") - def update_and_route_both(self, gw_s, gw_flux, height, sw_flux): + def update_and_route_both(self, gw_s, gw_flux, discharge, sw_flux): # The groundwater is adjusted by the fluxes #channel_infiltation = lfr.where(height > pot_channel_infiltation, pot_channel_infiltation, height) gw_s = gw_s + gw_flux #+ channel_infiltation*infil_to_gw_s # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) + # Discharge is affected by the surfacewater fluxes, and seepage is added - height = height + ((sw_flux + seepage)/self.channel_area) #- channel_infiltation + #height = height + ((sw_flux + seepage)/self.channel_area) #- channel_infiltation #discharge = lfr.pow(height, self.c) / self.coefficient - discharge = self.height_to_dis(height) + #discharge = self.height_to_dis(height) # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + inflow = (sw_flux + seepage)/self.channel_area + inflow = lfr.where(inflow < self.stdmin, self.stdmin, inflow) # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. - discharge = lfr.kinematic_wave(self.ldd, discharge, self.inflow,\ + discharge = lfr.kinematic_wave(self.ldd, discharge, inflow,\ self.alpha, self.beta, self.timestep,\ self.channel_length) - height = self.dis_to_height(discharge) + #height = self.dis_to_height(discharge) - #height = lfr.pow(self.coefficient*discharge, 0.6) + height = lfr.pow(self.coefficient*discharge, 0.6) # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage - gw_s = gw_s - (seepage / self.porosity) + gw_s = gw_s - (seepage / self.porosity) - return height, discharge, gw_s + return height, discharge, gw_s, seepage def update_and_route_surface(self, height, sw_flux): # Discharge is affected by the surfacewater fluxes, and seepage is added height = height + ((sw_flux)/self.channel_area) #- channel_infiltation - #discharge = lfr.pow(height, self.c) / self.coefficient + discharge = lfr.pow(height, self.c) / self.coefficient - discharge = self.height_to_dis(height) + #discharge = self.height_to_dis(height) # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) @@ -189,9 +199,9 @@ def update_and_route_surface(self, height, sw_flux): self.alpha, self.beta, self.timestep,\ self.channel_length) - #height = lfr.pow(self.coefficient*discharge, 0.6) + height = lfr.pow(self.coefficient*discharge, 0.6) - height = self.dis_to_height(discharge) + #height = self.dis_to_height(discharge) # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage return height, discharge @@ -243,17 +253,18 @@ def height_to_dis(self, height): @lfr.runtime_scope def dynamic_model(self, configuration): - dt = int(configuration.modelSettings['iterationsBeforeReport']) + dt_data = int(configuration.modelSettings['iterationsBeforeData']) + dt_report = int(configuration.modelSettings['iterationsBeforeReport']) start_date = utilityFunctions.string_to_datetime(configuration.modelSettings['startDate'], ", ") end_date = utilityFunctions.string_to_datetime(configuration.modelSettings['endDate'], ", ") - dT = int((end_date - start_date).total_seconds() / (dt*self.timestep)) + dT = int((end_date - start_date).total_seconds() / (dt_data*self.timestep)) # Loading initial conditions height = self.ini_sur_h gw_s = self.ini_gw_s int_s = self.ini_int_s gw_height = self.imperm_lay_height + gw_s/self.cell_area - discharge = lfr.pow(height, self.c) / self.coefficient + discharge = self.ini_dis # Open file to write maximum discharge values to for post simulation validation. with open(self.output_dir + "/maximumDischarge.csv", "w", newline="") as f: @@ -262,16 +273,16 @@ def dynamic_model(self, configuration): # Start model for dT large periods for i in range(dT): # Time in minutes is the small iteration multiplied with the timestep (both in seconds) divided by 60 seconds. - date = start_date + datetime.timedelta(seconds = i * (dt*self.timestep)) + date = start_date + datetime.timedelta(seconds = i * (dt_data*self.timestep)) # Load flux and storage values precipitation = self.retrieve_data.csv_timeseries_to_flux(configuration.generalSettings['inputDir'] + configuration.dataSettings['precipitationData'], - self.refactor, date) # m/s + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(configuration.generalSettings['inputDir'] + configuration.dataSettings['evapotranspirationData'], - self.refactor, date) # m/s + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) int_s, precipitation, evapotranspiration_surface = self.calculate_flux.interception(int_s, self.max_int_s, @@ -282,16 +293,13 @@ def dynamic_model(self, configuration): evapotranspiration_surface, evapotranspiration_soil = self.calculate_flux.evapotranspiration(precipitation, evapotranspiration_surface) - direct_infiltration, pot_channel_infiltation = self.calculate_flux.infiltration(gw_s, - self.max_gw_s, - self.Ks, - self.permeability, - self.porosity, - precipitation, - evapotranspiration_surface) - - # The infiltration happens only in the region that is used by the channel and therefore this factor should be accounted for - pot_channel_infiltation = pot_channel_infiltation * self.channel_rat # is in m/s + direct_infiltration = self.calculate_flux.infiltration(gw_s, + self.max_gw_s, + self.Ks, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) # Groundwater LDD, gradient and flow flux gw_ldd = lfr.d8_flow_direction(gw_height) @@ -300,32 +308,35 @@ def dynamic_model(self, configuration): gw_flow = self.Ks * gw_grad * self.timestep * (gw_height - self.imperm_lay_height) * self.resolution # Groundwater velocity in m2/s # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. - gw_flow = lfr.where(gw_flow * dt > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/dt, gw_flow) + gw_flow = lfr.where(gw_flow * dt_data > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/dt_data, gw_flow) gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) - gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), 0.000001) + gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), gw_flow) # Add all vertical processes for the surfacewater and all processes groundwater - gw_flux = lfr.upstream(gw_ldd, gw_flow) - gw_flow + gw_flux = direct_infiltration - evapotranspiration_soil + lfr.upstream(gw_ldd, gw_flow) - gw_flow sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters ############## ROUTING AND UPDATE PER TIMESTEP ############## if configuration.modelSettings['routing'] == 'both': - for j in range(dt): + for j in range(dt_data): # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly # Then route the water lateraly and repeat this for the amount of iterations required. - height, discharge, gw_s = self.update_and_route_both(gw_s, gw_flux, height, sw_flux) + height, discharge, gw_s, seepage = self.update_and_route_both(gw_s, gw_flux, discharge, sw_flux) # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() print("outflow: ", outflow) # Write value to csv for later validation - writer.writerow([i*60 + j, outflow]) + writer.writerow([i*dt_data + j, outflow]) + + variables = {"discharge": discharge, "int_s": int_s, "gw_s": gw_s, "seepage": seepage, "infiltration": direct_infiltration + } elif configuration.modelSettings['routing'] == 'surface': - for j in range(dt): + for j in range(dt_data): # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly # Then route the water lateraly and repeat this for the amount of iterations required. height, discharge = self.update_and_route_surface(height, sw_flux) @@ -338,17 +349,24 @@ def dynamic_model(self, configuration): # Write value to csv for later validation writer.writerow([i*60 + j, outflow]) + variables = {"discharge": discharge, "int_s": int_s, "height": height + } + + elif configuration.modelSettings['routing'] == 'subsurface': - for j in range(dt): + for j in range(dt_data): # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly # Then route the water lateraly and repeat this for the amount of iterations required. gw_s, seepage = self.update_and_route_subsurface(gw_s, gw_flux) - outflow = lfr.maximum(lfr.zonal_sum(seepage, self.dem>0)).get() + outflow = lfr.maximum(lfr.zonal_sum(seepage, self.dem>0)).get() / self.timestep total_volume = lfr.minimum(lfr.zonal_sum(gw_s, self.dem>0)).get() # Write value to csv for later validation writer.writerow([i*60 + j, outflow]) + + variables = {"gw_s": gw_s, "seepage": seepage, "groundwater_height": gw_height, + } else: raise("Unknown 'routing' mode, change configuration setting to one of the following: both, surface, subsurface") @@ -361,10 +379,12 @@ def dynamic_model(self, configuration): # Save / Report data + if i % (dt_report/dt_data) == 0: + print(i) + report.dynamic(date, variables) + + print(f"Done: {i+1}/{dT}") - variables = {"discharge": discharge, "int_s": int_s, "height": height, "gw_s": gw_s, "gw_flow": gw_flow, "gw_flux": gw_flux - } - report.dynamic(date, variables) return 0 diff --git a/model/reporting.py b/model/reporting.py index b5f59dc..958a9fb 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -15,7 +15,8 @@ class Report: def __init__(self, configuration): self.standard_LUE = StandardArraysLUE(configuration) - self.timestep = configuration.modelSettings['timestep'] + self.timestep = int(configuration.modelSettings['timestep']) + self.iterations = int(configuration.modelSettings['iterationsBeforeReport']) self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings["scenario"] self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] @@ -58,7 +59,7 @@ def balance_report(self, configuration): else: mean_evapotranspiration = 0 - end_date = end_date - datetime.timedelta(minutes=1) + end_date = end_date - datetime.timedelta(seconds=self.iterations * self.timestep) # Load initial files, if they cannot be loaded, assume zero (same as model does) try: From 01d35bd969f01c4c2e40fc686524fae012ad116a Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Fri, 2 Jun 2023 11:16:57 +0200 Subject: [PATCH 09/22] Implemented adjusted Horton --- model/HBM.py | 14 ++++++++++++-- model/RetrieveData.py | 15 +++++++++++---- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/model/HBM.py b/model/HBM.py index cf6bfc5..b52be98 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -65,7 +65,7 @@ def __init__(self, configuration, report): self.cell_area = self.resolution * self.resolution # Retrieve the soil properties - self.Ks, self.porosity, self.wilting_point = self.retrieve_data.soil_csv( + self.Ks, self.porosity, self.wilting_point, self.Ki, self.k = self.retrieve_data.soil_csv( configuration.generalSettings['inputDir'] + configuration.dataSettings['soilData'], soil_type) # soil characteristic # Retrieve the land-use properties @@ -293,12 +293,22 @@ def dynamic_model(self, configuration): evapotranspiration_surface, evapotranspiration_soil = self.calculate_flux.evapotranspiration(precipitation, evapotranspiration_surface) - direct_infiltration = self.calculate_flux.infiltration(gw_s, + """ direct_infiltration = self.calculate_flux.infiltration(gw_s, self.max_gw_s, self.Ks, self.permeability, self.porosity, precipitation, + evapotranspiration_surface) """ + + direct_infiltration = self.calculate_flux.adjusted_Horton(gw_s, + self.max_gw_s, + self.Ks, + self.Ki, + self.k, + self.permeability, + self.porosity, + precipitation, evapotranspiration_surface) # Groundwater LDD, gradient and flow flux diff --git a/model/RetrieveData.py b/model/RetrieveData.py index 6b577d7..86b5343 100644 --- a/model/RetrieveData.py +++ b/model/RetrieveData.py @@ -43,6 +43,7 @@ def soil_csv(self, data_file, soil_type): Ks = self.std_arr_lue.one() * 0.05 porosity = self.std_arr_lue.one() * 0.35 wilting_point = self.std_arr_lue.one() * 0.15 + k = self.std_arr_lue.one() * 1.6 # Read pandas data table data_table = pd.read_csv(data_file) @@ -50,24 +51,28 @@ def soil_csv(self, data_file, soil_type): # Split into ID and Ks value ID = data_table["ID"] Ks_value = data_table["Ks"] + Ki_value = data_table["Ki"] # When the ID of the table meets an ID within the partitioned array, assign value for count, ID in enumerate(ID): Ks = lfr.where(soil_type == ID, Ks_value[count], Ks) # To m/s + Ki = lfr.where(soil_type == ID, Ki_value[count], Ki) + Ki = Ki / 86400 Ks = Ks / 86400 - return Ks, porosity, wilting_point + return Ks, porosity, wilting_point, Ki, k - def land_characteristics_csv(self, data_file, land_use): + def land_characteristics_csv(self, data_file, land_use, cell_area): """Reads the land characteristics from a csv file Gives standard values to the entire array using the set array extent. Continues to read values from data file to the corresponding map IDs. Args: - data_file (path): the path to the csv data file containing the land characteristics information - land_use (lpa*): array containing the id that matches every cell to its \ + data_file (path): the path to the csv data file containing the land characteristics information + land_use (lpa*): array containing the id that matches every cell to its \ corresponding characteristic values. + cell_area (): cell size Returns: mannings (lpa*): Mannings Coefficient in a lue array @@ -106,6 +111,8 @@ def land_characteristics_csv(self, data_file, land_use): throughfall_fraction = lfr.where(land_use == ID, throughfall_fraction_value[count], throughfall_fraction) # crop_factor = lfr.where(land_use == ID, crop_factor_value[count], crop_factor) + interception_storage_max = interception_storage_max * cell_area + return mannings, permeability, interception_storage_max, throughfall_fraction def rounddown_datetime(self, dt): From 34a4b4d2969c64a063ca5d7fcb767951ae943184 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Mon, 12 Jun 2023 08:57:56 +0200 Subject: [PATCH 10/22] Adjustments for reinfiltration testing --- config/config.ini | 2 +- model/CalculateFlux.py | 30 ++++++++++++++++-------------- model/HBM.py | 29 +++++++++++------------------ model/RetrieveData.py | 7 ++++--- 4 files changed, 32 insertions(+), 36 deletions(-) diff --git a/config/config.ini b/config/config.ini index f4d92a4..afb6f38 100644 --- a/config/config.ini +++ b/config/config.ini @@ -47,7 +47,7 @@ validCellsPercentageD = 24.7 validCellsPercentageG = 25.04 [dataSettings] -iniGroundWaterStorage = /1_gw_s_2022-04-07-1215.tiff +iniGroundWaterStorage = /gw_test/60_gw_s_2022-04-14-1230.tiff iniWaterHeight = /1_height_2022-04-07-1215.tiff iniDischarge = /1_discharge_2022-04-07-1215.tiff iniInterceptionStorage = diff --git a/model/CalculateFlux.py b/model/CalculateFlux.py index a641c6e..dd5aa78 100644 --- a/model/CalculateFlux.py +++ b/model/CalculateFlux.py @@ -25,6 +25,7 @@ def __init__(self, configuration): self.include_percolation = configuration.generalSettings['includePercolation'] self.ca = int(configuration.modelSettings['resolution'])**2 self.T = 0 + self.k = 1.6 def interception(self, int_stor, int_stor_max, pre, rev, th_f): @@ -121,7 +122,7 @@ def infiltration(self, sgw, max_sgw, Ks, prm, por, pre, ev_s): direct_infiltration = lfr.where(enough_water_inf, pot_infiltration, pre - ev_s) # Either the pot_infiltration will fully infiltrate return direct_infiltration # or the available water at the surface will. - def adjusted_Horton(self, sgw, max_sgw, Ks, Ki, k, prm, por, pre, ev_s): + def adjusted_Horton(self, sgw, max_sgw, Ks, Ki, prm, por, pre, ev_s): """Created a version of the Horton infiltration that deals with dry periods. It returns to high infiltration rates over time with a relatively simple function. One of the drawbacks is that by using 'time' the intensity does not matter. @@ -143,17 +144,16 @@ def adjusted_Horton(self, sgw, max_sgw, Ks, Ki, k, prm, por, pre, ev_s): infiltration (lpa*): the infiltration rate based on horton (m3/s) """ - wet = pre > ev_s - + wet = lfr.maximum(pre).get() > lfr.maximum(ev_s).get() + + # Horton infiltration + pot_infiltration_hort = Ks + (Ki - Ks)*lfr.pow(self.std_arr_lue.one()*2.71828,(-1*self.k*self.T)) + pot_infiltration = pot_infiltration_hort * prm * self.ca # m3 that can infiltrate the soil + pot_infiltration = lfr.where((max_sgw -sgw)*por / (self.timestep * self.iterationsData) > pot_infiltration, + pot_infiltration, + (max_sgw -sgw)*por / (self.timestep * self.iterationsData)) + if wet: - # Horton infiltration - pot_infiltration_hort = Ks + (Ki - Ks)*2.71828^(-k*self.T) - pot_infiltration = pot_infiltration_hort * prm * self.ca # m3 that can infiltrate the soil - pot_infiltration = lfr.where((max_sgw -sgw)*por / (self.timestep * self.iterationsData) > pot_infiltration, - pot_infiltration, - (max_sgw -sgw)*por / (self.timestep * self.iterationsData)) - - # Add time to rainfall event self.T = self.T + (self.iterationsData * self.timestep / 3600) @@ -166,6 +166,8 @@ def adjusted_Horton(self, sgw, max_sgw, Ks, Ki, k, prm, por, pre, ev_s): else: infiltration = self.std_arr_lue.zero() self.T = self.T - (self.iterationsData * self.timestep / 3600) - self.T = lfr.where(self.T < 0, 0, self.T) - - return infiltration \ No newline at end of file + self.T = max(self.T, 0) + + pot_reinfiltration = pot_infiltration - infiltration + + return infiltration, pot_reinfiltration \ No newline at end of file diff --git a/model/HBM.py b/model/HBM.py index b52be98..0ed5572 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -65,7 +65,7 @@ def __init__(self, configuration, report): self.cell_area = self.resolution * self.resolution # Retrieve the soil properties - self.Ks, self.porosity, self.wilting_point, self.Ki, self.k = self.retrieve_data.soil_csv( + self.Ks, self.porosity, self.wilting_point, self.Ki = self.retrieve_data.soil_csv( configuration.generalSettings['inputDir'] + configuration.dataSettings['soilData'], soil_type) # soil characteristic # Retrieve the land-use properties @@ -148,21 +148,17 @@ def __init__(self, configuration, report): print("\n") - def update_and_route_both(self, gw_s, gw_flux, discharge, sw_flux): + def update_and_route_both(self, gw_s, gw_flux, discharge, sw_flux, pot_reinfiltration): # The groundwater is adjusted by the fluxes - #channel_infiltation = lfr.where(height > pot_channel_infiltation, pot_channel_infiltation, height) - gw_s = gw_s + gw_flux #+ channel_infiltation*infil_to_gw_s + height = lfr.pow(self.coefficient*discharge, 0.6) + reinfiltration = lfr.where(height > pot_reinfiltration, pot_reinfiltration, height) + gw_s = gw_s + gw_flux + (reinfiltration / self.porosity) # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) - - # Discharge is affected by the surfacewater fluxes, and seepage is added - #height = height + ((sw_flux + seepage)/self.channel_area) #- channel_infiltation - - #discharge = lfr.pow(height, self.c) / self.coefficient - - #discharge = self.height_to_dis(height) + height = height - reinfiltration + discharge = lfr.pow(height, self.c) / self.coefficient # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) @@ -174,14 +170,12 @@ def update_and_route_both(self, gw_s, gw_flux, discharge, sw_flux): self.alpha, self.beta, self.timestep,\ self.channel_length) - #height = self.dis_to_height(discharge) - height = lfr.pow(self.coefficient*discharge, 0.6) # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage gw_s = gw_s - (seepage / self.porosity) - return height, discharge, gw_s, seepage + return height, discharge, gw_s, seepage, reinfiltration def update_and_route_surface(self, height, sw_flux): # Discharge is affected by the surfacewater fluxes, and seepage is added @@ -301,11 +295,10 @@ def dynamic_model(self, configuration): precipitation, evapotranspiration_surface) """ - direct_infiltration = self.calculate_flux.adjusted_Horton(gw_s, + direct_infiltration, pot_reinfiltration = self.calculate_flux.adjusted_Horton(gw_s, self.max_gw_s, self.Ks, self.Ki, - self.k, self.permeability, self.porosity, precipitation, @@ -333,7 +326,7 @@ def dynamic_model(self, configuration): for j in range(dt_data): # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly # Then route the water lateraly and repeat this for the amount of iterations required. - height, discharge, gw_s, seepage = self.update_and_route_both(gw_s, gw_flux, discharge, sw_flux) + height, discharge, gw_s, seepage, reinfiltration = self.update_and_route_both(gw_s, gw_flux, discharge, sw_flux, pot_reinfiltration) # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() @@ -342,7 +335,7 @@ def dynamic_model(self, configuration): # Write value to csv for later validation writer.writerow([i*dt_data + j, outflow]) - variables = {"discharge": discharge, "int_s": int_s, "gw_s": gw_s, "seepage": seepage, "infiltration": direct_infiltration + variables = {"discharge": discharge, "int_s": int_s, "gw_s": gw_s, "seepage": seepage, "infiltration": direct_infiltration, "reinfiltration": reinfiltration, "gw_flow": gw_flow } elif configuration.modelSettings['routing'] == 'surface': diff --git a/model/RetrieveData.py b/model/RetrieveData.py index 86b5343..d5cee26 100644 --- a/model/RetrieveData.py +++ b/model/RetrieveData.py @@ -40,10 +40,11 @@ def soil_csv(self, data_file, soil_type): lpa*: lue partitioned array """ # Assign standard values for de Wupsel - Ks = self.std_arr_lue.one() * 0.05 + Ks = self.std_arr_lue.one() * 0.1 porosity = self.std_arr_lue.one() * 0.35 wilting_point = self.std_arr_lue.one() * 0.15 k = self.std_arr_lue.one() * 1.6 + Ki = self.std_arr_lue.one() * 0.3 # Read pandas data table data_table = pd.read_csv(data_file) @@ -60,7 +61,7 @@ def soil_csv(self, data_file, soil_type): Ki = Ki / 86400 Ks = Ks / 86400 - return Ks, porosity, wilting_point, Ki, k + return Ks, porosity, wilting_point, Ki def land_characteristics_csv(self, data_file, land_use, cell_area): """Reads the land characteristics from a csv file @@ -153,7 +154,7 @@ def csv_timeseries_to_flux(self, data_file, refactor, date): data = pd.read_csv(data_file, sep=",", names=['date_time', 'data_value']) data.set_index('date_time', inplace=True) data_value = data.loc[f'{date_time}']['data_value'] - data_value = data_value * refactor * self.std_arr_lue.one() # Convert to m/s rate from mm/h + data_value = data_value * refactor * self.std_arr_lue.one() # Convert to m3/s rate from mm/h except: data_value = self.std_arr_lue.zero() From f9a4cfb092debf94ba17f8cdbc673e5e9e845d74 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Mon, 12 Jun 2023 13:27:25 +0200 Subject: [PATCH 11/22] Made GIFs --- config/config.ini | 8 ++++---- model/tools/MakeGIF.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/config/config.ini b/config/config.ini index afb6f38..58d24aa 100644 --- a/config/config.ini +++ b/config/config.ini @@ -65,7 +65,7 @@ variables = discharge, seepage, groundWaterHeight, Qgw, Sgw, swFlux, gwFlux, i [gifSettings] variables = discharge, gw_s -fps = 30 -vmin = 0, 20 -vmax = 0.2, 50 -nrRasters = 149 \ No newline at end of file +fps = 10 +vmin = 0, 45 +vmax = 0.5, 50 +nrRasters = 191 \ No newline at end of file diff --git a/model/tools/MakeGIF.py b/model/tools/MakeGIF.py index 59355e1..7b38875 100644 --- a/model/tools/MakeGIF.py +++ b/model/tools/MakeGIF.py @@ -15,7 +15,7 @@ def __init__(self): pass def slice_pathname(pathname, idx, date): - date = date + datetime.timedelta(minutes=idx) + date = date + datetime.timedelta(minutes=15*idx) date_time = date.strftime("%Y-%m-%d-%H%M") return "{}_{}.tiff".format(pathname, date_time) @@ -34,7 +34,7 @@ def create_animation(raster_pathname, nr_rasters, animation_pathname, vmin, vmax image = rasterio.plot.show( data, ax=axis, - cmap="magma", + cmap="Blues", norm = colors.LogNorm(vmin, vmax), ) @@ -44,8 +44,8 @@ def create_animation(raster_pathname, nr_rasters, animation_pathname, vmin, vmax data = np.frombuffer(buffer.getvalue(), dtype=np.uint8) nr_cols, nr_rows = figure.canvas.get_width_height() image = data.reshape(nr_rows, nr_cols, -1) - - writer.append_data(image) + + writer.append_data(image[50:700][50:700]) plt.close() print(f'Processed {i} out of {nr_rasters} rasters.') From b445278094418dba7a80d13c92a3a0ffd648b3df Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Tue, 13 Jun 2023 10:38:03 +0200 Subject: [PATCH 12/22] Improved --- model/CalculateFlux.py | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/model/CalculateFlux.py b/model/CalculateFlux.py index dd5aa78..1cc2920 100644 --- a/model/CalculateFlux.py +++ b/model/CalculateFlux.py @@ -24,7 +24,7 @@ def __init__(self, configuration): self.include_interception = configuration.generalSettings['includeInterception'] self.include_percolation = configuration.generalSettings['includePercolation'] self.ca = int(configuration.modelSettings['resolution'])**2 - self.T = 0 + self.T = self.std_arr_lue.zero() self.k = 1.6 @@ -144,7 +144,7 @@ def adjusted_Horton(self, sgw, max_sgw, Ks, Ki, prm, por, pre, ev_s): infiltration (lpa*): the infiltration rate based on horton (m3/s) """ - wet = lfr.maximum(pre).get() > lfr.maximum(ev_s).get() + wet = pre > ev_s # Horton infiltration pot_infiltration_hort = Ks + (Ki - Ks)*lfr.pow(self.std_arr_lue.one()*2.71828,(-1*self.k*self.T)) @@ -152,22 +152,19 @@ def adjusted_Horton(self, sgw, max_sgw, Ks, Ki, prm, por, pre, ev_s): pot_infiltration = lfr.where((max_sgw -sgw)*por / (self.timestep * self.iterationsData) > pot_infiltration, pot_infiltration, (max_sgw -sgw)*por / (self.timestep * self.iterationsData)) - - if wet: - # Add time to rainfall event - self.T = self.T + (self.iterationsData * self.timestep / 3600) - - # Check if the amount of precipitation available for infiltration is larger than the potential infiltration - enough_rain = pre - ev_s > pot_infiltration - - # Calculate actual infiltration - infiltration = lfr.where(enough_rain, pot_infiltration, pre - ev_s) + + enough_rain = pre - ev_s > pot_infiltration + + self.T = lfr.where(wet, self.T + (self.iterationsData * self.timestep / 3600), + self.T - (self.iterationsData * self.timestep / 3600)) + + self.T = lfr.where(self.T < 0, 0, self.T) + + infiltration1 = lfr.where(enough_rain, pot_infiltration, pre - ev_s) + infiltration2 = lfr.where(infiltration1 < 0, 0, infiltration1) - else: - infiltration = self.std_arr_lue.zero() - self.T = self.T - (self.iterationsData * self.timestep / 3600) - self.T = max(self.T, 0) + - pot_reinfiltration = pot_infiltration - infiltration + pot_reinfiltration = pot_infiltration - infiltration2 - return infiltration, pot_reinfiltration \ No newline at end of file + return infiltration2, pot_reinfiltration \ No newline at end of file From 0debe4953a6aa95392308bc8bd7c23d8d80991c3 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Tue, 13 Jun 2023 10:38:12 +0200 Subject: [PATCH 13/22] Test 1x1 res --- config/config.ini | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/config/config.ini b/config/config.ini index 58d24aa..5172360 100644 --- a/config/config.ini +++ b/config/config.ini @@ -5,7 +5,7 @@ outputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/output/ network = False useAPI = False -scenario = De Hupsel5 +scenario = De Hupsel makeGIF = False @@ -40,9 +40,9 @@ groundWaterBase = 23.25 porosity = 0.35 channel_width = 1 -arrayExtent = 1000 +arrayExtent = 5000 partitionExtent = 1000 -resolution = 5 +resolution = 1 validCellsPercentageD = 24.7 validCellsPercentageG = 25.04 From ffd202344bf5f34c2ddda264eb5b4eb81fa19311 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Tue, 13 Jun 2023 10:38:34 +0200 Subject: [PATCH 14/22] Ported model to new framework --- model/HBM_v2.py | 416 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 416 insertions(+) create mode 100644 model/HBM_v2.py diff --git a/model/HBM_v2.py b/model/HBM_v2.py new file mode 100644 index 0000000..b61348f --- /dev/null +++ b/model/HBM_v2.py @@ -0,0 +1,416 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 15 2023 + +@author: steven.hosper +""" +# The HydrologicBaseModel +import lue.framework as lfr +import math as math +import os +import datetime +import sys +import time +import csv + +# Submodules +from configuration_v2 import Configuration +from reporting import Report +from StandardArraysLUE import StandardArraysLUE +from RetrieveData import RetrieveData +from CalculateFlux import CalculateFlux +from utilityFunctionsHBM import utilityFunctions + +# Other tools +import tools.MakeGIF + +# Timer to add some measure of functionality to the program +start_time = time.time() + +usage = """\ +Run the main model of the hydrologic base model. + +Usage: + {command} + +Options: + {command} : --hpx:thread = integer; + The integer is the amount of cores used during the model run. +""".format( + command=os.path.basename(sys.argv[0]) +) + + +class MyModel(lfr.Model): + def __init__(self, configuration, report): + lfr.Model.__init__(self) + self.configuration = configuration + self.standard_LUE = StandardArraysLUE(configuration) + self.retrieve_data = RetrieveData(configuration) + self.calculate_flux = CalculateFlux(configuration) + self.report = report + + def initialize(self): + # Initialize submodules + + + # Set directories + self.input_dir = self.configuration.generalSettings['inputDir'] + self.configuration.generalSettings['scenario'] + + # Timesteps (all should have same unit, seconds in this case) + self.timestep = float(self.configuration.modelSettings['timestep']) + self.update_timestep = float(self.configuration.modelSettings['iterationsBeforeData']) + self.report_timestep = float(self.configuration.modelSettings['iterationsBeforeReport']) + + partition_shape = 2 * (self.configuration.modelSettings['partitionExtent'],) + + # Initialize data required from memory files + # Get all constants + self.dem = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['dem'], partition_shape) # DEM map of the study area + land_use = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['landUseMap'], partition_shape) # Land-use, example: road + soil_type = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['soilMap'], partition_shape) # example: sand or clay + self.resolution = float(self.configuration.modelSettings['resolution']) + self.cell_area = self.resolution * self.resolution + + # Retrieve the soil properties + self.Ks, self.porosity, self.wilting_point, self.Ki = self.retrieve_data.soil_csv( + self.configuration.generalSettings['inputDir'] + self.configuration.dataSettings['soilData'], soil_type) # soil characteristic + + # Retrieve the land-use properties + self.mannings, self.permeability, self.max_int_s, self.throughfall_frac = \ + self.retrieve_data.land_characteristics_csv( + self.configuration.generalSettings['inputDir'] + self.configuration.dataSettings['landUseData'], land_use, + self.cell_area) # land-use characteristics + + self.gw_base = float(self.configuration.modelSettings['groundWaterBase']) * self.standard_LUE.one() + # self.ldd = lfr.d8_flow_direction(self.dem) + self.ldd = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['ldd'], partition_shape) + + + #### SETTING CONSTANTS #### + ## GENERAL ## + self.routing = self.configuration.modelSettings['routing'] + self.start_date = utilityFunctions.string_to_datetime(self.configuration.modelSettings['startDate'], ", ") + self.stdmin = float(self.configuration.modelSettings['standard_minimal_value']) + slope = self.standard_LUE.one() * 0.008 + self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(slope, 0.05, 0.00001) + self.channel_width = int(self.configuration.modelSettings['channel_width']) + self.coefficient = self.mannings / (0.008 * self.channel_width) + aquifer_height = float(self.configuration.modelSettings['impermeableLayerBelowDEM']) + self.imperm_lay_height = self.dem - aquifer_height + self.water_below_dem = float(self.configuration.modelSettings['waterBelowDEM']) + self.max_gw_s = aquifer_height * self.cell_area # Full storage of porosity + self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point + # Refactorings value from mm/hour to m/s times the cell area. + self.refactor = (self.cell_area / 1000) / 3600 + self.d = (5**(2/3))/(4**(2/3)) + + # Channel length and area + self.channel_length = self.resolution * self.standard_LUE.one() + self.channel_area = self.channel_width * self.channel_length + + # Kinematic Surface Water Routing Constants + self.alpha = 1.5 + self.beta = 0.6 + self.c = 5/3 + + # Load initial groundWaterStorage, if no raster is supplied, use the waterBelowDEM in combination with DEM to create a initialGroundWaterStorage layer. + try: + self.gw_s = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniGroundWaterStorage'], partition_shape) + except: + print("Did not find a initial groundwater height file, looked at: {}".format(self.configuration.dataSettings['iniGroundWaterStorage'])) + self.gw_s = lfr.where(self.dem > (self.gw_base + self.water_below_dem), + ((self.dem - self.water_below_dem)-self.imperm_lay_height) * self.cell_area, + (self.gw_base - self.imperm_lay_height) * self.cell_area) + self.gw_s = lfr.where(self.gw_s > self.max_gw_s, self.max_gw_s, self.gw_s) + + # Load initial discharge, if no raster is supplied, set to zero. + try: + self.height = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniWaterHeight'], partition_shape) + except: + print("Did not find a initial discharge file, looked at: {}".format(self.configuration.dataSettings['iniWaterHeight'])) + self.height = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + try: + self.discharge = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniDischarge'], partition_shape) + except: + print("Did not find a initial discharge file, looked at: {}".format(self.configuration.dataSettings['iniDischarge'])) + self.discharge = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + # Initial InterceptionStorage and groundWaterStorage + try: + self.int_s = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniInterceptionStorage'], partition_shape) + except: + print("Did not find a initial interception storage file, looked at: {}".format(self.configuration.dataSettings['iniInterceptionStorage'])) + self.int_s = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + # REPORTING INITIAL CONDITIONS + variables = {"ini_gw_s": self.gw_s, "ini_int_s": self.int_s, "ini_sur_h": self.height, + } + + self.report.initial(variables) + + + def update_fluxes(self, gw_s, date, update_timestep): + # Load flux and storage values + precipitation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['precipitationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['evapotranspirationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + self.int_s, precipitation, evapotranspiration_surface = self.calculate_flux.interception(self.int_s, + self.max_int_s, + precipitation, + ref_evaporation, + self.throughfall_frac) + + evapotranspiration_surface, evapotranspiration_soil = self.calculate_flux.evapotranspiration(precipitation, + evapotranspiration_surface) + + """ direct_infiltration = self.calculate_flux.infiltration(gw_s, + self.max_gw_s, + self.Ks, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) """ + + direct_infiltration, pot_reinfiltration = self.calculate_flux.adjusted_Horton(gw_s, + self.max_gw_s, + self.Ks, + self.Ki, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) + + # Groundwater LDD, gradient and flow flux + gw_height = self.imperm_lay_height + self.gw_s / self.cell_area + gw_ldd = lfr.d8_flow_direction(gw_height) + del_h_gw = gw_height - lfr.downstream(gw_ldd, gw_height) + gw_grad = (del_h_gw) / self.resolution + gw_flow = self.Ks * gw_grad * self.timestep * (gw_height - self.imperm_lay_height) * self.resolution # Groundwater velocity in m2/s + + # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. + gw_flow = lfr.where(gw_flow * update_timestep > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/update_timestep, gw_flow) + gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) + gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), gw_flow) + + # Add all vertical processes for the surfacewater and all processes groundwater + gw_flux = direct_infiltration - evapotranspiration_soil + lfr.upstream(gw_ldd, gw_flow) - gw_flow + sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters + + return gw_flux, sw_flux, pot_reinfiltration + + def report_output(self): + pass + + def route_both(self, gw_s, gw_flux, discharge, sw_flux, pot_reinfiltration): + """Route surface and subsurface. + + We use the available information in regards to surface and subsurface + layers in order to determine the values of each variable in the next + state. Both layers are updated and routed. + + Args: + gw_s (LUE partitioned array): the groundwater storage in m3 of soil filled with water (m3). + gw_flux (LUE partitioned array): the groundwater flux based on input and flow (m3/s). + discharge (LUE partitioned array): the flow volume that travels through each cell (m3/s). + sw_flux (LUE partitioned array): the surface water flux (m3/s). + pot_reinfiltration (LUE partitioned array): the potential reinfiltration of runoff (m3/s). + + Returns: + height (LUE partitioned array): the water height in each cell (m). + updated_discharge (LUE partitioned array): the updated flow volume that travels through each cell (m3/s). + updated_gw_s (LUE partitioned array): the updated groundwater storage in m3 of soil filled with water. (m3) + seepage (LUE partitioned array): the amount of water leaving the soil (m3/s). + reinfiltration (LUE partitioned array): the amount of runoff that infiltrates the soil (m3/s). + """ + # The groundwater is adjusted by the fluxes + height = lfr.pow(self.coefficient*discharge, 0.6) + reinfiltration = lfr.where(height > pot_reinfiltration, pot_reinfiltration, height) + gw_s = gw_s + gw_flux + (reinfiltration / self.porosity) + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) + + height = height - reinfiltration + discharge = lfr.pow(height, self.c) / self.coefficient + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + inflow = (sw_flux + seepage)/self.channel_area + inflow = lfr.where(inflow < self.stdmin, self.stdmin, inflow) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, inflow,\ + self.alpha, self.beta, self.timestep,\ + self.channel_length) + + height = lfr.pow(self.coefficient*discharge, 0.6) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + gw_s = gw_s - (seepage / self.porosity) + + return height, discharge, gw_s, seepage, reinfiltration + + def route_surface(self, height, sw_flux): + """Route only the surface + + The height and surface water flux are used to calculate a new + water table height, and route water in downstream direction. + + Args: + height (LUE partitioned array): the old water height (m). + sw_flux (LUE partitioned array): the surface water flux (m3/s). + + Returns: + updated_height (LUE partitioned array): the water height in each cell (m). + discharge (LUE partitioned array): the flow volume that travels through each cell (m3/s). + """ + # Discharge is affected by the surfacewater fluxes, and seepage is added + new_height = height + ((sw_flux)/self.channel_area) #- channel_infiltation + + discharge = lfr.pow(new_height, self.c) / self.coefficient + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, self.stdmin,\ + self.alpha, self.beta, self.timestep,\ + self.channel_length) + + updated_height = lfr.pow(self.coefficient*discharge, 0.6) + + return updated_height, discharge + + def route_subsurface(self, gw_s, gw_flux): + """Route only the subsurface + + The subsurface is routed based on the supplied groundwater storage, + and the given groundwater flux. The new groundwater storage and the + seepage are returned by the function. + + Args: + gw_s (LUE partitioned array): The groundwater storage in m3 of soil filled with water. + gw_flux (LUE partitioned array): The groundwater flux based on input and flow (m3/s). + + Returns: + updated_gw_s (LUE partitioned array): The updated groundwater storage in m3 of soil filled with water. + seepage (LUE partitioned array): The amount of water leaving the soil (m3/s). + """ + # The groundwater is adjusted by the fluxes + new_gw_s = gw_s + gw_flux + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(new_gw_s > self.max_gw_s, (new_gw_s - self.max_gw_s)*self.porosity, 0) + + updated_gw_s = new_gw_s - (seepage / self.porosity) + + return updated_gw_s, seepage + + def simulate(self, time_step): + date = self.start_date + datetime.timedelta(seconds=time_step) + + # Get new fluxes + if time_step % self.update_timestep: + self.gw_flux, self.sw_flux, self.pot_reinfiltration = self.update_fluxes(self.gw_s, date, self.update_timestep) + + # Update all variables and route. + if self.routing == "both": + self.height, self.discharge, self.gw_s, self.seepage, self.reinfiltration = self.route_both( + self.gw_s, + self.gw_flux, + self.discharge, + self.sw_flux, + self.pot_reinfiltration + ) + + variables = {"discharge": self.discharge, "int_s": self.int_s, + "gw_s": self.gw_s, "seepage": self.seepage, + "reinfiltration": self.reinfiltration + } + + elif self.routing == 'surface': + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + self.height, self.discharge = self.update_and_route_surface(self.height, + self.sw_flux + ) + + variables = {"discharge": self.discharge, "int_s": self.int_s, + "height": self.height + } + + + elif self.routing == 'subsurface': + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + self.gw_s, self.seepage = self.update_and_route_subsurface(self.gw_s, self.gw_flux) + + variables = {"gw_s": self.gw_s, "seepage": self.seepage, + "groundwater_height": self.gw_height, + } + + else: + raise("Unknown 'routing' mode, change configuration setting to one of the following: both, surface, subsurface") + + # Report output + if time_step % self.report_timestep == 0: + self.report.dynamic(date, variables) + + + + def finalize(self): + # Create the balance report + self.report.balance_report(self.configuration) + + # Process the results into a gif + if self.configuration.generalSettings['makeGIF'] == 'True': + print(f"Creating a GIF for: {self.configuration.gifSettings['variables']}.") + tools.MakeGIF.run(self.configuration) + +class MyProgressor(lfr.Progressor): + def __init__(self): + lfr.Progressor.__init__(self) + + + def initialize(self): + sys.stdout.write("[") + sys.stdout.flush() + + def simulate(self, time_step): + sys.stdout.write(".") + sys.stdout.flush() + + def finalize(self): + sys.stdout.write("]\n") + sys.stdout.flush() + +def calculate_timesteps(start_date: str, end_date: str, sep, timestep_size): + sd = utilityFunctions.string_to_datetime(start_date, sep) + ed = utilityFunctions.string_to_datetime(end_date, sep) + timesteps = int((ed - sd).total_seconds() / int(timestep_size)) + return timesteps + +@lfr.runtime_scope +def main(): + configuration = Configuration("F:/Projecten intern (2023)/Stage Steven Hosper/Model/v1/config/config.ini") + report = Report(configuration) + model = MyModel(configuration, report) + progressor = MyProgressor() + nr_time_steps = calculate_timesteps(configuration.modelSettings['startDate'], + configuration.modelSettings['endDate'], + ", ", + configuration.modelSettings['timestep'] + ) + lfr.run_deterministic(model, progressor, nr_time_steps, rate_limit=15) + +if __name__ == "__main__": + sys.exit(main()) + +print("--- %s seconds ---" % (time.time() - start_time)) From 10b4eebeae58ad96e42b43cac11a9706f36a8266 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Tue, 13 Jun 2023 10:38:58 +0200 Subject: [PATCH 15/22] Saved files to v2 folder for tests. --- model/HBM.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/model/HBM.py b/model/HBM.py index 0ed5572..ba87ee2 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -51,7 +51,7 @@ def __init__(self, configuration, report): # Set directories self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] - self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings['scenario'] + self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings['scenario'] + " v2" partition_shape = 2 * (configuration.modelSettings['partitionExtent'],) @@ -345,8 +345,8 @@ def dynamic_model(self, configuration): height, discharge = self.update_and_route_surface(height, sw_flux) # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) - outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() - total_volume = lfr.minimum(lfr.zonal_sum(height, self.dem>0)).get() + outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)) + total_volume = lfr.minimum(lfr.zonal_sum(height, self.dem>0)) print("outflow: ", outflow, "total: ", total_volume) # Write value to csv for later validation @@ -363,7 +363,7 @@ def dynamic_model(self, configuration): gw_s, seepage = self.update_and_route_subsurface(gw_s, gw_flux) outflow = lfr.maximum(lfr.zonal_sum(seepage, self.dem>0)).get() / self.timestep - total_volume = lfr.minimum(lfr.zonal_sum(gw_s, self.dem>0)).get() + total_volume = lfr.minimum(lfr.zonal_sum(gw_s, self.dem>0)) # Write value to csv for later validation writer.writerow([i*60 + j, outflow]) From 771c7b0809ba70ce73b4fea7fbabec11a7f76b76 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Tue, 13 Jun 2023 10:39:14 +0200 Subject: [PATCH 16/22] Reported to v2 for tests --- model/reporting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/reporting.py b/model/reporting.py index 958a9fb..1eb07f4 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -17,7 +17,7 @@ def __init__(self, configuration): self.standard_LUE = StandardArraysLUE(configuration) self.timestep = int(configuration.modelSettings['timestep']) self.iterations = int(configuration.modelSettings['iterationsBeforeReport']) - self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings["scenario"] + self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings["scenario"] + " v2" self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] def initial(self, variables): From 89ddbd5dbafaef076465ea4399d85144d865aa0f Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Wed, 14 Jun 2023 08:28:36 +0200 Subject: [PATCH 17/22] Computational tests --- .vscode/launch.json | 2 +- config/config_computational_1.ini | 71 +++++++++++++++++++++++++++++++ config/config_computational_2.ini | 71 +++++++++++++++++++++++++++++++ config/config_computational_3.ini | 71 +++++++++++++++++++++++++++++++ 4 files changed, 214 insertions(+), 1 deletion(-) create mode 100644 config/config_computational_1.ini create mode 100644 config/config_computational_2.ini create mode 100644 config/config_computational_3.ini diff --git a/.vscode/launch.json b/.vscode/launch.json index 4454da8..ea30711 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -2,7 +2,7 @@ // Use IntelliSense to learn about possible attributes. // Hover to view descriptions of existing attributes. // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 - "version": "0.3.6", + "version": "0.3.7", "configurations": [ { "name": "Python: Current File", diff --git a/config/config_computational_1.ini b/config/config_computational_1.ini new file mode 100644 index 0000000..f580215 --- /dev/null +++ b/config/config_computational_1.ini @@ -0,0 +1,71 @@ +[generalSettings] +# Directories +inputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/data/ +outputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/output/ + +network = False +useAPI = False +scenario = De Hupsel + +makeGIF = False + +# Include Processes +includePrecipitation = True +includeEvapotranspiration = True +includeInfiltration = True +includeInterception = True +includePercolation = False + +[apiSettings] +username = __key__ +password = Cy0BNm8p.vpytC2vYPT9g7OKdgxvqggyV0k9zzJVy +precipAPI = 730d6675-35dd-4a35-aa9b-bfb8155f9ca7 +evapAPI = e262dc03-f12b-4082-a4f4-7d0534e31fa4 +demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 + + +[modelSettings] +# Date = y, m, d, h, m, s +startDate = 2023, 4, 6, 12, 30, 0 +endDate = 2023, 4, 6, 12, 32, 30 +iterationsBeforeData = 300 +iterationsBeforeReport = 900 +timestep = 1 +standard_minimal_value = 1E-20 +routing = both + +waterBelowDEM = 0.0 +impermeableLayerBelowDEM= 2.00 +groundWaterBase = 23.25 +porosity = 0.35 +channel_width = 1 + +arrayExtent = 5000 +partitionExtent = 5000 +resolution = 1 +validCellsPercentageD = 24.7 +validCellsPercentageG = 25.04 + +[dataSettings] +iniGroundWaterStorage = +iniWaterHeight = +iniDischarge = +iniInterceptionStorage = +dem = /dem_f64_pcr_shaped.tiff +ldd = /ldd_f64_pcr_shaped.tiff +soilMap = /bodem.tiff +landUseMap = /landgebruik.tiff +soilData = /soil_conversion.csv +landUseData = /landuse_conversion_v2.csv +precipitationData = /precipitation.csv +evapotranspirationData = /evapotranspiration.csv + +[reportSettings] +variables = discharge, seepage, groundWaterHeight, Qgw, Sgw, swFlux, gwFlux, infiltration, evapotranspirationSoil + +[gifSettings] +variables = discharge, gw_s +fps = 10 +vmin = 0, 45 +vmax = 0.5, 50 +nrRasters = 191 \ No newline at end of file diff --git a/config/config_computational_2.ini b/config/config_computational_2.ini new file mode 100644 index 0000000..445ef05 --- /dev/null +++ b/config/config_computational_2.ini @@ -0,0 +1,71 @@ +[generalSettings] +# Directories +inputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/data/ +outputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/output/ + +network = False +useAPI = False +scenario = De Hupsel + +makeGIF = False + +# Include Processes +includePrecipitation = True +includeEvapotranspiration = True +includeInfiltration = True +includeInterception = True +includePercolation = False + +[apiSettings] +username = __key__ +password = Cy0BNm8p.vpytC2vYPT9g7OKdgxvqggyV0k9zzJVy +precipAPI = 730d6675-35dd-4a35-aa9b-bfb8155f9ca7 +evapAPI = e262dc03-f12b-4082-a4f4-7d0534e31fa4 +demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 + + +[modelSettings] +# Date = y, m, d, h, m, s +startDate = 2023, 4, 6, 12, 30, 0 +endDate = 2023, 4, 6, 12, 45, 0 +iterationsBeforeData = 300 +iterationsBeforeReport = 900 +timestep = 1 +standard_minimal_value = 1E-20 +routing = both + +waterBelowDEM = 0.0 +impermeableLayerBelowDEM= 2.00 +groundWaterBase = 23.25 +porosity = 0.35 +channel_width = 1 + +arrayExtent = 5000 +partitionExtent = 5000 +resolution = 1 +validCellsPercentageD = 24.7 +validCellsPercentageG = 25.04 + +[dataSettings] +iniGroundWaterStorage = +iniWaterHeight = +iniDischarge = +iniInterceptionStorage = +dem = /dem_f64_pcr_shaped.tiff +ldd = /ldd_f64_pcr_shaped.tiff +soilMap = /bodem.tiff +landUseMap = /landgebruik.tiff +soilData = /soil_conversion.csv +landUseData = /landuse_conversion_v2.csv +precipitationData = /precipitation.csv +evapotranspirationData = /evapotranspiration.csv + +[reportSettings] +variables = discharge, seepage, groundWaterHeight, Qgw, Sgw, swFlux, gwFlux, infiltration, evapotranspirationSoil + +[gifSettings] +variables = discharge, gw_s +fps = 10 +vmin = 0, 45 +vmax = 0.5, 50 +nrRasters = 191 \ No newline at end of file diff --git a/config/config_computational_3.ini b/config/config_computational_3.ini new file mode 100644 index 0000000..445ef05 --- /dev/null +++ b/config/config_computational_3.ini @@ -0,0 +1,71 @@ +[generalSettings] +# Directories +inputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/data/ +outputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/output/ + +network = False +useAPI = False +scenario = De Hupsel + +makeGIF = False + +# Include Processes +includePrecipitation = True +includeEvapotranspiration = True +includeInfiltration = True +includeInterception = True +includePercolation = False + +[apiSettings] +username = __key__ +password = Cy0BNm8p.vpytC2vYPT9g7OKdgxvqggyV0k9zzJVy +precipAPI = 730d6675-35dd-4a35-aa9b-bfb8155f9ca7 +evapAPI = e262dc03-f12b-4082-a4f4-7d0534e31fa4 +demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 + + +[modelSettings] +# Date = y, m, d, h, m, s +startDate = 2023, 4, 6, 12, 30, 0 +endDate = 2023, 4, 6, 12, 45, 0 +iterationsBeforeData = 300 +iterationsBeforeReport = 900 +timestep = 1 +standard_minimal_value = 1E-20 +routing = both + +waterBelowDEM = 0.0 +impermeableLayerBelowDEM= 2.00 +groundWaterBase = 23.25 +porosity = 0.35 +channel_width = 1 + +arrayExtent = 5000 +partitionExtent = 5000 +resolution = 1 +validCellsPercentageD = 24.7 +validCellsPercentageG = 25.04 + +[dataSettings] +iniGroundWaterStorage = +iniWaterHeight = +iniDischarge = +iniInterceptionStorage = +dem = /dem_f64_pcr_shaped.tiff +ldd = /ldd_f64_pcr_shaped.tiff +soilMap = /bodem.tiff +landUseMap = /landgebruik.tiff +soilData = /soil_conversion.csv +landUseData = /landuse_conversion_v2.csv +precipitationData = /precipitation.csv +evapotranspirationData = /evapotranspiration.csv + +[reportSettings] +variables = discharge, seepage, groundWaterHeight, Qgw, Sgw, swFlux, gwFlux, infiltration, evapotranspirationSoil + +[gifSettings] +variables = discharge, gw_s +fps = 10 +vmin = 0, 45 +vmax = 0.5, 50 +nrRasters = 191 \ No newline at end of file From 987321bfd6a457e2f91bb5da098b6fd8a4662fb5 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Wed, 21 Jun 2023 09:58:30 +0200 Subject: [PATCH 18/22] Updates & computational tests --- config/config_computational_1.ini | 4 +- config/config_computational_3.ini | 2 +- model/HBM.py | 4 +- model/HBM_computational_tests.py | 429 +++++++++++++++++++++++++++ model/HBM_v2.py | 8 +- model/HBM_v2_computational_tests.py | 437 ++++++++++++++++++++++++++++ model/reporting.py | 2 +- 7 files changed, 875 insertions(+), 11 deletions(-) create mode 100644 model/HBM_computational_tests.py create mode 100644 model/HBM_v2_computational_tests.py diff --git a/config/config_computational_1.ini b/config/config_computational_1.ini index f580215..a01921e 100644 --- a/config/config_computational_1.ini +++ b/config/config_computational_1.ini @@ -27,7 +27,7 @@ demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 [modelSettings] # Date = y, m, d, h, m, s startDate = 2023, 4, 6, 12, 30, 0 -endDate = 2023, 4, 6, 12, 32, 30 +endDate = 2023, 4, 6, 12, 30, 50 iterationsBeforeData = 300 iterationsBeforeReport = 900 timestep = 1 @@ -41,7 +41,7 @@ porosity = 0.35 channel_width = 1 arrayExtent = 5000 -partitionExtent = 5000 +partitionExtent = 1000 resolution = 1 validCellsPercentageD = 24.7 validCellsPercentageG = 25.04 diff --git a/config/config_computational_3.ini b/config/config_computational_3.ini index 445ef05..7fe04d9 100644 --- a/config/config_computational_3.ini +++ b/config/config_computational_3.ini @@ -41,7 +41,7 @@ porosity = 0.35 channel_width = 1 arrayExtent = 5000 -partitionExtent = 5000 +partitionExtent = 1000 resolution = 1 validCellsPercentageD = 24.7 validCellsPercentageG = 25.04 diff --git a/model/HBM.py b/model/HBM.py index ba87ee2..676faf6 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -22,7 +22,7 @@ from utilityFunctionsHBM import utilityFunctions # Other tools -import tools.MakeGIF +#import tools.MakeGIF # Timer to add some measure of functionality to the program start_time = time.time() @@ -424,6 +424,6 @@ def dynamic_model(self, configuration): # Process the results into a gif if configuration.generalSettings['makeGIF'] == 'True': print(f"Creating a GIF for: {configuration.gifSettings['variables']}.") - tools.MakeGIF.run(configuration) + #tools.MakeGIF.run(configuration) print("--- %s seconds ---" % (time.time() - start_time)) diff --git a/model/HBM_computational_tests.py b/model/HBM_computational_tests.py new file mode 100644 index 0000000..676faf6 --- /dev/null +++ b/model/HBM_computational_tests.py @@ -0,0 +1,429 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 15 2023 + +@author: steven.hosper +""" +# The HydrologicBaseModel +import lue.framework as lfr +import math as math +import os +import datetime +import sys +import time +import csv + +# Submodules +from configuration_v2 import Configuration +from reporting import Report +from StandardArraysLUE import StandardArraysLUE +from RetrieveData import RetrieveData +from CalculateFlux import CalculateFlux +from utilityFunctionsHBM import utilityFunctions + +# Other tools +#import tools.MakeGIF + +# Timer to add some measure of functionality to the program +start_time = time.time() + +usage = """\ +Run the main model of the hydrologic base model. + +Usage: + {command} + +Options: + {command} : --hpx:thread = integer; + The integer is the amount of cores used during the model run. +""".format( + command=os.path.basename(sys.argv[0]) +) + +class mainModel: + def __init__(self, configuration, report): + print("Initializing the program...") + # Initialize submodules + self.standard_LUE = StandardArraysLUE(configuration) + self.retrieve_data = RetrieveData(configuration) + self.calculate_flux = CalculateFlux(configuration) + self.report = report + + # Set directories + self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] + self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings['scenario'] + " v2" + + partition_shape = 2 * (configuration.modelSettings['partitionExtent'],) + + # Initialize data required from memory files + # Get all constants + self.dem = lfr.from_gdal(self.input_dir + configuration.dataSettings['dem'], partition_shape) # DEM map of the study area + self.dem = lfr.where(self.dem < 0.1, 35, self.dem) + land_use = lfr.from_gdal(self.input_dir + configuration.dataSettings['landUseMap'], partition_shape) # Land-use, example: road + soil_type = lfr.from_gdal(self.input_dir + configuration.dataSettings['soilMap'], partition_shape) # example: sand or clay + self.resolution = float(configuration.modelSettings['resolution']) + self.cell_area = self.resolution * self.resolution + + # Retrieve the soil properties + self.Ks, self.porosity, self.wilting_point, self.Ki = self.retrieve_data.soil_csv( + configuration.generalSettings['inputDir'] + configuration.dataSettings['soilData'], soil_type) # soil characteristic + + # Retrieve the land-use properties + self.mannings, self.permeability, self.max_int_s, self.throughfall_frac = \ + self.retrieve_data.land_characteristics_csv( + configuration.generalSettings['inputDir'] + configuration.dataSettings['landUseData'], land_use, + self.cell_area) # land-use characteristics + + self.gw_base = float(configuration.modelSettings['groundWaterBase']) * self.standard_LUE.one() + # self.ldd = lfr.d8_flow_direction(self.dem) + self.ldd = lfr.from_gdal(self.input_dir + configuration.dataSettings['ldd'], partition_shape) + + + #### SETTING CONSTANTS #### + ## GENERAL ## + self.stdmin = float(configuration.modelSettings['standard_minimal_value']) + self.slope = self.standard_LUE.one() * 0.008 + self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(self.slope, 0.05, 0.00001) + self.channel_width = int(configuration.modelSettings['channel_width']) + self.coefficient = self.mannings / (0.008 * self.channel_width) + self.imperm_below_dem = float(configuration.modelSettings['impermeableLayerBelowDEM']) + self.imperm_lay_height = self.dem - self.imperm_below_dem + self.water_below_dem = float(configuration.modelSettings['waterBelowDEM']) + self.max_gw_s = self.imperm_below_dem * self.cell_area # Full storage of porosity + self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point + # Refactorings value from mm/hour to m/s times the cell area. + self.refactor = (self.cell_area / 1000) / 3600 + self.d = (5**(2/3))/(4**(2/3)) + + # Channel length and area + self.channel_length = self.resolution * self.standard_LUE.one() + self.channel_area = self.channel_width * self.channel_length + self.channel_rat = self.channel_area / self.cell_area + self.infil_to_gw_s = self.channel_area / self.porosity + # self.notBoundaryCells = generate.boundaryCell() # Currently not working + + # Kinematic Surface Water Routing Constants + self.alpha = 1.5 + self.beta = 0.6 + self.timestep = 1.0 * float(configuration.modelSettings['timestep']) + self.c = 5/3 + + # Static, really small value because inflow = 0 is not accepted + self.inflow = self.standard_LUE.one()*self.stdmin + + # Load initial groundWaterStorage, if no raster is supplied, use the waterBelowDEM in combination with DEM to create a initialGroundWaterStorage layer. + try: + self.ini_gw_s = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniGroundWaterStorage'], partition_shape) + except: + print("Did not find a initial groundwater height file, looked at: {}".format(configuration.dataSettings['iniGroundWaterStorage'])) + self.ini_gw_s = lfr.where(self.dem > (self.gw_base + self.water_below_dem), + ((self.dem - self.water_below_dem)-self.imperm_lay_height) * self.cell_area, + (self.gw_base - self.imperm_below_dem) * self.cell_area) + self.ini_gw_s = lfr.where(self.ini_gw_s > self.max_gw_s, self.max_gw_s, self.ini_gw_s) + + # Load initial discharge, if no raster is supplied, set to zero. + try: + self.ini_sur_h = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniWaterHeight'], partition_shape) + except: + print("Did not find a initial discharge file, looked at: {}".format(configuration.dataSettings['iniWaterHeight'])) + self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + try: + self.ini_dis = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniDischarge'], partition_shape) + except: + print("Did not find a initial discharge file, looked at: {}".format(configuration.dataSettings['iniDischarge'])) + self.ini_dis = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + # Initial InterceptionStorage and groundWaterStorage + try: + self.ini_int_s = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniInterceptionStorage'], partition_shape) + except: + print("Did not find a initial interception storage file, looked at: {}".format(configuration.dataSettings['iniInterceptionStorage'])) + self.ini_int_s = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + # REPORTING INITIAL CONDITIONS + variables = {"ini_gw_s": self.ini_gw_s, "ini_int_s": self.ini_int_s, "ini_sur_h": self.ini_sur_h, + } + self.report.initial(variables) + + print("\n") + + def update_and_route_both(self, gw_s, gw_flux, discharge, sw_flux, pot_reinfiltration): + # The groundwater is adjusted by the fluxes + height = lfr.pow(self.coefficient*discharge, 0.6) + reinfiltration = lfr.where(height > pot_reinfiltration, pot_reinfiltration, height) + gw_s = gw_s + gw_flux + (reinfiltration / self.porosity) + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) + + height = height - reinfiltration + discharge = lfr.pow(height, self.c) / self.coefficient + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + inflow = (sw_flux + seepage)/self.channel_area + inflow = lfr.where(inflow < self.stdmin, self.stdmin, inflow) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, inflow,\ + self.alpha, self.beta, self.timestep,\ + self.channel_length) + + height = lfr.pow(self.coefficient*discharge, 0.6) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + gw_s = gw_s - (seepage / self.porosity) + + return height, discharge, gw_s, seepage, reinfiltration + + def update_and_route_surface(self, height, sw_flux): + # Discharge is affected by the surfacewater fluxes, and seepage is added + height = height + ((sw_flux)/self.channel_area) #- channel_infiltation + + discharge = lfr.pow(height, self.c) / self.coefficient + + #discharge = self.height_to_dis(height) + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, self.inflow,\ + self.alpha, self.beta, self.timestep,\ + self.channel_length) + + height = lfr.pow(self.coefficient*discharge, 0.6) + + #height = self.dis_to_height(discharge) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + return height, discharge + + def update_and_route_subsurface(self, gw_s, gw_flux): + # The groundwater is adjusted by the fluxes + gw_s = gw_s + gw_flux + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) + + gw_s = gw_s - (seepage / self.porosity) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + + return gw_s, seepage + + def dis_to_height(self, discharge): + """Use the open channel flow formula to convert discharge to height + + The Open Channel Flow calculation for rectangular channels using the known width, slope and mannings. + + Args: + discharge (lue partitioned array): discharge array that is a result of the kinematic wave formule. + + Returns: + height (lue partitioned array): height array that can be used to add vertical and lateral in- and outflows. + + """ + height = (1/self.d)*self.mannings * discharge / self.slope_sqrd + + return height + + def height_to_dis(self, height): + """Use the open channel flow formula to convert discharge to height + + The Open Channel Flow calculation for rectangular channels using the known width, slope and mannings. + + Args: + height (lue partitioned array): height array after all in- and outfluxes have been added. + + Returns: + discharge (lue partitioned array): discharge array that can be used as an input for the kinematic wave formula. + + """ + discharge = self.d * (height * self.slope_sqrd)/self.mannings + + return discharge + + @lfr.runtime_scope + def dynamic_model(self, configuration): + dt_data = int(configuration.modelSettings['iterationsBeforeData']) + dt_report = int(configuration.modelSettings['iterationsBeforeReport']) + start_date = utilityFunctions.string_to_datetime(configuration.modelSettings['startDate'], ", ") + end_date = utilityFunctions.string_to_datetime(configuration.modelSettings['endDate'], ", ") + dT = int((end_date - start_date).total_seconds() / (dt_data*self.timestep)) + + # Loading initial conditions + height = self.ini_sur_h + gw_s = self.ini_gw_s + int_s = self.ini_int_s + gw_height = self.imperm_lay_height + gw_s/self.cell_area + discharge = self.ini_dis + + # Open file to write maximum discharge values to for post simulation validation. + with open(self.output_dir + "/maximumDischarge.csv", "w", newline="") as f: + writer = csv.writer(f, delimiter=';') + + # Start model for dT large periods + for i in range(dT): + # Time in minutes is the small iteration multiplied with the timestep (both in seconds) divided by 60 seconds. + date = start_date + datetime.timedelta(seconds = i * (dt_data*self.timestep)) + + # Load flux and storage values + precipitation = self.retrieve_data.csv_timeseries_to_flux(configuration.generalSettings['inputDir'] + + configuration.dataSettings['precipitationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(configuration.generalSettings['inputDir'] + + configuration.dataSettings['evapotranspirationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + int_s, precipitation, evapotranspiration_surface = self.calculate_flux.interception(int_s, + self.max_int_s, + precipitation, + ref_evaporation, + self.throughfall_frac) + + evapotranspiration_surface, evapotranspiration_soil = self.calculate_flux.evapotranspiration(precipitation, + evapotranspiration_surface) + + """ direct_infiltration = self.calculate_flux.infiltration(gw_s, + self.max_gw_s, + self.Ks, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) """ + + direct_infiltration, pot_reinfiltration = self.calculate_flux.adjusted_Horton(gw_s, + self.max_gw_s, + self.Ks, + self.Ki, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) + + # Groundwater LDD, gradient and flow flux + gw_ldd = lfr.d8_flow_direction(gw_height) + del_h_gw = gw_height - lfr.downstream(gw_ldd, gw_height) + gw_grad = (del_h_gw) / self.resolution + gw_flow = self.Ks * gw_grad * self.timestep * (gw_height - self.imperm_lay_height) * self.resolution # Groundwater velocity in m2/s + + # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. + gw_flow = lfr.where(gw_flow * dt_data > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/dt_data, gw_flow) + gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) + gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), gw_flow) + + # Add all vertical processes for the surfacewater and all processes groundwater + gw_flux = direct_infiltration - evapotranspiration_soil + lfr.upstream(gw_ldd, gw_flow) - gw_flow + sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters + + + ############## ROUTING AND UPDATE PER TIMESTEP ############## + + if configuration.modelSettings['routing'] == 'both': + for j in range(dt_data): + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + height, discharge, gw_s, seepage, reinfiltration = self.update_and_route_both(gw_s, gw_flux, discharge, sw_flux, pot_reinfiltration) + + # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) + outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() + print("outflow: ", outflow) + + # Write value to csv for later validation + writer.writerow([i*dt_data + j, outflow]) + + variables = {"discharge": discharge, "int_s": int_s, "gw_s": gw_s, "seepage": seepage, "infiltration": direct_infiltration, "reinfiltration": reinfiltration, "gw_flow": gw_flow + } + + elif configuration.modelSettings['routing'] == 'surface': + for j in range(dt_data): + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + height, discharge = self.update_and_route_surface(height, sw_flux) + + # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) + outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)) + total_volume = lfr.minimum(lfr.zonal_sum(height, self.dem>0)) + print("outflow: ", outflow, "total: ", total_volume) + + # Write value to csv for later validation + writer.writerow([i*60 + j, outflow]) + + variables = {"discharge": discharge, "int_s": int_s, "height": height + } + + + elif configuration.modelSettings['routing'] == 'subsurface': + for j in range(dt_data): + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + gw_s, seepage = self.update_and_route_subsurface(gw_s, gw_flux) + + outflow = lfr.maximum(lfr.zonal_sum(seepage, self.dem>0)).get() / self.timestep + total_volume = lfr.minimum(lfr.zonal_sum(gw_s, self.dem>0)) + + # Write value to csv for later validation + writer.writerow([i*60 + j, outflow]) + + variables = {"gw_s": gw_s, "seepage": seepage, "groundwater_height": gw_height, + } + + else: + raise("Unknown 'routing' mode, change configuration setting to one of the following: both, surface, subsurface") + break + + + + # Adjust the GW Table for the LDD creation of the next timestep. + gw_height = self.imperm_lay_height + gw_s/self.cell_area + + + # Save / Report data + if i % (dt_report/dt_data) == 0: + print(i) + report.dynamic(date, variables) + + + print(f"Done: {i+1}/{dT}") + return 0 + + + +# Initialize HPX runtime and run model, on the root locality ------------------- +# General configuration options, which are valid on all +# platforms. Platform-specific options can be passed on the command line. +cfg = [ + # Make sure hpx_main is always executed + "hpx.run_hpx_main!=1", + # Allow for unknown command line options + "hpx.commandline.allow_unknown!=1", + # Disable HPX' short options + "hpx.commandline.aliasing!=0", + # Don't print diagnostics during forced terminate + "hpx.diagnostics_on_terminate!=1", + # Make AGAS clean up resources faster than by default + "hpx.agas.max_pending_refcnt_requests!=50", +] + +lfr.start_hpx_runtime(cfg) + +# The root locality will distribute the work over all other +# localities. Never perform Python code on the other localities than the +# root locality unless you know what you are doing. +if lfr.on_root_locality(): + # Run the main model + configuration = Configuration("F:/Projecten intern (2023)/Stage Steven Hosper/Model/v1/config/config.ini") + report = Report(configuration) + main = mainModel(configuration, report) + main.dynamic_model(configuration) + report.balance_report(configuration) + + # Process the results into a gif + if configuration.generalSettings['makeGIF'] == 'True': + print(f"Creating a GIF for: {configuration.gifSettings['variables']}.") + #tools.MakeGIF.run(configuration) + +print("--- %s seconds ---" % (time.time() - start_time)) diff --git a/model/HBM_v2.py b/model/HBM_v2.py index b61348f..e52d98b 100644 --- a/model/HBM_v2.py +++ b/model/HBM_v2.py @@ -22,7 +22,7 @@ from utilityFunctionsHBM import utilityFunctions # Other tools -import tools.MakeGIF +#import tools.MakeGIF # Timer to add some measure of functionality to the program start_time = time.time() @@ -51,9 +51,6 @@ def __init__(self, configuration, report): self.report = report def initialize(self): - # Initialize submodules - - # Set directories self.input_dir = self.configuration.generalSettings['inputDir'] + self.configuration.generalSettings['scenario'] @@ -101,6 +98,7 @@ def initialize(self): self.water_below_dem = float(self.configuration.modelSettings['waterBelowDEM']) self.max_gw_s = aquifer_height * self.cell_area # Full storage of porosity self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point + # Refactorings value from mm/hour to m/s times the cell area. self.refactor = (self.cell_area / 1000) / 3600 self.d = (5**(2/3))/(4**(2/3)) @@ -372,7 +370,7 @@ def finalize(self): # Process the results into a gif if self.configuration.generalSettings['makeGIF'] == 'True': print(f"Creating a GIF for: {self.configuration.gifSettings['variables']}.") - tools.MakeGIF.run(self.configuration) + #tools.MakeGIF.run(self.configuration) class MyProgressor(lfr.Progressor): def __init__(self): diff --git a/model/HBM_v2_computational_tests.py b/model/HBM_v2_computational_tests.py new file mode 100644 index 0000000..104c8d8 --- /dev/null +++ b/model/HBM_v2_computational_tests.py @@ -0,0 +1,437 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 15 2023 + +@author: steven.hosper +""" +# The HydrologicBaseModel +import lue.framework as lfr +import math as math +import os +import datetime +import sys +import time +import csv +import numpy as np + +# Submodules +from configuration_v2 import Configuration +from reporting import Report +from StandardArraysLUE import StandardArraysLUE +from RetrieveData import RetrieveData +from CalculateFlux import CalculateFlux +from utilityFunctionsHBM import utilityFunctions + +# Other tools +#import tools.MakeGIF + +# Timer to add some measure of functionality to the program +start_time = time.time() + +usage = """\ +Run the main model of the hydrologic base model. + +Usage: + {command} + +Options: + {command} : --hpx:thread = integer; + The integer is the amount of cores used during the model run. +""".format( + command=os.path.basename(sys.argv[0]) +) + + +class MyModel(lfr.Model): + def __init__(self, configuration, report): + lfr.Model.__init__(self) + self.configuration = configuration + self.standard_LUE = StandardArraysLUE(configuration) + self.retrieve_data = RetrieveData(configuration) + self.calculate_flux = CalculateFlux(configuration) + self.report = report + + def initialize(self): + # Set directories + self.input_dir = self.configuration.generalSettings['inputDir'] + self.configuration.generalSettings['scenario'] + + # Timesteps (all should have same unit, seconds in this case) + self.timestep = float(self.configuration.modelSettings['timestep']) + self.update_timestep = float(self.configuration.modelSettings['iterationsBeforeData']) + self.report_timestep = float(self.configuration.modelSettings['iterationsBeforeReport']) + + partition_shape = 2 * (self.configuration.modelSettings['partitionExtent'],) + + # Initialize data required from memory files + # Get all constants + dem = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['dem'], partition_shape) # DEM map of the study area + land_use = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['landUseMap'], partition_shape) # Land-use, example: road + soil_type = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['soilMap'], partition_shape) # example: sand or clay + resolution = float(self.configuration.modelSettings['resolution']) + self.cell_area = resolution * resolution + + # Retrieve the soil properties + self.Ks, self.porosity, self.wilting_point, self.Ki = self.retrieve_data.soil_csv( + self.configuration.generalSettings['inputDir'] + self.configuration.dataSettings['soilData'], soil_type) # soil characteristic + + # Retrieve the land-use properties + mannings, self.permeability, self.max_int_s, self.throughfall_frac = \ + self.retrieve_data.land_characteristics_csv( + self.configuration.generalSettings['inputDir'] + self.configuration.dataSettings['landUseData'], land_use, + self.cell_area) # land-use characteristics + + self.gw_base = float(self.configuration.modelSettings['groundWaterBase']) * self.standard_LUE.one() + # self.ldd = lfr.d8_flow_direction(dem) + self.ldd = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['ldd'], partition_shape) + + + #### SETTING CONSTANTS #### + ## GENERAL ## + self.start_date = utilityFunctions.string_to_datetime(self.configuration.modelSettings['startDate'], ", ") + self.stdmin = self.standard_LUE.one()* float(self.configuration.modelSettings['standard_minimal_value']) + #slope = self.standard_LUE.one() * 0.008 + #self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(slope, 0.05, 0.00001) + self.channel_width = int(self.configuration.modelSettings['channel_width']) + self.coefficient = mannings / (0.008 * self.channel_width) + aquifer_height = float(self.configuration.modelSettings['impermeableLayerBelowDEM']) + self.imperm_lay_height = dem - aquifer_height + self.water_below_dem = float(self.configuration.modelSettings['waterBelowDEM']) + self.max_gw_s = aquifer_height * self.cell_area # Full storage of porosity + self.min_gw_s = self.max_gw_s * (self.wilting_point / self.porosity) # Minimum storage because of wilting point + + # Refactorings value from mm/hour to m/s times the cell area. + self.refactor = (self.cell_area / 1000) / 3600 + self.d = (5**(2/3))/(4**(2/3)) + + # Channel length and area + self.channel_length = resolution * self.standard_LUE.one() + self.channel_area = self.channel_width * self.channel_length + + # Load initial groundWaterStorage, if no raster is supplied, use the waterBelowDEM in combination with DEM to create a initialGroundWaterStorage layer. + try: + self.gw_s = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniGroundWaterStorage'], partition_shape) + except: + print("Did not find a initial groundwater height file, looked at: {}".format(self.configuration.dataSettings['iniGroundWaterStorage'])) + self.gw_s = lfr.where(dem > (self.gw_base + self.water_below_dem), + ((dem - self.water_below_dem)-self.imperm_lay_height) * self.cell_area, + (self.gw_base - self.imperm_lay_height) * self.cell_area) + self.gw_s = lfr.where(self.gw_s > self.max_gw_s, self.max_gw_s, self.gw_s) + + # Load initial discharge, if no raster is supplied, set to zero. + try: + self.height = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniWaterHeight'], partition_shape) + except: + #print("Did not find a initial discharge file, looked at: {}".format(self.configuration.dataSettings['iniWaterHeight'])) + self.height = lfr.where(dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + try: + self.discharge = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniDischarge'], partition_shape) + except: + #print("Did not find a initial discharge file, looked at: {}".format(self.configuration.dataSettings['iniDischarge'])) + self.discharge = lfr.where(dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + # Initial InterceptionStorage and groundWaterStorage + try: + self.int_s = lfr.from_gdal(self.input_dir + self.configuration.dataSettings['iniInterceptionStorage'], partition_shape) + except: + #print("Did not find a initial interception storage file, looked at: {}".format(self.configuration.dataSettings['iniInterceptionStorage'])) + self.int_s = lfr.where(dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + + # REPORTING INITIAL CONDITIONS + variables = {"ini_gw_s": self.gw_s, "ini_int_s": self.int_s, "ini_sur_h": self.height, + } + + self.report.initial(variables) + + + def update_fluxes(self, gw_s, date, update_timestep): + # Load flux and storage values + precipitation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['precipitationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['evapotranspirationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + self.int_s, precipitation, evapotranspiration_surface = self.calculate_flux.interception(self.int_s, + self.max_int_s, + precipitation, + ref_evaporation, + self.throughfall_frac) + + evapotranspiration_surface, evapotranspiration_soil = self.calculate_flux.evapotranspiration(precipitation, + evapotranspiration_surface) + + """ direct_infiltration = self.calculate_flux.infiltration(gw_s, + self.max_gw_s, + self.Ks, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) """ + + direct_infiltration, pot_reinfiltration = self.calculate_flux.adjusted_Horton(gw_s, + self.max_gw_s, + self.Ks, + self.Ki, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) + + # Groundwater LDD, gradient and flow flux + gw_height = self.imperm_lay_height + self.gw_s / self.cell_area + gw_ldd = lfr.d8_flow_direction(gw_height) + del_h_gw = gw_height - lfr.downstream(gw_ldd, gw_height) + gw_grad = (del_h_gw) / self.channel_length + gw_flow = self.Ks * gw_grad * self.timestep * (gw_height - self.imperm_lay_height) * self.channel_length # Groundwater velocity in m2/s + + # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. + gw_flow = lfr.where(gw_flow * update_timestep > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/update_timestep, gw_flow) + gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) + gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), gw_flow) + + # Add all vertical processes for the surfacewater and all processes groundwater + gw_flux = direct_infiltration - evapotranspiration_soil + lfr.upstream(gw_ldd, gw_flow) - gw_flow + sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters + + return gw_flux, sw_flux, pot_reinfiltration + + def route_both(self, gw_s, gw_flux, discharge, sw_flux, pot_reinfiltration): + """Route surface and subsurface. + + We use the available information in regards to surface and subsurface + layers in order to determine the values of each variable in the next + state. Both layers are updated and routed. + + Args: + gw_s (LUE partitioned array): the groundwater storage in m3 of soil filled with water (m3). + gw_flux (LUE partitioned array): the groundwater flux based on input and flow (m3/s). + discharge (LUE partitioned array): the flow volume that travels through each cell (m3/s). + sw_flux (LUE partitioned array): the surface water flux (m3/s). + pot_reinfiltration (LUE partitioned array): the potential reinfiltration of runoff (m3/s). + + Returns: + height (LUE partitioned array): the water height in each cell (m). + updated_discharge (LUE partitioned array): the updated flow volume that travels through each cell (m3/s). + updated_gw_s (LUE partitioned array): the updated groundwater storage in m3 of soil filled with water. (m3) + seepage (LUE partitioned array): the amount of water leaving the soil (m3/s). + reinfiltration (LUE partitioned array): the amount of runoff that infiltrates the soil (m3/s). + """ + # Kinematic Surface Water Routing Constants + alpha = 1.5 + beta = 0.6 + c = 5/3 + + # The groundwater is adjusted by the fluxes + height = lfr.pow(self.coefficient*discharge, 0.6) + reinfiltration = lfr.where(height > pot_reinfiltration, pot_reinfiltration, height) + gw_s = gw_s + gw_flux + (reinfiltration / self.porosity) + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(gw_s > self.max_gw_s, (gw_s - self.max_gw_s)*self.porosity, 0) + + height = height - reinfiltration + discharge = lfr.pow(height, c) / self.coefficient + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + inflow = (sw_flux + seepage)/self.channel_area + inflow = lfr.where(inflow < self.stdmin, self.stdmin, inflow) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, inflow,\ + alpha, beta, self.timestep,\ + self.channel_length) + + height = lfr.pow(self.coefficient*discharge, 0.6) + + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + gw_s = gw_s - (seepage / self.porosity) + + return height, discharge, gw_s, seepage, reinfiltration + + def route_surface(self, height, sw_flux): + """Route only the surface + + The height and surface water flux are used to calculate a new + water table height, and route water in downstream direction. + + Args: + height (LUE partitioned array): the old water height (m). + sw_flux (LUE partitioned array): the surface water flux (m3/s). + + Returns: + updated_height (LUE partitioned array): the water height in each cell (m). + discharge (LUE partitioned array): the flow volume that travels through each cell (m3/s). + """ + # Kinematic Surface Water Routing Constants + alpha = 1.5 + beta = 0.6 + c = 5/3 + + # Discharge is affected by the surfacewater fluxes, and seepage is added + new_height = height + ((sw_flux)/self.channel_area) #- channel_infiltation + + discharge = lfr.pow(new_height, c) / self.coefficient + + # Because the kinematic wave has difficulties working with zero's, we have opted for a very small value. This will impact model results. + discharge = lfr.where(discharge < self.stdmin, self.stdmin, discharge) + + # Water routing based on the kinematic wave function, currently alpha is a float. Hopefully mannings raster can be used in the future. + discharge = lfr.kinematic_wave(self.ldd, discharge, self.stdmin, + alpha, beta, self.timestep, + self.channel_length) + + updated_height = lfr.pow(self.coefficient*discharge, 0.6) + + return updated_height, discharge + + def route_subsurface(self, gw_s, gw_flux): + """Route only the subsurface + + The subsurface is routed based on the supplied groundwater storage, + and the given groundwater flux. The new groundwater storage and the + seepage are returned by the function. + + Args: + gw_s (LUE partitioned array): The groundwater storage in m3 of soil filled with water. + gw_flux (LUE partitioned array): The groundwater flux based on input and flow (m3/s). + + Returns: + updated_gw_s (LUE partitioned array): The updated groundwater storage in m3 of soil filled with water. + seepage (LUE partitioned array): The amount of water leaving the soil (m3/s). + """ + # The groundwater is adjusted by the fluxes + new_gw_s = gw_s + gw_flux + + # If the groundwater table surpases the digital elevation map, groundwater is turned into runoff. + seepage = lfr.where(new_gw_s > self.max_gw_s, (new_gw_s - self.max_gw_s)*self.porosity, 0) + + updated_gw_s = new_gw_s - (seepage / self.porosity) + + return updated_gw_s, seepage + + def simulate(self, time_step): + date = self.start_date + datetime.timedelta(seconds=time_step) + + # Get new fluxes + if time_step % self.update_timestep: + self.gw_flux, self.sw_flux, self.pot_reinfiltration = self.update_fluxes(self.gw_s, date, self.update_timestep) + + # Update all variables and route. + if self.configuration.modelSettings['routing'] == "both": + self.height, self.discharge, self.gw_s, self.seepage, self.reinfiltration = self.route_both( + self.gw_s, + self.gw_flux, + self.discharge, + self.sw_flux, + self.pot_reinfiltration + ) + + variables = {"discharge": self.discharge, "int_s": self.int_s, + "gw_s": self.gw_s, "seepage": self.seepage, + "reinfiltration": self.reinfiltration + } + + elif self.configuration.modelSettings['routing'] == 'surface': + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + self.height, self.discharge = self.route_surface(self.height, + self.sw_flux + ) + + variables = {"discharge": self.discharge, "int_s": self.int_s, + "height": self.height + } + + + elif self.configuration.modelSettings['routing'] == 'subsurface': + # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly + # Then route the water lateraly and repeat this for the amount of iterations required. + self.gw_s, self.seepage = self.route_subsurface(self.gw_s, self.gw_flux) + + variables = {"gw_s": self.gw_s, "seepage": self.seepage + } + + else: + raise("Unknown 'routing' mode, change configuration setting to one of the following: both, surface, subsurface") + + # Report output + if time_step % self.report_timestep == 0: + self.report.dynamic(date, variables) + + return self.discharge.future() + + + + def finalize(self): + # Create the balance report + #self.report.balance_report(self.configuration) + + # Process the results into a gif + if self.configuration.generalSettings['makeGIF'] == 'True': + print(f"Creating a GIF for: {self.configuration.gifSettings['variables']}.") + #tools.MakeGIF.run(self.configuration) + +class MyProgressor(lfr.Progressor): + def __init__(self): + lfr.Progressor.__init__(self) + + + def initialize(self): + sys.stdout.write("[") + sys.stdout.flush() + + def simulate(self, time_step): + sys.stdout.write(".") + sys.stdout.flush() + + def finalize(self): + sys.stdout.write("]\n") + sys.stdout.flush() + +def calculate_timesteps(start_date: str, end_date: str, sep, timestep_size): + sd = utilityFunctions.string_to_datetime(start_date, sep) + ed = utilityFunctions.string_to_datetime(end_date, sep) + timesteps = int((ed - sd).total_seconds() / int(timestep_size)) + return timesteps + +@lfr.runtime_scope +def main(): + configuration = Configuration("F:/Projecten intern (2023)/Stage Steven Hosper/Model/v1/config/config_computational_1.ini") + report = Report(configuration) + model = MyModel(configuration, report) + progressor = MyProgressor() + nr_time_steps = calculate_timesteps(configuration.modelSettings['startDate'], + configuration.modelSettings['endDate'], + ", ", + configuration.modelSettings['timestep'] + ) + lfr.run_deterministic(model, progressor, nr_time_steps, rate_limit=5) + +if __name__ == "__main__": + main() + +print("--- %s seconds ---" % (time.time() - start_time)) + +if __name__ == "__main__": + main() + +print("--- %s seconds ---" % (time.time() - start_time)) + +if __name__ == "__main__": + main() + +print("--- %s seconds ---" % (time.time() - start_time)) + +if __name__ == "__main__": + main() + +print("--- %s seconds ---" % (time.time() - start_time)) + +if __name__ == "__main__": + main() + +print("--- %s seconds ---" % (time.time() - start_time)) diff --git a/model/reporting.py b/model/reporting.py index 1eb07f4..698c2ff 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -6,8 +6,8 @@ """ import lue.framework as lfr import datetime +#from osgeo import gdal import numpy as np -from osgeo import gdal import pandas as pd from StandardArraysLUE import StandardArraysLUE From 1361cd9ce6a82471befba2cc6d15f34b35656cae Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Wed, 21 Jun 2023 14:51:49 +0200 Subject: [PATCH 19/22] Update fluxes split to optimize. --- model/HBM_v2.py | 103 ++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 90 insertions(+), 13 deletions(-) diff --git a/model/HBM_v2.py b/model/HBM_v2.py index e52d98b..f3e95d4 100644 --- a/model/HBM_v2.py +++ b/model/HBM_v2.py @@ -149,7 +149,7 @@ def initialize(self): self.report.initial(variables) - def update_fluxes(self, gw_s, date, update_timestep): + def update_all_fluxes(self, gw_s, date, update_timestep): # Load flux and storage values precipitation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + self.configuration.dataSettings['precipitationData'], @@ -202,10 +202,84 @@ def update_fluxes(self, gw_s, date, update_timestep): sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters return gw_flux, sw_flux, pot_reinfiltration - - def report_output(self): - pass - + + def update_surface_fluxes(self, date): + # Load flux and storage values + gw_s = self.standard_LUE.one() * 1 + + precipitation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['precipitationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['evapotranspirationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + self.int_s, precipitation, evapotranspiration_surface = self.calculate_flux.interception(self.int_s, + self.max_int_s, + precipitation, + ref_evaporation, + self.throughfall_frac) + + evapotranspiration_surface, evapotranspiration_soil = self.calculate_flux.evapotranspiration(precipitation, + evapotranspiration_surface) + + direct_infiltration = self.calculate_flux.infiltration(gw_s, + self.max_gw_s, + self.Ks, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) + + sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters + + return sw_flux + + def update_subsurface_fluxes(self, gw_s, date, update_timestep): + # Load flux and storage values + precipitation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['precipitationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(self.configuration.generalSettings['inputDir'] + + self.configuration.dataSettings['evapotranspirationData'], + self.refactor, date) # m3/s into a cell (refactor incorporates cell area) + + self.int_s, precipitation, evapotranspiration_surface = self.calculate_flux.interception(self.int_s, + self.max_int_s, + precipitation, + ref_evaporation, + self.throughfall_frac) + + evapotranspiration_surface, evapotranspiration_soil = self.calculate_flux.evapotranspiration(precipitation, + evapotranspiration_surface) + + direct_infiltration = self.calculate_flux.infiltration(gw_s, + self.max_gw_s, + self.Ks, + self.permeability, + self.porosity, + precipitation, + evapotranspiration_surface) + + # Groundwater LDD, gradient and flow flux + gw_height = self.imperm_lay_height + self.gw_s / self.cell_area + gw_ldd = lfr.d8_flow_direction(gw_height) + del_h_gw = gw_height - lfr.downstream(gw_ldd, gw_height) + gw_grad = (del_h_gw) / self.resolution + gw_flow = self.Ks * gw_grad * self.timestep * (gw_height - self.imperm_lay_height) * self.resolution # Groundwater velocity in m2/s + + # If the groundwater flow because of the impermeable layer is larger than the amount of water available, than it should be set so only the stored water will move. + gw_flow = lfr.where(gw_flow * update_timestep > gw_s - self.min_gw_s, (gw_s - self.min_gw_s)/update_timestep, gw_flow) + gw_flow = lfr.where(gw_s < self.min_gw_s, self.stdmin, gw_flow) + gw_flow = lfr.where(gw_ldd == 5, self.standard_LUE.zero(), gw_flow) + + # Add all vertical processes for the surfacewater and all processes groundwater + gw_flux = direct_infiltration - evapotranspiration_soil + lfr.upstream(gw_ldd, gw_flow) - gw_flow + + return gw_flux + def route_both(self, gw_s, gw_flux, discharge, sw_flux, pot_reinfiltration): """Route surface and subsurface. @@ -314,12 +388,13 @@ def route_subsurface(self, gw_s, gw_flux): def simulate(self, time_step): date = self.start_date + datetime.timedelta(seconds=time_step) - # Get new fluxes - if time_step % self.update_timestep: - self.gw_flux, self.sw_flux, self.pot_reinfiltration = self.update_fluxes(self.gw_s, date, self.update_timestep) - # Update all variables and route. if self.routing == "both": + # Get new fluxes + if time_step % self.update_timestep: + self.gw_flux, self.sw_flux, self.pot_reinfiltration = self.update_all_fluxes(self.gw_s, date, self.update_timestep) + + self.height, self.discharge, self.gw_s, self.seepage, self.reinfiltration = self.route_both( self.gw_s, self.gw_flux, @@ -334,8 +409,9 @@ def simulate(self, time_step): } elif self.routing == 'surface': - # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly - # Then route the water lateraly and repeat this for the amount of iterations required. + if time_step % self.update_timestep: + self.gw_flux, self.sw_flux, self.pot_reinfiltration = self.update_surface_fluxes(date) + self.height, self.discharge = self.update_and_route_surface(self.height, self.sw_flux ) @@ -346,8 +422,9 @@ def simulate(self, time_step): elif self.routing == 'subsurface': - # Use the gw_flux and sw_flux to adjust the surface- and groundwater height accordingly - # Then route the water lateraly and repeat this for the amount of iterations required. + if time_step % self.update_timestep: + self.gw_flux, self.sw_flux, self.pot_reinfiltration = self.update_subsurface_fluxes(self.gw_s, date, self.update_timestep) + self.gw_s, self.seepage = self.update_and_route_subsurface(self.gw_s, self.gw_flux) variables = {"gw_s": self.gw_s, "seepage": self.seepage, From 0116a24d38eeb36cca638fca01c374c18d4b3b04 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Wed, 21 Jun 2023 15:02:49 +0200 Subject: [PATCH 20/22] commandline testing --- model/HBM_v2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/HBM_v2.py b/model/HBM_v2.py index f3e95d4..8882b0c 100644 --- a/model/HBM_v2.py +++ b/model/HBM_v2.py @@ -474,7 +474,7 @@ def calculate_timesteps(start_date: str, end_date: str, sep, timestep_size): @lfr.runtime_scope def main(): - configuration = Configuration("F:/Projecten intern (2023)/Stage Steven Hosper/Model/v1/config/config.ini") + configuration = Configuration(sys.argv[1]) report = Report(configuration) model = MyModel(configuration, report) progressor = MyProgressor() From 23f17d879f52de6332a7e84de83851aaf442b8ca Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Wed, 21 Jun 2023 15:07:59 +0200 Subject: [PATCH 21/22] More commandline improvements. --- model/HBM_v2.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/model/HBM_v2.py b/model/HBM_v2.py index 8882b0c..945607a 100644 --- a/model/HBM_v2.py +++ b/model/HBM_v2.py @@ -474,6 +474,8 @@ def calculate_timesteps(start_date: str, end_date: str, sep, timestep_size): @lfr.runtime_scope def main(): + if len(sys.argv) == 1: + sys.exit("Missing configuration file.") configuration = Configuration(sys.argv[1]) report = Report(configuration) model = MyModel(configuration, report) @@ -486,6 +488,7 @@ def main(): lfr.run_deterministic(model, progressor, nr_time_steps, rate_limit=15) if __name__ == "__main__": - sys.exit(main()) + main() + sys.exit(print("--- %s seconds ---" % (time.time() - start_time))) + -print("--- %s seconds ---" % (time.time() - start_time)) From 78852e6de39f2dbc7507ff7cbef1482efe6aa940 Mon Sep 17 00:00:00 2001 From: StevenHosper Date: Thu, 6 Jul 2023 11:30:19 +0200 Subject: [PATCH 22/22] finalizing --- config/config.ini | 30 ++++++++++++------------------ model/HBM.py | 17 ++++++++--------- model/HBM_v2.py | 17 ++++++++--------- model/RetrieveData.py | 12 ++++++------ model/reporting.py | 4 ++-- 5 files changed, 36 insertions(+), 44 deletions(-) diff --git a/config/config.ini b/config/config.ini index 5172360..a3ee4dc 100644 --- a/config/config.ini +++ b/config/config.ini @@ -5,34 +5,28 @@ outputDir = C:/Users/steven.hosper/Desktop/Mapje Stage/output/ network = False useAPI = False -scenario = De Hupsel +scenario = De Hupsel5 makeGIF = False +balanceCheck= False # Include Processes -includePrecipitation = True -includeEvapotranspiration = True +includePrecipitation = False +includeEvapotranspiration = False includeInfiltration = True includeInterception = True includePercolation = False -[apiSettings] -username = __key__ -password = Cy0BNm8p.vpytC2vYPT9g7OKdgxvqggyV0k9zzJVy -precipAPI = 730d6675-35dd-4a35-aa9b-bfb8155f9ca7 -evapAPI = e262dc03-f12b-4082-a4f4-7d0534e31fa4 -demAPI = a60ad336-c95b-4fb6-b852-96fc352ee808 - [modelSettings] # Date = y, m, d, h, m, s startDate = 2023, 4, 6, 12, 30, 0 -endDate = 2023, 4, 8, 12, 30, 0 +endDate = 2023, 4, 6, 13, 45, 0 iterationsBeforeData = 300 -iterationsBeforeReport = 900 +iterationsBeforeReport = 300 timestep = 1 standard_minimal_value = 1E-20 -routing = both +routing = surface waterBelowDEM = 0.0 impermeableLayerBelowDEM= 2.00 @@ -40,16 +34,18 @@ groundWaterBase = 23.25 porosity = 0.35 channel_width = 1 -arrayExtent = 5000 +arrayExtent = 1000 partitionExtent = 1000 resolution = 1 + +# Only important for the final balance calculation validCellsPercentageD = 24.7 validCellsPercentageG = 25.04 [dataSettings] iniGroundWaterStorage = /gw_test/60_gw_s_2022-04-14-1230.tiff -iniWaterHeight = /1_height_2022-04-07-1215.tiff -iniDischarge = /1_discharge_2022-04-07-1215.tiff +iniWaterHeight = +iniDischarge = iniInterceptionStorage = dem = /dem_f64_pcr_shaped.tiff ldd = /ldd_f64_pcr_shaped.tiff @@ -60,8 +56,6 @@ landUseData = /landuse_conversion_v2.csv precipitationData = /precipitation.csv evapotranspirationData = /evapotranspiration.csv -[reportSettings] -variables = discharge, seepage, groundWaterHeight, Qgw, Sgw, swFlux, gwFlux, infiltration, evapotranspirationSoil [gifSettings] variables = discharge, gw_s diff --git a/model/HBM.py b/model/HBM.py index 676faf6..569c293 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -51,7 +51,7 @@ def __init__(self, configuration, report): # Set directories self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] - self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings['scenario'] + " v2" + self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings['scenario'] partition_shape = 2 * (configuration.modelSettings['partitionExtent'],) @@ -126,7 +126,7 @@ def __init__(self, configuration, report): self.ini_sur_h = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniWaterHeight'], partition_shape) except: print("Did not find a initial discharge file, looked at: {}".format(configuration.dataSettings['iniWaterHeight'])) - self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.zero(), self.standard_LUE.zero()) + self.ini_sur_h = lfr.where(self.dem > 0, self.standard_LUE.one() * 0.001, self.standard_LUE.zero()) try: self.ini_dis = lfr.from_gdal(self.input_dir + configuration.dataSettings['iniDischarge'], partition_shape) @@ -225,7 +225,7 @@ def dis_to_height(self, discharge): height (lue partitioned array): height array that can be used to add vertical and lateral in- and outflows. """ - height = (1/self.d)*self.mannings * discharge / self.slope_sqrd + height = (((2/3) * discharge)**(1/self.beta))/ self.channel_width return height @@ -241,7 +241,7 @@ def height_to_dis(self, height): discharge (lue partitioned array): discharge array that can be used as an input for the kinematic wave formula. """ - discharge = self.d * (height * self.slope_sqrd)/self.mannings + discharge = self.alpha * ((self.channel_width * height) ** self.beta) return discharge @@ -317,7 +317,7 @@ def dynamic_model(self, configuration): # Add all vertical processes for the surfacewater and all processes groundwater gw_flux = direct_infiltration - evapotranspiration_soil + lfr.upstream(gw_ldd, gw_flow) - gw_flow - sw_flux = precipitation - evapotranspiration_surface - direct_infiltration # Is now in cubic meters + sw_flux = self.standard_LUE.zero() ############## ROUTING AND UPDATE PER TIMESTEP ############## @@ -345,12 +345,12 @@ def dynamic_model(self, configuration): height, discharge = self.update_and_route_surface(height, sw_flux) # Get the maximum value of the discharge raster (to limit the amount of tasks created by HPX) - outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)) - total_volume = lfr.minimum(lfr.zonal_sum(height, self.dem>0)) + outflow = lfr.minimum(lfr.zonal_sum(discharge, self.ldd == 5)).get() + total_volume = lfr.maximum(lfr.zonal_sum(height, height>0)).get() print("outflow: ", outflow, "total: ", total_volume) # Write value to csv for later validation - writer.writerow([i*60 + j, outflow]) + writer.writerow([i*60 + j, outflow, total_volume]) variables = {"discharge": discharge, "int_s": int_s, "height": height } @@ -383,7 +383,6 @@ def dynamic_model(self, configuration): # Save / Report data if i % (dt_report/dt_data) == 0: - print(i) report.dynamic(date, variables) diff --git a/model/HBM_v2.py b/model/HBM_v2.py index 945607a..76f386e 100644 --- a/model/HBM_v2.py +++ b/model/HBM_v2.py @@ -48,7 +48,7 @@ def __init__(self, configuration, report): self.standard_LUE = StandardArraysLUE(configuration) self.retrieve_data = RetrieveData(configuration) self.calculate_flux = CalculateFlux(configuration) - self.report = report + self.report = report def initialize(self): # Set directories @@ -88,7 +88,7 @@ def initialize(self): ## GENERAL ## self.routing = self.configuration.modelSettings['routing'] self.start_date = utilityFunctions.string_to_datetime(self.configuration.modelSettings['startDate'], ", ") - self.stdmin = float(self.configuration.modelSettings['standard_minimal_value']) + self.stdmin = float(self.configuration.modelSettings['standard_minimal_value']) * self.standard_LUE.one() slope = self.standard_LUE.one() * 0.008 self.slope_sqrd = utilityFunctions.calculate_sqrd_slope(slope, 0.05, 0.00001) self.channel_width = int(self.configuration.modelSettings['channel_width']) @@ -387,7 +387,7 @@ def route_subsurface(self, gw_s, gw_flux): def simulate(self, time_step): date = self.start_date + datetime.timedelta(seconds=time_step) - + # Update all variables and route. if self.routing == "both": # Get new fluxes @@ -410,9 +410,9 @@ def simulate(self, time_step): elif self.routing == 'surface': if time_step % self.update_timestep: - self.gw_flux, self.sw_flux, self.pot_reinfiltration = self.update_surface_fluxes(date) + self.sw_flux = self.update_surface_fluxes(date) - self.height, self.discharge = self.update_and_route_surface(self.height, + self.height, self.discharge = self.route_surface(self.height, self.sw_flux ) @@ -474,9 +474,7 @@ def calculate_timesteps(start_date: str, end_date: str, sep, timestep_size): @lfr.runtime_scope def main(): - if len(sys.argv) == 1: - sys.exit("Missing configuration file.") - configuration = Configuration(sys.argv[1]) + configuration = Configuration("F:/Projecten intern (2023)/Stage Steven Hosper/Model/v1/config/config.ini") report = Report(configuration) model = MyModel(configuration, report) progressor = MyProgressor() @@ -485,7 +483,8 @@ def main(): ", ", configuration.modelSettings['timestep'] ) - lfr.run_deterministic(model, progressor, nr_time_steps, rate_limit=15) + + lfr.run_deterministic(model, progressor, nr_time_steps, rate_limit=5) if __name__ == "__main__": main() diff --git a/model/RetrieveData.py b/model/RetrieveData.py index d5cee26..9db5594 100644 --- a/model/RetrieveData.py +++ b/model/RetrieveData.py @@ -87,7 +87,7 @@ def land_characteristics_csv(self, data_file, land_use, cell_area): # Standard values dummy = self.std_arr_lue.one() permeability = dummy * 0.8 - mannings = dummy * 0.045 + mannings = dummy * 0.03 interception_storage_max = dummy * 0.001 throughfall_fraction = dummy * 0.90 LAI = dummy * 0.5 @@ -104,12 +104,12 @@ def land_characteristics_csv(self, data_file, land_use, cell_area): crop_factor_value = data["Crop_type"] # When an ID matches the array ID, assign value - for count, ID in enumerate(ID): - mannings = lfr.where(land_use == ID, mannings_friction[count], mannings) - permeability = lfr.where(land_use == ID, permeability_value[count], permeability) - interception_storage_max = lfr.where(land_use == ID, interception_storage_max_value[count], interception_storage_max) + #for count, ID in enumerate(ID): + # mannings = lfr.where(land_use == ID, mannings_friction[count], mannings) + # permeability = lfr.where(land_use == ID, permeability_value[count], permeability) + # interception_storage_max = lfr.where(land_use == ID, interception_storage_max_value[count], interception_storage_max) # LAI = lfr.where(land_use == ID, LAI_value[count], LAI) - throughfall_fraction = lfr.where(land_use == ID, throughfall_fraction_value[count], throughfall_fraction) + # throughfall_fraction = lfr.where(land_use == ID, throughfall_fraction_value[count], throughfall_fraction) # crop_factor = lfr.where(land_use == ID, crop_factor_value[count], crop_factor) interception_storage_max = interception_storage_max * cell_area diff --git a/model/reporting.py b/model/reporting.py index 698c2ff..27e2837 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -6,7 +6,7 @@ """ import lue.framework as lfr import datetime -#from osgeo import gdal +from osgeo import gdal import numpy as np import pandas as pd from StandardArraysLUE import StandardArraysLUE @@ -17,7 +17,7 @@ def __init__(self, configuration): self.standard_LUE = StandardArraysLUE(configuration) self.timestep = int(configuration.modelSettings['timestep']) self.iterations = int(configuration.modelSettings['iterationsBeforeReport']) - self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings["scenario"] + " v2" + self.output_dir = configuration.generalSettings['outputDir'] + configuration.generalSettings["scenario"] self.input_dir = configuration.generalSettings['inputDir'] + configuration.generalSettings['scenario'] def initial(self, variables):