using generic LddStructureChan and IsUpOfStructureChan

Author: Cinzia Mazzetti

src/lisflood/Lisflood_dynamic.py

@@ -191,7 +191,7 @@ def splitlanduse(array1, array2=None, array3=None):
         # if option['simulatePolders']:
         # ChannelToPolderM3=ChannelToPolderM3Old;
 
-        # # cmcheck moved to routing.py
+        # # moved to routing.py
         # # Calculate total water storage in river channel
         # if option['InitLisflood'] or (not(option['SplitRouting'])):
         #     # kinematic routing
@@ -224,13 +224,13 @@ def splitlanduse(array1, array2=None, array3=None):
 
         if option['InitLisflood'] or option['repAverageDis']:
             self.CumQ += self.ChanQ
-            #cmcheck we should use ChanQAvg here not ChanQ
+            #cmcheck - we should use ChanQAvg here not ChanQ
             self.avgdis = self.CumQ/self.TimeSinceStart
             # to calculate average discharge over the entire simulation
 
         self.DischargeM3Out += np.where(self.AtLastPointC ,self.ChanQ * self.DtSec,0)
         # Cumulative outflow out of map
-        # cmcheck we should use ChanQAvg here not ChanQ
+        # cmcheck - we should use ChanQAvg here not ChanQ
 
         # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         # Calculate water level

src/lisflood/hydrological_modules/evapowater.py

@@ -60,8 +60,8 @@ def initial(self):
         maskinfo = MaskInfo.instance()
         if option['openwaterevapo']:
             LakeMask = loadmap('LakeMask', pcr=True)
-            # cmcheck
-            lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresKinematic, 5)
+            # lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresKinematic, 5)
+            lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresChan, 5)
             LddEva = lddrepair(lmask)
             lddC = compressArray(LddEva)
             inAr = decompress(np.arange(maskinfo.info.mapC[0], dtype="int32"))
@@ -133,8 +133,6 @@ def dynamic(self):
             UpstreamEva = self.var.EWRef * self.var.MMtoM3 * self.var.WaterFraction
             # evaporation for loop is amount of water per timestep [cu m]
             # Volume of potential evaporation from water surface  per time step (conversion to [m3])
-
-            # cmcheck
             # ChanMIter = self.var.ChanM3Kin.copy()
             ChanMIter = self.var.ChanM3.copy()
             # for Iteration loop: First value is amount of water in the channel

src/lisflood/hydrological_modules/indicatorcalc.py

@@ -132,7 +132,8 @@ def dynamic(self):
                 # EXTERNAL FLOW
                 wreg = decompress(self.var.WUseRegionC)
                      # to pcraster map because of the following expression!
-                self.var.RegionMonthExternalInflowM3 = compressArray(areatotal(cover(ifthen(self.var.WaterRegionInflowPoints != 0, upstream(self.var.LddStructuresKinematic,decompress(self.var.MonthDisM3))),0),wreg))
+                # self.var.RegionMonthExternalInflowM3 = compressArray(areatotal(cover(ifthen(self.var.WaterRegionInflowPoints != 0, upstream(self.var.LddStructuresKinematic,decompress(self.var.MonthDisM3))),0),wreg))
+                self.var.RegionMonthExternalInflowM3 = compressArray(areatotal(cover(ifthen(self.var.WaterRegionInflowPoints != 0, upstream(self.var.LddStructuresChan,decompress(self.var.MonthDisM3))),0),wreg))
                 self.var.RegionMonthAbstractionRequiredAllSourcesM3 =  np.take(np.bincount(self.var.WUseRegionC,weights=self.var.MonthAbstractionRequiredAllSourcesM3),self.var.WUseRegionC)
                 self.var.RegionMonthAbstractionRequiredSurfaceGroundWaterM3 = np.take(np.bincount(self.var.WUseRegionC,weights=self.var.MonthAbstractionRequiredSurfaceGroundWaterM3),self.var.WUseRegionC)
                 self.var.RegionMonthAbstractionRequiredSurfaceWaterM3 = np.take(np.bincount(self.var.WUseRegionC,weights=self.var.MonthAbstractionRequiredSurfaceWaterM3),self.var.WUseRegionC)

src/lisflood/hydrological_modules/inflow.py

@@ -125,7 +125,7 @@ def dynamic(self):
             # Get inflow hydrograph at each inflow point [m3/s]
             QIn = compressArray(QIn)
             QIn[np.isnan(QIn)] = 0
