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.ini b/config/config.ini index 4f97352..a3ee4dc 100644 --- a/config/config.ini +++ b/config/config.ini @@ -8,6 +8,7 @@ useAPI = False scenario = De Hupsel5 makeGIF = False +balanceCheck= False # Include Processes includePrecipitation = False @@ -16,50 +17,49 @@ 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, 7, 13, 30, 0 -iterationsBeforeReport = 60 +endDate = 2023, 4, 6, 13, 45, 0 +iterationsBeforeData = 300 +iterationsBeforeReport = 300 timestep = 1 +standard_minimal_value = 1E-20 +routing = surface waterBelowDEM = 0.0 impermeableLayerBelowDEM= 2.00 groundWaterBase = 23.25 porosity = 0.35 +channel_width = 1 arrayExtent = 1000 partitionExtent = 1000 -resolution = 5 -validCellsPercentage = 35.62 +resolution = 1 + +# Only important for the final balance calculation +validCellsPercentageD = 24.7 +validCellsPercentageG = 25.04 [dataSettings] -iniGroundWaterStorage = +iniGroundWaterStorage = /gw_test/60_gw_s_2022-04-14-1230.tiff iniWaterHeight = +iniDischarge = 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 [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/config/config_computational_1.ini b/config/config_computational_1.ini new file mode 100644 index 0000000..a01921e --- /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, 30, 50 +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 = 1000 +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..7fe04d9 --- /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 = 1000 +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/model/CalculateFlux.py b/model/CalculateFlux.py index 5cbd834..1cc2920 100644 --- a/model/CalculateFlux.py +++ b/model/CalculateFlux.py @@ -16,14 +16,16 @@ 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 = self.std_arr_lue.zero() + self.k = 1.6 def interception(self, int_stor, int_stor_max, pre, rev, th_f): @@ -35,16 +37,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 +54,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 +66,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 +99,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 +114,57 @@ 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, 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 + + # 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)) + + 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) + + + + pot_reinfiltration = pot_infiltration - infiltration2 + + return infiltration2, pot_reinfiltration \ No newline at end of file diff --git a/model/HBM.py b/model/HBM.py index 942583c..569c293 100644 --- a/model/HBM.py +++ b/model/HBM.py @@ -19,9 +19,10 @@ from StandardArraysLUE import StandardArraysLUE from RetrieveData import RetrieveData from CalculateFlux import CalculateFlux +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() @@ -40,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'] @@ -59,30 +61,56 @@ def __init__(self, configuration): 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( + 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) # 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) self.ldd = lfr.from_gdal(self.input_dir + configuration.dataSettings['ldd'], partition_shape) - # Set constants - self.resolution = float(configuration.modelSettings['resolution']) - self.cell_area = self.resolution * self.resolution - self.slope = lfr.slope(self.dem, self.resolution) + + + #### 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 + 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) @@ -95,62 +123,142 @@ 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.one() * 0.001, 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 = 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): - pass + 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) - @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]) - dT = int((end_date - start_date).seconds / dt) + height = height - reinfiltration + discharge = lfr.pow(height, self.c) / self.coefficient - # 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 + # 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) - # 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) + height = lfr.pow(self.coefficient*discharge, 0.6) - # 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 + # Any water that is moved from groundwater to discharge has to be removed from the groundwaterStorage + gw_s = gw_s - (seepage / self.porosity) - # Refactorings value from mm/hour to m/h times the cell area. - refactor = (self.cell_area / 1000) / 3600 + 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 - # Kinematic Surface Water Routing Constants - alpha = 1.5 - beta = 0.6 - timestep = 1.0 * float(configuration.modelSettings['timestep']) - c = 5/3 + 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) - # Static, really small value because inflow = 0 is not accepted - inflow = self.standard_LUE.one()*1E-20 + # 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 = (((2/3) * discharge)**(1/self.beta))/ self.channel_width + + 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.alpha * ((self.channel_width * height) ** self.beta) + + 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: @@ -159,16 +267,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_data*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) # m3/s into a cell (refactor incorporates cell area) ref_evaporation = self.retrieve_data.csv_timeseries_to_flux(configuration.generalSettings['inputDir'] + configuration.dataSettings['evapotranspirationData'], - 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, @@ -179,72 +287,106 @@ def dynamic_model(self, configuration, report): 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) + """ direct_infiltration = 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 * channel_rat # is in m/s + 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 * 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) - gw_flow = lfr.where(gw_s < self.min_gw_s, 1E-20, 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(), 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 + gw_flux = direct_infiltration - evapotranspiration_soil + lfr.upstream(gw_ldd, gw_flow) - gw_flow + sw_flux = self.standard_LUE.zero() - 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) + + ############## 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)).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, total_volume]) - # 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) + variables = {"discharge": discharge, "int_s": int_s, "height": height + } + - # Write value to csv for later validation - writer.writerow([i*60 + j, outflow]) + 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: + report.dynamic(date, variables) + + print(f"Done: {i+1}/{dT}") - variables = {"discharge": discharge, "int_s": int_s, "height": height, "gw_s": gw_s, - } - report.dynamic(date, variables) return 0 @@ -274,13 +416,13 @@ 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 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 new file mode 100644 index 0000000..76f386e --- /dev/null +++ b/model/HBM_v2.py @@ -0,0 +1,493 @@ +# -*- 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): + # 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']) * 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']) + 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_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'], + 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 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. + + 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) + + # 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, + 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': + if time_step % self.update_timestep: + self.sw_flux = self.update_surface_fluxes(date) + + 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.routing == 'subsurface': + 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, + "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=5) + +if __name__ == "__main__": + main() + sys.exit(print("--- %s seconds ---" % (time.time() - start_time))) + + 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/RetrieveData.py b/model/RetrieveData.py index 6b577d7..9db5594 100644 --- a/model/RetrieveData.py +++ b/model/RetrieveData.py @@ -40,9 +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) @@ -50,24 +52,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 - 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 @@ -81,7 +87,7 @@ def land_characteristics_csv(self, data_file, land_use): # 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 @@ -98,14 +104,16 @@ def land_characteristics_csv(self, data_file, land_use): 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 + return mannings, permeability, interception_storage_max, throughfall_fraction def rounddown_datetime(self, dt): @@ -146,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() diff --git a/model/reporting.py b/model/reporting.py index 8542009..27e2837 100644 --- a/model/reporting.py +++ b/model/reporting.py @@ -6,16 +6,25 @@ """ import lue.framework as lfr import datetime -import numpy as np from osgeo import gdal +import numpy as np import pandas as pd +from StandardArraysLUE import StandardArraysLUE # Reporting for the HydrologicBaseModel class Report: def __init__(self, configuration): - self.timestep = configuration.modelSettings['timestep'] + 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.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") @@ -24,14 +33,13 @@ def dynamic(self, date, variables): variable, dateTime )) - return 0 def balance_report(self, configuration): start_date = self.string_to_datetime(configuration.modelSettings['startDate'], seperator= ", ") 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 +49,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"]) @@ -51,31 +59,49 @@ 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: + 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_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']) - 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"))) + 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: + 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: + 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_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") @@ -98,13 +124,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): @@ -119,14 +142,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 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.') 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