Merge pull request #186 from ec-jrc/feature/initialization

Feature/initialization (includes MCT)

Author: doc78

VERSION

@@ -1 +1 @@
-4.3.1.9998
+4.3.1.500

docs/3_step5_model-initialisation/index.md

@@ -38,3 +38,4 @@ tomli>=1.2.1
 toolz>=0.10.0
 xarray>=0.20.2
 zipp>=0.6.0
+scipy

requirements.txt

@@ -228,12 +228,14 @@ def splitlanduse(array1, array2=None, array3=None):
             WaterLevelDyn = -9999
             # Set water level dynamic wave to dummy value (needed
 
-        if option['InitLisflood'] or option['repAverageDis']:
-            # self.CumQ += self.ChanQ
+        if option['InitLisflood']:
+            if  (self.TimeSinceStart > np.round(self.NumDaysSpinUp/self.DtDay)) :
+                self.CumQ += self.ChanQAvg
+                self.avgdis = self.CumQ/(self.TimeSinceStart + self.TimeSinceStartPrerunChunkInit[0] - np.round(self.NumDaysSpinUp/self.DtDay))
+                # to calculate average discharge over the entire simulation
+        elif option['repAverageDis']:
             self.CumQ += self.ChanQAvg
-            #cmcheck - we should use ChanQAvg here not ChanQ
-            self.avgdis = self.CumQ/self.TimeSinceStart
-            # to calculate average discharge over the entire simulation
+            self.avgdis = self.CumQ/(self.TimeSinceStart)
 
         #self.DischargeM3Out += np.where(self.AtLastPointC ,self.ChanQ * self.DtSec,0)
         self.DischargeM3Out += np.where(self.AtLastPointC, self.ChanQAvg * self.DtSec, 0)

src/lisflood/Lisflood_dynamic.py

@@ -10,8 +10,8 @@
 
 __version__ = version
 __authors__ = "Ad de Roo, Emiliano Gelati, Peter Burek, Johan van der Knijff, Niko Wanders"
-__date__ = "23/01/2024"
-__copyright__ = "Copyright 2019-2024, European Commission - Joint Research Centre"
-__maintainer__ = "Stefania Grimaldi, Cinzia Mazzetti, Carlo Russo, Damien Decremer, Corentin Carton De Wiart, Valerio Lorini, Ad de Roo"
+__date__ = "27/02/2025"
+__copyright__ = "Copyright 2019-2025, European Commission - Joint Research Centre"
+__maintainer__ = "Stefania Grimaldi, Cinzia Mazzetti, Jesus Casado Rodriguez, Carlo Russo, Damien Decremer, Corentin Carton De Wiart"
 __status__ = "Operation"
 __institution__ = "European Commission DG Joint Research Centre (JRC) - E1, D2 Units"

src/lisflood/__init__.py

@@ -78,7 +78,7 @@ def __call__(self, *args, **kwargs):
                 return_data = self.cache[key]
         else:
             return_data = self.cache[key]
-            self.found[self.name] += 1
+            self.found[self.name] += 1          
         return return_data
 
     @classmethod

src/lisflood/global_modules/decorators.py