-            #cmcheck
+            #cmcheck - inflow
             self.var.QInM3 = QIn * self.var.DtSec
             # Convert to [m3] per time step
             self.var.TotalQInM3 += self.var.QInM3
@@ -145,7 +145,7 @@ def dynamic_inloop(self, NoRoutingExecuted):
         settings = LisSettings.instance()
         option = settings.options
 
-        if option['inflow']: #cmcheck
+        if option['inflow']: #cmcheck - inflow
 
             self.var.QInDt = (self.var.QInM3Old + (NoRoutingExecuted + 1) * self.var.QDelta) * self.var.InvNoRoutSteps
             # flow from inlets per sub step

src/lisflood/hydrological_modules/routing.py

@@ -109,19 +109,20 @@ def initial(self):
         # upstream of gauge locations
         # Calculate inverse, so we can multiply in dynamic (faster than divide)
 
-        self.var.IsChannelPcr = boolean(loadmap('Channels', pcr=True))  #pcr
+        self.var.IsChannelPcr = boolean(loadmap('Channels', pcr=True))  #pcr map
         self.var.IsChannel = np.bool8(compressArray(self.var.IsChannelPcr))     #bool
         self.var.IsChannelPcr = boolean(decompress(self.var.IsChannel)) #pcr
-
         # Identify grid cells containing a river channel
+
+        self.var.IsChannelKinematic = self.var.IsChannel.copy()     #bool
+        # Identify channel pixels using kinematic or split routing
+        # (identical to IsChannel, unless MCT wave is used, see below)
+
         self.var.IsStructureChan = np.bool8(maskinfo.in_zero())        #bool
         # Initialise map that identifies special inflow/outflow structures (reservoirs, lakes) within the
         # channel routing. Set to (dummy) value of zero modified in reservoir and lake
         # routines (if those are used)
 
-        self.var.IsChannelKinematic = self.var.IsChannel.copy()     #bool
-        # Identify kinematic wave channel pixels
-        # (identical to IsChannel, unless MCT wave is used, see below)
         self.var.IsStructureKinematic = np.bool8(maskinfo.in_zero())        #bool
         # Initialise map that identifies special inflow/outflow structures (reservoirs, lakes) within the
         # kinematic wave channel routing. Set to (dummy) value of zero modified in reservoir and lake
@@ -231,7 +232,7 @@ def initial(self):
         # Area (sq m) of bank full discharge cross section [m2]
         # (trapezoid area equation)
 
-        #cmcheck
+        #cmcheck - TotalCrossSectionAreaHalfBankFull is not 1/2 TotalCrossSectionAreaBankFull it's trapezoid
         # ChanUpperWidthHalfBankFull = self.var.ChanBottomWidth + 2 * self.var.ChanSdXdY * 0.5 * ChanDepthThreshold
         # TotalCrossSectionAreaHalfBankFull = 0.5 * \
         #     0.5 * ChanDepthThreshold * (ChanUpperWidthHalfBankFull + self.var.ChanBottomWidth)
@@ -394,7 +395,7 @@ def initial(self):
         # Cumulative inflow volume from inflow hydrographs [m3]
         self.var.sumDis = maskinfo.in_zero()
         self.var.sumIn = maskinfo.in_zero()
-        # cmcheck non so se sostituita da self.var.sumInWB
+        # cmcheck - non so se sostituita da self.var.sumInWB
         self.var.sumInWB = maskinfo.in_zero()
 
 
@@ -470,7 +471,6 @@ def initialSecond(self):
                                            self.var.Beta, self.var.ChanLength, self.var.DtRouting,
                                           alpha_floodplains=self.var.ChannelAlpha2, flagnancheck=flags['nancheck'])
 
-        # cmcheck
         # ************************************************************
         # ***** WATER BALANCE                             ************
         # ************************************************************
@@ -488,15 +488,15 @@ def initialSecond(self):
 
         if not option['InitLisflood'] and option['repMBTs']:
            self.var.StorageStepINIT = self.var.ChanM3
-           DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
-           # cmchek to fix IsUpsOfStructureKinematicC for MCT
+           # DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+           DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
            if not(option['SplitRouting']):
             # self.var.StorageStepINIT = self.var.ChanM3Kin
             # self.var.StorageStepINIT = self.var.ChanM3
             if option['simulateReservoirs']:
                self.var.StorageStepINIT += self.var.ReservoirStorageIniM3
