Merge pull request #114 from ec-jrc/release_proposal_4.1.0

Release 4.1.0

README.md

@@ -15,51 +15,6 @@ Other useful resources
 | Lisflood Usecases   |                                                           | https://github.com/ec-jrc/lisflood-usecases                  |
 
 
-## Notes for Release 4.0.0
-
-LISFLOOD-OS v4.0.0 includes the following changes:
-
-1) pixel-by-pixel computation of water infiltration into the soil within the soilloop.py module.
-
-2) changes in the allocation of water abstraction.
-
-3) patches in the modules riceirrigation.py and frost.py to avoid non-realistic results in specific conditions.
-
-4) improvements in the computation of the mass balance.
-
-5) bug fixes in the use of the transientlandusechange option, in the use of the option simulatePF, in the writing of the monthly outputs.
-
-Moreover, to ensure compatibility with LISFLOOD-EPIC, modules source code is compatible with XArray variables, and XArray variables will be allocated when EPIC module is ON.  
-LISFLOOD-EPIC has been developed by JRC-D2. When compared to LISFLOOD v3.2.0, LISFLOOD-EPIC enables a more detailed representation of crop growth and irrigation. References:  ASR - Assessing groundwater irrigation sustainability in the Euro-Mediterranean region with an integrated agro-hydrologic model (copernicus.org) , Gelati, E. et al.: Integrated model of irrigation, crop growth and water resources (LISFLOOD-EPIC): description and application in the Euro-Mediterranean area, in preparation, 2020.  
-Important note: EPIC modules are not *yet* available from the OS-LISFLOOD repository.  
-
-1) <ins>computation of water infiltration into the soil within the soilloop.py module.</ins>  
-LISFLOOD v3.2.0: all the pixels of a land use fraction are allocated to a vector and the number of loops is defined by the most critical pixel. More specifically, all the pixels within the computational domain are allocated the number of iterations that are necessary to allow the numerical stability of the most critical pixel. The number of iterations consequently depends on the specific computational domain. This can lead to different results when modelling sub-catchmennts and entire basins with highly heterogenous soils.  
-LISFLOOD v4.0.0: pixel by pixel computation and parallelization of the computations by using the python package numba.  
-Benefits of v4.0.0: each pixel is allocated the number of iterations that allow its numerical stability, this guarantees the consistency of results when modelling sub-catchments and the entire basin with highly heterogeneous soils. Decreased computational time due to parallel computations.  
-
-2) <ins>changes in the allocation of water abstraction.</ins>  
-LISFLOOD v3.2.0: direct abstraction of the consumptive water use.  
-LISFLOOD v4.0.0: abstraction of the demanded water volume for each use (industrial, domestic, energy, livestock); the consumptive water use leaves the system; the unused water volume is returned to the channels.  
-Minor note:  small differences in the computation of the leakages.  
-
-3) <ins>patches in the modules riceirrigation.py and frost.py to avoid non-realistic results in specific conditions.</ins>  
-riceirrigation.py: LISFLOOD v3.2.0: in presence of a very thick third soil layer (10^2 m) the drainage of the rice fields causes non realistic large flow discharge values in the channels.  
-LISFLOOD v4.0.0: only soil layers 1a and 1b are used by the rice module (for both the saturation and the drainage phases). By definition, the second soil layer includes the roots depth.  
-frost.py: LISFLOOD v3.2.0: there is no upper boundary to the frost index value. Rain events on deeply frozen soils are totally converted into runoff and can cause false alarms.  
-LISFLOOD v4.0.0: the maximum frost index value is set to 57.  
-
-4) <ins>improvements in the computation of the mass balance.</ins>  
-LISFLOOD v3.2.0: the mass balance is computed by considering the input and output volumes up to the current computational step.  Small numerical errors accumulate over time and give the erroneous impression of large mass balance errors.
-LISFLOOD v4.0.0: computation of the mass balance by considering the single computational step.
-
-5) <ins>bug fixes in the use of the transientlandusechange option, in the use of the option simulatePF, in the writing of the monthly outputs.</ins>  
-Correct use of the rice fraction by landusechange.py DYNAMIC; use of the correct soil water content variable for the computation of pF1 (ws1a replaced ws1); correct reporting of the monthly time steps by the optional module indicatorcalc.py.  
-
-Lisflood 4.0.0 is the version used for EFAS 5 and GloFAS 4 calibration.  
-
-Note: the users are recommended to download the [reference settings xml](https://github.com/ec-jrc/lisflood-code/tree/master/src/lisfloodSettings_reference.xml) file and adapt it by inserting their own paths and modelling choices).
-
 
 ## Quick start
 
@@ -185,11 +140,15 @@ Command above will also install the executable `lisflood` in the virtualenv, so
 lisflood /absolute_path/to/my/local/folder/LF_ETRS89_UseCase/settings/cold.xml
 ```
 
-### IMPORTANT NOTE
+### IMPORTANT NOTE 1
 
 Please note that there are known issues when installing the LISFLOOD code on Windows (source code and pypi package). We cannot provide Windows support and strongly recommend using LINUX to install the LISFLOOD code.
 Windows users are recommended to execute LISFLOOD with a Docker image.
 
+### IMPORTANT NOTE 2
+
+The users are recommended to download the [reference settings xml](https://github.com/ec-jrc/lisflood-code/tree/master/src/lisfloodSettings_reference.xml) file and adapt it by inserting their own paths and modelling choices.
+
 
 
 ## Collaborate

VERSION

@@ -1 +1 @@
-4.0.1
+4.1.0

requirements.txt

@@ -5,8 +5,8 @@ coverage>=6.0
 Cython>=0.29.14
 dask>=2021.10.0 
 future>=0.18.0
-GDAL>=3.3.0
-importlib-metadata>=0.23
+GDAL>=3.2.0
+importlib-metadata<5.0.0
 iniconfig>=1.1.1
 lisflood-utilities>=0.12.19
 llvmlite>=0.37.0

src/lisflood/global_modules/default_options.py

@@ -702,7 +702,7 @@
                                                          restrictoption=['nonInit'], monthly=False,
                                                          yearly=False),
                        'PrefFlowMaps': ReportedMap(name='PrefFlowMaps', output_var='PrefFlowPixel',
-                                                   unit='mm', end=[], steps=[], all=['repE2O2,repPrefFlowMaps'],
+                                                   unit='mm', end=[], steps=[], all=['repE2O2','repPrefFlowMaps'],
                                                    restrictoption=['nonInit'], monthly=False,
                                                    yearly=False),
                        'RainMaps': ReportedMap(name='RainMaps', output_var='Rain', unit='mm',

src/lisflood/hydrological_modules/lakes.py

@@ -57,7 +57,7 @@ def initial(self):
         binding = settings.binding
         maskinfo = MaskInfo.instance()
 
-        if option['simulateLakes']:
+        if option['simulateLakes'] and not option['InitLisflood']:
 
             LakeSitesC = loadmap('LakeSites')
             LakeSitesC[LakeSitesC < 1] = 0
@@ -99,6 +99,8 @@ def initial(self):
             LakeArea = pcraster.lookupscalar(str(binding['TabLakeArea']), LakeSitePcr)
             LakeAreaC = compressArray(LakeArea)
             self.var.LakeAreaCC = np.compress(LakeSitesC > 0, LakeAreaC)
+
+            self.var.LakeSitesC2 = LakeSitesC   # additional var 
 
             # Surface area of each lake [m2]
             LakeA = pcraster.lookupscalar(str(binding['TabLakeA']), LakeSitePcr)
@@ -176,7 +178,8 @@ def initial(self):
 
             self.var.LakeStorageM3CC = LakeStorageIniM3CC.copy()
             self.var.LakeStorageM3BalanceCC = LakeStorageIniM3CC.copy()
-
+
+
             self.var.LakeStorageIniM3 = maskinfo.in_zero()
             self.var.LakeLevel = maskinfo.in_zero()
             self.var.LakeInflowOld = maskinfo.in_zero()
@@ -186,6 +189,7 @@ def initial(self):
             np.put(self.var.LakeLevel, self.var.LakeIndex, self.var.LakeLevelCC)
             np.put(self.var.LakeInflowOld, self.var.LakeIndex, self.var.LakeInflowOldCC)
             np.put(self.var.LakeOutflow, self.var.LakeIndex, self.var.LakeOutflowCC)
+
 
             self.var.EWLakeCUMM3 = maskinfo.in_zero()
             self.var.EWLakeWBM3 = maskinfo.in_zero()
@@ -204,7 +208,10 @@ def dynamic_inloop(self, NoRoutingExecuted):
         option = settings.options
         maskinfo = MaskInfo.instance()
         if not(option['InitLisflood']) and option['simulateLakes']:    # only with no InitLisflood
-            #self.var.LakeInflow = cover(ifthen(defined(self.var.LakeSites), upstream(self.var.LddStructuresKinematic, self.var.ChanQ)), scalar(0.0))
+
+            if NoRoutingExecuted==0:
+                self.var.LakeStorageM3CC=np.compress(self.var.LakeSitesC2 > 0, self.var.LakeStorageM3)
+
             self.var.LakeInflowCC = np.bincount(self.var.downstruct, weights=self.var.ChanQ)[self.var.LakeIndex]
             # Lake inflow in [m3/s]
 
@@ -269,6 +276,7 @@ def dynamic_inloop(self, NoRoutingExecuted):
                     self.var.sumLakeOutCC += QLakeOutM3DtCC
                     # summing up over all sub timesteps
 
+
             if NoRoutingExecuted == (self.var.NoRoutSteps-1):
 
                 # expanding the size after last sub timestep
@@ -283,6 +291,7 @@ def dynamic_inloop(self, NoRoutingExecuted):
                 np.put(self.var.LakeInflowOld, self.var.LakeIndex, self.var.LakeInflowOldCC)
                 np.put(self.var.LakeOutflow, self.var.LakeIndex, self.var.LakeOutflowCC)
 
+
                 if option['repsimulateLakes']:
                     np.put(self.var.LakeInflowM3S, self.var.LakeIndex, self.var.sumLakeInCC / self.var.DtSec)
                     np.put(self.var.LakeOutflowM3S, self.var.LakeIndex, self.var.sumLakeOutCC / self.var.DtSec)

src/lisflood/hydrological_modules/reservoir.py

@@ -58,7 +58,7 @@ def initial(self):
         settings = LisSettings.instance()
         option = settings.options
         maskinfo = MaskInfo.instance()
-        if option['simulateReservoirs']:
+        if option['simulateReservoirs'] and not option['InitLisflood']:
 
             # NoSubStepsRes=max(1,roundup(self.var.DtSec/loadmap('DtSecReservoirs')))
             # Number of sub-steps based on value of DtSecReservoirs,
@@ -181,7 +181,7 @@ def dynamic_inloop(self, NoRoutingExecuted):
         settings = LisSettings.instance()
         option = settings.options
         maskinfo = MaskInfo.instance()
-        if option['simulateReservoirs']:
+        if option['simulateReservoirs'] and not option['InitLisflood']:
             InvDtSecDay = 1 / float(86400)
             # InvDtSecDay=self.var.InvDtSec
             # ReservoirInflow = cover(ifthen(defined(self.var.ReservoirSites), upstream(
@@ -197,7 +197,13 @@ def dynamic_inloop(self, NoRoutingExecuted):
 
             QResInM3Dt = ReservoirInflowCC * self.var.DtRouting
             # Reservoir inflow in [m3] per timestep (routing step)
-
+
+            # print('RESERVOIRS MODULE IN')
+            # print(np.sum(self.var.ReservoirStorageM3))
+            if NoRoutingExecuted==0:
+                self.var.ReservoirStorageM3CC=np.compress(self.var.ReservoirSitesC > 0, self.var.ReservoirStorageM3)  ########################
+            # print(np.sum(self.var.ReservoirStorageM3CC))
+
             self.var.ReservoirStorageM3CC += QResInM3Dt
             # New reservoir storage [m3] = plus inflow for this sub step
             self.var.ReservoirFillCC = self.var.ReservoirStorageM3CC / self.var.TotalReservoirStorageM3CC
@@ -296,6 +302,8 @@ def dynamic_inloop(self, NoRoutingExecuted):
                     self.var.sumResOutCC += QResOutM3DtCC
                     # summing up over all sub timesteps
 
+
+
             if NoRoutingExecuted == (self.var.NoRoutSteps-1):
 
                 # expanding the size after last sub timestep
@@ -304,6 +312,11 @@ def dynamic_inloop(self, NoRoutingExecuted):
                 np.put(self.var.ReservoirStorageM3, self.var.ReservoirIndex, self.var.ReservoirStorageM3CC)
                 np.put(self.var.ReservoirFill, self.var.ReservoirIndex, self.var.ReservoirFillCC)
 
+                # print('RESERVOIRS MODULE OUT')
+                # print(np.sum(self.var.ReservoirStorageM3))
+                # print(np.sum(self.var.ReservoirStorageM3CC))
+
+
                 if option['repsimulateReservoirs']:
                     np.put(self.var.ReservoirInflowM3S, self.var.ReservoirIndex, self.var.sumResInCC / self.var.DtSec)
                     np.put(self.var.ReservoirOutflowM3S, self.var.ReservoirIndex, self.var.sumResOutCC / self.var.DtSec)

src/lisflood/hydrological_modules/routing.py

@@ -389,6 +389,8 @@ def initialSecond(self):
 
                 self.var.Chan2M3Kin = self.var.CrossSection2Area * self.var.ChanLength + self.var.Chan2M3Start
                 self.var.ChanM3Kin = self.var.ChanM3 - self.var.Chan2M3Kin + self.var.Chan2M3Start
+
+                self.var.ChanM3Kin = np.where((self.var.ChanM3Kin < 0.0) & (self.var.ChanM3Kin > -0.0000001),0.0,self.var.ChanM3Kin)
 
                 self.var.Chan2QKin = (self.var.Chan2M3Kin * self.var.InvChanLength * self.var.InvChannelAlpha2) ** (self.var.InvBeta)
                 self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta)

src/lisflood/hydrological_modules/waterabstraction.py

@@ -73,7 +73,11 @@ def initial(self):
             self.var.FractionGroundwaterUsed = np.minimum(
                 np.maximum(loadmap('FractionGroundwaterUsed'), maskinfo.in_zero()), 1.0)
             self.var.FractionNonConventionalWaterUsed = loadmap('FractionNonConventionalWaterUsed')
-            self.var.FractionLakeReservoirWaterUsed = loadmap('FractionLakeReservoirWaterUsed')
+
+            if not option['InitLisflood']:
+               self.var.FractionLakeReservoirWaterUsed = loadmap('FractionLakeReservoirWaterUsed')
+            else:
+               self.var.FractionLakeReservoirWaterUsed = maskinfo.in_zero()  
             self.var.EFlowThreshold = loadmap('EFlowThreshold')
             # EFlowThreshold is map with m3/s discharge, e.g. the 10th percentile discharge of the baseline run
 
@@ -431,22 +435,22 @@ def dynamic(self):
             # 9. Lakes and reservoirs abstraction
             # ************************************************************
             # 9.1 Max abstractable volumes from lakes and reservoirs: pixel values and integral over water use regions 
-            if option['simulateReservoirs']:
+            if option['simulateReservoirs'] and not(option['InitLisflood']):
                 PotentialAbstractionFromReservoirsM3 = np.minimum(0.02 * self.var.ReservoirStorageM3,
                                                                   0.01 * self.var.TotalReservoirStorageM3C) * self.var.DtDay  #  self.var.DtDay is required when DtSec is not 86400
                 PotentialAbstractionFromReservoirsM3 = np.where(np.isnan(PotentialAbstractionFromReservoirsM3), 0,
                                                                 PotentialAbstractionFromReservoirsM3)
             else:
                 PotentialAbstractionFromReservoirsM3 = maskinfo.in_zero()
 
-            if option['simulateLakes']:
+            if option['simulateLakes'] and not(option['InitLisflood']):
                 PotentialAbstractionFromLakesM3 = 0.10 * self.var.LakeStorageM3 * self.var.DtDay #  self.var.DtDay is required when DtSec is not 86400
                 PotentialAbstractionFromLakesM3 = np.where(np.isnan(PotentialAbstractionFromLakesM3), 0,
                                                            PotentialAbstractionFromLakesM3)
             else:
                 PotentialAbstractionFromLakesM3 = maskinfo.in_zero()
 
-            if option['simulateReservoirs'] or option['simulateLakes']:
+            if option['simulateReservoirs'] or option['simulateLakes'] and not(option['InitLisflood']):
                 PotentialAbstractionFromLakesAndReservoirsM3 = PotentialAbstractionFromLakesM3 + PotentialAbstractionFromReservoirsM3
                 # potential total m3 that can be extracted from all lakes and reservoirs in a pixel
             else:
@@ -459,30 +463,23 @@ def dynamic(self):
             # 9.2 Water regions' required and actual abstraction from lakes (Lak) and reservoirs (Res)
             areatotal_withdrawal_LakRes_required_M3 = self.var.FractionLakeReservoirWaterUsed * areatotal_withdrawal_SW_required
 
+            #print('this areatotal_withdrawal_LakRes_required_M3')
+            #print(np.sum(areatotal_withdrawal_LakRes_required_M3))
             # total amount of m3 abstracted from all lakes and reservoirs in the water regions
             self.var.areatotal_withdrawal_LakRes_actual_M3 = np.minimum(areatotal_withdrawal_LakRes_required_M3, AreatotalPotentialAbstractionFromLakesAndReservoirsM3)
+            #print(' self.var.areatotal_withdrawal_LakRes_actual_M3')
+            #print(np.sum( self.var.areatotal_withdrawal_LakRes_actual_M3))
             FractionAbstractedByLakesReservoirs = np.where(is_SW_withdrawal_required_WUR, self.var.areatotal_withdrawal_LakRes_actual_M3 / areatotal_withdrawal_SW_required,maskinfo.in_zero())
             # 9.3 Distribute actual abstractions among lakes and reservoirs inside each water region
             FractionLakesReservoirsEmptying = np.where(AreatotalPotentialAbstractionFromLakesAndReservoirsM3 > 0,
                                                    self.var.areatotal_withdrawal_LakRes_actual_M3 / AreatotalPotentialAbstractionFromLakesAndReservoirsM3, maskinfo.in_zero())
             # 9.4 Update the storages of lakes and reservoirs
             self.var.LakeAbstractionM3 = PotentialAbstractionFromLakesM3 * FractionLakesReservoirsEmptying
-            if option['simulateLakes']:
+            if option['simulateLakes'] and not(option['InitLisflood']):
                 self.var.LakeStorageM3 = self.var.LakeStorageM3 - self.var.LakeAbstractionM3             
-                if option['InitLisflood']:
-                   self.var.LakeAreaCCindex = maskinfo.in_zero()
-                   np.put(self.var.LakeAreaCCindex,self.var.LakeIndex,self.var.LakeAreaCC)
-                   self.var.LakeLevelindex = np.where(self.var.LakeStorageM3>0, self.var.LakeStorageM3 / self.var.LakeAreaCCindex, 0.0)
-                   self.var.LakeLevel = self.var.LakeLevelindex
-
             self.var.ReservoirAbstractionM3 = PotentialAbstractionFromReservoirsM3 * FractionLakesReservoirsEmptying
-            if option['simulateReservoirs']:          
+            if option['simulateReservoirs'] and not(option['InitLisflood']):          
                 self.var.ReservoirStorageM3 = self.var.ReservoirStorageM3 - self.var.ReservoirAbstractionM3
-                if option['InitLisflood']:
-                   self.var.TotalReservoirStorageM3CCindex = maskinfo.in_zero()
-                   np.put(self.var.TotalReservoirStorageM3CCindex,self.var.ReservoirIndex,self.var.TotalReservoirStorageM3CC)
-                   self.var.ReservoirFillindex =np.where(self.var.ReservoirStorageM3>0, self.var.ReservoirStorageM3 / self.var.TotalReservoirStorageM3CCindex,0.0)                   
-                   self.var.ReservoirFill = self.var.ReservoirFillindex
             # subtract abstracted water from lakes and reservoir storage
 
 

src/lisfloodSettings_reference.xml

@@ -1711,6 +1711,9 @@ This is necessary when using projected coordinates (x,y)
 </comment>
 </textvar>
 
+<textvar name="NetCDFTimeChunks" value="$(NetCDFTimeChunks)"/>
+<textvar name="MapsCaching" value="$(MapsCaching)"/>
+
 <textvar name="MaskMap" value="$(MaskMap)">
 <comment>
 Clone map