@@ -17,6 +17,7 @@
 """
 from __future__ import (absolute_import, print_function, unicode_literals)
 
+
 from future.backports import OrderedDict
 from future.utils import with_metaclass
 from nine import (iteritems, str, range, map, nine)
@@ -42,6 +43,7 @@
 from .errors import LisfloodError, LisfloodWarning, LisfloodFileError
 from .decorators import cached
 from .default_options import default_options
+##from .add1 import loadmap # Carlo, not sure....
 
 
 project_dir = os.path.normpath(os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../..'))
@@ -394,6 +396,12 @@ def __init__(self, settings_file, sys_args=[]):
         self.step_start_dt = inttodate(self.step_start_int - 1, ref_date_start, binding=self.binding)
         self.step_end_dt = inttodate(self.step_end_int - 1, ref_date_start, binding=self.binding)
         self.maskpath = self.binding['MaskMap']
+        self.NumDaysSpinUp_d = self.binding['NumDaysSpinUp'] 
+        self.NumDaysSpinUp_d = int(self.NumDaysSpinUp_d)
+        DtSec_d1 = self.binding['DtSec']
+        DtSec_d = int(DtSec_d1)
+        self.DtDay_d = DtSec_d / 86400.
+
 
     def build_reportedmaps_dicts(self):
         self.report_timeseries = self._report_tss()
@@ -448,7 +456,7 @@ def _check_simulation_dates(self):
 
         int_start, str_start = datetoint(self.binding['StepStart'], self.binding)
         int_end, str_end = datetoint(self.binding['StepEnd'], self.binding)
-
+        self.numsteps = int_end
         # test if start and end > begin
         if (int_start < 0) or (int_end < 0) or ((int_end - int_start) < 0):
             str_begin = begin.strftime("%d/%m/%Y %H:%M")

src/lisflood/global_modules/default_options.py

@@ -95,6 +95,27 @@ def initial(self):
         self.var.LZ = np.where(LZInitValue == -9999, LZSteady, LZInitValue)
         # Initialise lower store with steady-state value
         # if LZInitValue is set to -9999
+        
+        if option['InitLisflood']:
+            self.var.LZInflowCUM = loadmap('LZInflowCUMInit') 
+            self.var.TimeSinceStartPrerunChunkInit = loadmap('TimeSinceStartPrerunChunkInit')        
+            if isinstance(self.var.TimeSinceStartPrerunChunkInit, float):
+               MAP = []
+               MAP = maskinfo.in_zero() +   self.var.TimeSinceStartPrerunChunkInit
+               self.var.TimeSinceStartPrerunChunkInit= []
+               self.var.TimeSinceStartPrerunChunkInit = MAP       
+            # Cumulative lower zone inflow [mm]
+            # Needed for calculation of average LZ inflow (initialisation)
+            # Added the option to read LZInflowCUM from settings file: this is needed when the prerun is completed in separate chunks (e.g. 1990-200 and 2000-2010).
+        else:
+           self.var.NumDaysSpinUp = 0.0 # this allows to write lzavin in the model run, recommended ONLY when the run is done in 1 chunk
+           self.var.TimeSinceStartPrerunChunkInit = 0.0 # this allows to write lzavin in the model run, recommended ONLY when the run is done in 1 chunk
+           if isinstance(self.var.TimeSinceStartPrerunChunkInit, float):
+               MAP = []
+               MAP = maskinfo.in_zero() +   self.var.TimeSinceStartPrerunChunkInit
+               self.var.TimeSinceStartPrerunChunkInit= []
+               self.var.TimeSinceStartPrerunChunkInit = MAP              
+           self.var.LZInflowCUM = maskinfo.in_zero() # this allows to write lzavin in the model run, recommended ONLY when the run is done in 1 chunk        
 
         # Water in lower store [mm]
         self.var.LZThreshold = loadmap('LZThreshold')
@@ -119,9 +140,6 @@ def initial(self):
 
         self.var.GwLossCUM = maskinfo.in_zero()
         # Cumulative groundwater loss [mm]
-        self.var.LZInflowCUM = maskinfo.in_zero()
-        # Cumulative lower zone inflow [mm]
-        # Needed for calculation of average LZ inflow (initialisation)
 
         self.var.GwPercUZLZ = self.var.allocateVariableAllVegetation()
         self.var.GwLossLZ = maskinfo.in_zero()
@@ -135,8 +153,6 @@ def dynamic(self):
         # outflow from LZ to channel stops when LZ is below its threshold.
         # LZ can be below its threshold because of water abstractions
         self.var.LZOutflow = np.minimum(self.var.LowerZoneK * self.var.LZ, self.var.LZ - self.var.LZThreshold)
-        # Total UZ 
-        self.var.UZtotal = self.var.deffraction(self.var.UZ)
         # Outflow out of lower zone [mm per model timestep]
         self.var.LZOutflow = np.maximum(self.var.LZOutflow, 0)
         # No outflow if LZ is below threshold
@@ -163,19 +179,22 @@ def dynamic(self):
         self.var.LZ = self.var.LZ - self.var.GwLossLZ
         # (ground)water in lower response box [mm]
 
-        self.var.LZInflowCUM += (self.var.GwPercUZLZPixel - self.var.GwLossLZ)
-        # cumulative inflow into lower zone (can be used to improve
-        self.var.LZInflowCUM = np.maximum(self.var.LZInflowCUM, 0.0)
-        # LZInflowCUM would become negativ, if LZInit is set to high
-        # therefore this line is preventing LZInflowCUM getting negativ
-
         self.var.GwLossPixel = self.var.GwLossLZ
         # from GROUNDWATER TRANSFER to here
         # Compute pixel-average fluxes
         self.var.GwLossCUM += self.var.GwLossPixel
         self.var.GwLossWB = self.var.GwLossPixel
         # Accumulated groundwater loss over simulation period [mm]
-        self.var.LZAvInflow = (self.var.LZInflowCUM * self.var.InvDtDay) / self.var.TimeSinceStart
-        # Average inflow into lower zone over executed time steps [mm/day]
+        
+        if  (self.var.TimeSinceStart > np.round(self.var.NumDaysSpinUp/self.var.DtDay)) : 
+            self.var.LZInflowCUM += (self.var.GwPercUZLZPixel - self.var.GwLossLZ)
+            # cumulative inflow into lower zone (can be used to improve
+            self.var.LZInflowCUM = np.maximum(self.var.LZInflowCUM, 0.0)
+            # LZInflowCUM would become negativ, if LZInit is set to high
+            # therefore this line is preventing LZInflowCUM getting negativ
+            self.var.LZAvInflow = (self.var.LZInflowCUM * self.var.InvDtDay) / ( self.var.TimeSinceStart + self.var.TimeSinceStartPrerunChunkInit[0] - np.round(self.var.NumDaysSpinUp/self.var.DtDay))
+            # Average inflow into lower zone over executed time steps [mm/day]
+            self.var.TimeSinceStartPrerunChunk = self.var.TimeSinceStartPrerunChunkInit + self.var.TimeSinceStart - np.round(self.var.NumDaysSpinUp/self.var.DtDay)
+       
 
         self.var.LZOutflowToChannelPixel = self.var.LZOutflowToChannel

src/lisflood/global_modules/settings.py

@@ -57,9 +57,8 @@ def initial(self):
         binding = settings.binding
         maskinfo = MaskInfo.instance()
 
-        if option['simulateLakes'] and not option['InitLisflood']:
-
-            LakeSitesC = loadmap('LakeSites')
+        if option['simulateLakes']:
+            LakeSitesC = loadmap('LakeSites')               # moved here to use the caching feature during calibration
             LakeSitesC[LakeSitesC < 1] = 0
             LakeSitesC[self.var.IsChannel == 0] = 0
             # Get rid of any lakes that are not part of the channel network
@@ -75,7 +74,10 @@ def initial(self):
                 # rebuild lists of reported files with repsimulateLakes and simulateLakes = False
                 settings.build_reportedmaps_dicts()
                 return
-            # break if no lakes
+                # break if no lakes
+            LakeSitePcr = loadmap('LakeSites', pcr=True)    # moved here to use the caching feature during calibration
+
+        if option['simulateLakes'] and not option['InitLisflood']:            
 
             self.var.IsStructureKinematic = np.where(LakeSitesC > 0, np.bool8(1), self.var.IsStructureKinematic)
             # Add lake locations to structures map (used to modify LddKinematic
@@ -86,7 +88,6 @@ def initial(self):
 
             # PCRaster part
             # -----------------------
-            LakeSitePcr = loadmap('LakeSites', pcr=True)
             LakeSitePcr = pcraster.ifthen((pcraster.defined(LakeSitePcr) & pcraster.boolean(decompress(self.var.IsChannel))), LakeSitePcr)
             IsStructureLake = pcraster.boolean(LakeSitePcr)
             # additional structure map only for lakes to calculate water balance

src/lisflood/hydrological_modules/groundwater.py

@@ -102,6 +102,7 @@ def initial(self):
         self.var.DtSecChannel = loadmap('DtSecChannel')
         # Sub time step used for kinematic wave channel routing [seconds]
         # within the model,the smallest out of DtSecChannel and DtSec is used
+        self.var.NumDaysSpinUp = loadmap('NumDaysSpinUp')
 
         self.var.MMtoM = 0.001
         # Multiplier to convert wate depths in mm to meters

src/lisflood/hydrological_modules/lakes.py

@@ -36,7 +36,7 @@ def __init__(self, readmeteo_variable):
         # initialise xarray readers
         if option['readNetcdfStack']:
             self.forcings = {}
-            for data in ['PrecipitationMaps', 'TavgMaps', 'ET0Maps', 'E0Maps']:
+            for data in ['PrecipitationMaps', 'TavgMaps', 'ET0Maps', 'ES0Maps', 'E0Maps']:
                 self.forcings[data] = xarray_reader(data)
 
 # --------------------------------------------------------------------------
@@ -61,6 +61,7 @@ def dynamic(self):
             self.var.Precipitation = self.forcings['PrecipitationMaps'][step] * self.var.DtDay * self.var.PrScaling
             self.var.Tavg = self.forcings['TavgMaps'][step]
             self.var.ETRef = self.forcings['ET0Maps'][step] * self.var.DtDay * self.var.CalEvaporation
+            self.var.ESRef = self.forcings['ES0Maps'][step] * self.var.DtDay * self.var.CalEvaporation
             self.var.EWRef =self.forcings['E0Maps'][step] * self.var.DtDay * self.var.CalEvaporation
 
         else:
@@ -71,11 +72,11 @@ def dynamic(self):
             # average DAILY temperature (even if you are running the model on say an hourly time step) [degrees C]
             self.var.ETRef = readmapsparse(binding['ET0Maps'], self.var.currentTimeStep(), self.var.ETRef) * self.var.DtDay * self.var.CalEvaporation
             # daily reference evapotranspiration (conversion to [mm] per time step)
+            self.var.ESRef = readmapsparse(binding['ES0Maps'], self.var.currentTimeStep(), self.var.ESRef) * self.var.DtDay * self.var.CalEvaporation
             # potential evaporation rate from a bare soil surface (conversion to [mm] per time step)
             self.var.EWRef = readmapsparse(binding['E0Maps'], self.var.currentTimeStep(), self.var.EWRef) * self.var.DtDay * self.var.CalEvaporation
             # potential evaporation rate from water surface (conversion to [mm] per time step)
 
-        self.var.ESRef = (self.var.EWRef + self.var.ETRef)/2
 
         if option['TemperatureInKelvin']:
             self.var.Tavg -= 273.15

src/lisflood/hydrological_modules/miscInitial.py

@@ -90,13 +90,11 @@ def initial(self):
 
         settings = LisSettings.instance()
         option = settings.options
+        binding = settings.binding
         maskinfo = MaskInfo.instance()
-        if option['simulateReservoirs'] and not option['InitLisflood']:
-
-            binding = settings.binding
-
+        if option['simulateReservoirs']:
             # load reservoir locations and keep only those on the channel network
-            reservoirs = loadmap('ReservoirSites')
+            reservoirs = loadmap('ReservoirSites')                  # moved here to use the caching feature during calibration
             reservoirs[(reservoirs < 1) | (self.var.IsChannel == 0)] = 0
             self.var.ReservoirSitesC = reservoirs
             self.var.ReservoirSitesCC = np.compress(reservoirs > 0, reservoirs)
@@ -110,6 +108,9 @@ def initial(self):
                 # rebuild lists of reported files with simulateReservoirs and repsimulateReservoirs = False
                 settings.build_reportedmaps_dicts()
                 return
+            ReservoirSitePcr = loadmap('ReservoirSites', pcr=True)    # moved here to use the caching feature during calibration
+
+        if option['simulateReservoirs'] and not option['InitLisflood']:
 
             # Add reservoir locations to structures map 
             # (used to modify LddKinematic and to calculate LddStructuresKinematic)
@@ -119,11 +120,9 @@ def initial(self):
             self.var.IsStructureChan = np.where(self.var.ReservoirSitesC > 0, np.bool8(1), self.var.IsStructureChan)
             # Add reservoir locations to structures map (used to modify LddChan
             # and to calculate LddStructuresChan)
-
-            ReservoirSitePcr = loadmap('ReservoirSites', pcr=True)
+            
             self.var.ReservoirSites = ReservoirSitePcr
-            ReservoirSitePcr = ifthen((defined(ReservoirSitePcr) & boolean(decompress(self.var.IsChannel))), ReservoirSitePcr)
-            # Get rid of any reservoirs that are not part of the channel network
+            # filter out reservoirs that are not part of the channel network
             # (following logic of 'old' code the inflow into these reservoirs is
             # always zero, so either change this or leave them out!)
             ReservoirSitePcr = ifthen((defined(ReservoirSitePcr) & boolean(decompress(self.var.IsChannel))), ReservoirSitePcr)

src/lisflood/hydrological_modules/readmeteo.py

@@ -310,6 +310,8 @@ def initial(self):
         # Initialise discharge at kinematic wave pixels (note that InvBeta is simply 1/beta, computational efficiency!)
 
         self.var.CumQ = maskinfo.in_zero()
+        if option['InitLisflood']:
+           self.var.CumQ = loadmap('CumQInit')
         # Initialise sum of discharge to calculate average
 
         # ************************************************************
@@ -399,13 +401,12 @@ def initial(self):
             # For pixels in order 1 and beyond, upstream contribution is calculated during the calculation step.
             # Initialisation is needed when using Lakes because lakes use average discharge from previous step to calculate the inflow.
 
-            if option['simulateLakes'] or option['simulateReservoirs'] and not option['InitLisflood']:
-                # Initialising average discharge for lakes
-                PrevDischargeAvg = loadmap('PrevDischargeAvg')
-                # Outflow (x+dx) Q during previous routing sub-step for full cross-section (average over last routing sub-step)
-                # Used to calculated average Inflow (x) to reservoirs and lakes
-                self.var.ChanQAvgDt = np.where(PrevDischargeAvg == -9999, self.var.ChanQAvgDt, PrevDischargeAvg)  # np
-                self.var.ChanQKinAvgDt = np.where(PrevDischargeAvg == -9999, self.var.ChanQKinAvgDt, PrevDischargeAvg)  # np
+            # Initialising average discharge for lakes, reservoirs and transmission loss warm start
+            PrevDischargeAvg = loadmap('PrevDischargeAvg')
+            # Outflow (x+dx) Q during previous routing sub-step for full cross-section (average over last routing sub-step)
+            # Used to calculated average Inflow (x) to reservoirs and lakes
+            self.var.ChanQAvgDt = np.where(PrevDischargeAvg == -9999, self.var.ChanQAvgDt, PrevDischargeAvg)  # np
+            self.var.ChanQKinAvgDt = np.where(PrevDischargeAvg == -9999, self.var.ChanQKinAvgDt, PrevDischargeAvg)  # np
 
         # ************************************************************
         # ***** CUMULATIVE OUTPUT VARIABLES  *************************
@@ -635,7 +636,11 @@ def dynamic(self, NoRoutingExecuted):
             self.polder_module.dynamic_inloop()
 
         self.inflow_module.dynamic_inloop(NoRoutingExecuted)
-        self.transmission_module.dynamic_inloop()
+        self.transmission_module.dynamic_inloop(NoRoutingExecuted)
+
+        # ************************************************************
+        # ***** CHANNEL FLOW ROUTING: KINEMATIC WAVE  ****************
+        # ************************************************************
 
         if not(option['dynamicWave']):
 
@@ -882,6 +887,7 @@ def dynamic(self, NoRoutingExecuted):
                          StorageStep =  StorageStep + self.var.ReservoirStorageM3.copy()
                          # DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
                          DisStructureSR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
+                         DisStructureSR[self.var.AtLastPointC == 1 ] = 0 # this line avoids double-counting when a reservoir is located at the outlet of the cacthment
                          DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureSR), self.var.Catchments)
                          DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
 
@@ -890,6 +896,7 @@ def dynamic(self, NoRoutingExecuted):
                          StorageStep =  StorageStep + self.var.LakeStorageM3Balance.copy()
                          # DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
                          DisStructureSR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
+                         DisStructureSR[self.var.AtLastPointC == 1 ] = 0 # this line avoids double-counting when a reservoir is located at the outlet of the cacthment
                          DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureSR), self.var.Catchments)
                          DisLake = maskinfo.in_zero()
                          np.put(DisLake, self.var.LakeIndex, 0.5 * self.var.LakeInflowCC * self.var.DtRouting)