-               DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
-               # cmchek to fix IsUpsOfStructureKinematicC for MCT
+               # DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+               DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
             if option['simulateLakes']:
                self.var.StorageStepINIT += self.var.LakeStorageIniM3
                DisStructure += np.where(compressArray(self.var.IsUpsOfStructureLake), 0.5 * self.var.ChanQ * self.var.DtRouting, 0)
@@ -815,18 +815,17 @@ def dynamic(self, NoRoutingExecuted):
                        if option['simulateReservoirs']:
                          sum1 =self.var.ChanQ.copy()
                          StorageStep =  StorageStep + self.var.ReservoirStorageM3.copy()
-                         # cmcheck IsUpsOfStructureKinematicC
-                         DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
-                         # cmchek to fix IsUpsOfStructureKinematicC for MCT
+                         # DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+                         DisStructureR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
                          DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureR), self.var.Catchments)
                          DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
 
                       if not option['InitLisflood']:
                        if option['simulateLakes']:
                          sum1 =self.var.ChanQ.copy()
                          StorageStep =  StorageStep + self.var.LakeStorageM3Balance.copy()
-                         DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
-                         # cmchek to fix IsUpsOfStructureKinematicC for MCT
+                         # DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+                         DisStructureR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
                          DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureR), self.var.Catchments)
                          maskinfo = MaskInfo.instance()
                          DisLake = maskinfo.in_zero()
@@ -864,16 +863,16 @@ def dynamic(self, NoRoutingExecuted):
                          sum1=[]
                          sum1 =self.var.ChanQ.copy()
                          StorageStep =  StorageStep + self.var.ReservoirStorageM3.copy()
-                         DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
-                         # cmchek to fix IsUpsOfStructureKinematicC for MCT
+                         # DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+                         DisStructureSR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
                          DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureSR), self.var.Catchments)
                          DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
 
                       if option['simulateLakes']:
                          sum1 =self.var.ChanQ.copy()
                          StorageStep =  StorageStep + self.var.LakeStorageM3Balance.copy()
-                         DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
-                         # cmchek to fix IsUpsOfStructureKinematicC for MCT
+                         # DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+                         DisStructureSR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
                          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)

src/lisflood/hydrological_modules/structures.py

@@ -43,33 +43,41 @@ def __init__(self, structures_variable):
     def initial(self):
         """ initial part of the structures module
         """
-        self.var.LddStructuresKinematic = self.var.LddKinematic #pcr
+        self.var.LddStructuresKinematic = self.var.LddKinematic     #pcr map
         LddStructuresKinematicNp = compressArray(self.var.LddStructuresKinematic)
-        self.var.LddStructuresChan = self.var.LddChan   #pcr
+        # Unmodified version of LddKinematic is needed to connect inflow and outflow points
+        # of each structure (called LddStructuresKinematic now)
+
+        self.var.LddStructuresChan = self.var.LddChan   #pcr map
         LddStructuresChanNp = compressArray(self.var.LddStructuresChan)
+        # Unmodified version of LddChan is needed to connect inflow and outflow points
+        # of each structure (called LddStructuresChan now)
 
         settings = LisSettings.instance()
         option = settings.options
         if not option['InitLisflood']:
             # not done in Init Lisflood
-            IsUpsOfStructureKinematic = downstream(
+            IsUpsOfStructureKinematic = downstream(     #pcr map
                 self.var.LddKinematic,
                 cover(boolean(decompress(self.var.IsStructureKinematic)), boolean(0))
             )
-            IsUpsOfStructureChan = downstream(
+            # Downstream assigns to result the expression value of the neighbouring downstream cell
+            # Over is used to cover missing values on an expression with values taken from one or more different expression(s)
+            # Decompress is numpy2pcr
+            # Find location of pixels immediately upstream of a structure on the LddKinematic
+
+            IsUpsOfStructureChan = downstream(      #pcr map
                 self.var.LddChan,
                 cover(boolean(decompress(self.var.IsStructureChan)), boolean(0))
             )
+            # Find location of pixels immediately upstream of a structure on the LddChan
+
+            self.var.IsUpsOfStructureKinematicC = compressArray(IsUpsOfStructureKinematic)  #np compressed array
+            # Location of pixels immediately upstream of a structure on the LddKinematic
+            self.var.IsUpsOfStructureChanC = compressArray(IsUpsOfStructureChan)    #np compressed array
+            # Location of pixels immediately upstream of a structure on the LddChan
 
-            # Get all pixels just upstream of kinematic structure locations
-            self.var.IsUpsOfStructureKinematicC = compressArray(IsUpsOfStructureKinematic)  #np
-            # Unmodified version of LddKinematic is needed to connect inflow and outflow points
-            # of each structure (called LddStructuresKinematic now)
-            self.var.IsUpsOfStructureChanC = compressArray(IsUpsOfStructureChan)    #np
-            # Unmodified version of LddChan is needed to connect inflow and outflow points
-            # of each structure (called LddStructuresChan now)
-
-            self.var.LddKinematic = lddrepair(ifthenelse(IsUpsOfStructureKinematic, 5, self.var.LddKinematic))
-            # Cells just upstream of each structure are treated as pits in the kinematic wave channel routing
-            self.var.LddChan = lddrepair(ifthenelse(IsUpsOfStructureChan, 5, self.var.LddChan))
-            # Cells just upstream of each structure are treated as pits in the kinematic wave channel routing
+            self.var.LddKinematic = lddrepair(ifthenelse(IsUpsOfStructureKinematic, 5, self.var.LddKinematic))  #pcr map
+            # Update LddKinematic by adding a pit in the pixel immediately upstream of a structure
+            self.var.LddChan = lddrepair(ifthenelse(IsUpsOfStructureChan, 5, self.var.LddChan))     #pcr map
+            # Update LddChan by adding a pit in the pixel immediately upstream of a structure

src/lisflood/hydrological_modules/waterabstraction.py

@@ -159,14 +159,16 @@ def initial(self):
                 pitWuse2 = ifthen(pitWuseMax == self.var.UpArea, WUseRegion)
                 # search outlets in the inland water regions by using the maximum  upstream area as criterium
 
-                pitWuse3 = downstream(self.var.LddStructuresKinematic, WUseRegion)
+                # pitWuse3 = downstream(self.var.LddStructuresKinematic, WUseRegion)
+                pitWuse3 = downstream(self.var.LddStructuresChan, WUseRegion)
                 pitWuse3b = ifthen(pitWuse3 != WUseRegion, WUseRegion)
                 # search point where ldd leaves a water region
 
                 pitWuse = cover(pitWuse1b, pitWuse2, pitWuse3b, nominal(0))
                 # join all sources of pits
 
-                LddWaterRegion = lddrepair(ifthenelse(pitWuse == 0, self.var.LddStructuresKinematic, 5))
+                # LddWaterRegion = lddrepair(ifthenelse(pitWuse == 0, self.var.LddStructuresKinematic, 5))
+                LddWaterRegion = lddrepair(ifthenelse(pitWuse == 0, self.var.LddStructuresChan, 5))
                 # create a Ldd with pits at every water region outlet
                 # this results in a interrupted ldd, so water cannot be transfered to the next water region
                 lddC = compressArray(LddWaterRegion)
@@ -185,8 +187,10 @@ def initial(self):
                 # outflowpoints to calculate upstream inflow for balances and Water Exploitation Index
                 # both inland outflowpoints to downstream subbasin, and coastal outlets
 
+                # WaterRegionInflow1 = boolean(
+                #     upstream(self.var.LddStructuresKinematic, cover(scalar(self.var.WaterRegionOutflowPoints), 0)))
                 WaterRegionInflow1 = boolean(
-                    upstream(self.var.LddStructuresKinematic, cover(scalar(self.var.WaterRegionOutflowPoints), 0)))
+                    upstream(self.var.LddStructuresChan, cover(scalar(self.var.WaterRegionOutflowPoints), 0)))
                 self.var.WaterRegionInflowPoints = ifthen(WaterRegionInflow1, boolean(1))
                 # inflowpoints to calculate upstream inflow for balances and Water Exploitation Index
             else:

src/lisflood/hydrological_modules/waterbalance.py

@@ -89,7 +89,8 @@ def initial(self):
             #DisStructure = cover(ifthen(
             #    self.var.IsUpsOfStructureKinematic, self.var.ChanQ * self.var.DtRouting), scalar(0.0))
 
-            DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+            # DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+            DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
 
             # Discharge upstream of structure locations (coded as pits) in [m3/time step]
             # Needed for mass balance error calculations (see comment to calculation of WaterInit below)
@@ -232,7 +233,8 @@ def dynamic(self):
             # for this double counting)
             # new 12.11.09 PB
             # added cumulative transmission loss
-            DisStru = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+            # DisStru = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+            DisStru = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
             DischargeM3Structures = np.take(np.bincount(self.var.Catchments, weights=DisStru), self.var.Catchments)
 
             # on the last time step lakes and reservoirs calculated with the previous routing results