Feature/mct calpoint

src/lisflood/Lisflood_dynamic.py

@@ -209,7 +209,6 @@ def splitlanduse(array1, array2=None, array3=None):
         #     # Total channel storage [m3] = Volume in main channel (ChanM3Kin) + volume above bankfull (Chan2M3Kin - Chan2M3Start)
         #     # at t+dt
 
-
         self.TotalCrossSectionArea = self.ChanM3 * self.InvChanLength
         # Total river channel cross-section area at t+dt
 
@@ -240,7 +239,6 @@ def splitlanduse(array1, array2=None, array3=None):
         #self.DischargeM3Out += np.where(self.AtLastPointC ,self.ChanQ * self.DtSec,0)
         self.DischargeM3Out += np.where(self.AtLastPointC, self.ChanQAvg * self.DtSec, 0)
         # Cumulative outflow out of map
-        # cmcheck - we should use ChanQAvg here not ChanQ
 
         # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         # Calculate water level

src/lisflood/Lisflood_initial.py

@@ -49,6 +49,7 @@
 from .hydrological_modules.surface_routing import surface_routing
 from .hydrological_modules.reservoir import Reservoir
 from .hydrological_modules.lakes import lakes
+from .hydrological_modules.mctconfluence import mctconfluence
 from .hydrological_modules.polder import polder
 from .hydrological_modules.waterabstraction import waterabstraction
 from .hydrological_modules.indicatorcalc import indicatorcalc
@@ -138,6 +139,7 @@ def __init__(self):
         self.surface_routing_module = surface_routing(self)
         self.reservoir_module = Reservoir(self) # get_reservoir(option['reservoirHanazaki'])
         self.lakes_module = lakes(self)
+        self.mctconfluence_module = mctconfluence(self)
         self.polder_module = polder(self)
         self.waterabstraction_module = waterabstraction(self)
         self.indicatorcalc_module = indicatorcalc(self)
@@ -189,10 +191,12 @@ def __init__(self):
 
         self.snow_module.initial()
         self.frost_module.initial()
+
         self.leafarea_module.initial()
-        self.soilloop_module.initial()
 
+        self.soilloop_module.initial()
         self.soil_module.initial()
+
         self.routing_module.initial()
 
         self.groundwater_module.initial()
@@ -201,6 +205,8 @@ def __init__(self):
         self.inflow_module.initial()
         self.surface_routing_module.initial()
 
+        # At this point LddChan and LddKinematic do not have any structure reservoirs/lakes MCT confluence
+
         self.reservoir_module.initial()
         self.lakes_module.initial()
         self.polder_module.initial()
@@ -209,17 +215,28 @@ def __init__(self):
 
         self.structures_module.initial()
         # Structures such as reservoirs and lakes are modelled by interrupting the channel flow paths
+        # At this point LddKinematic and LddChan have pits upstream of (structures) reservoirs and lakes
+
+        if option.get('MCTRoutingInterface'):
+            self.mctconfluence_module.initial()
+            # initialising MCT confluence points and adding MCT confluence sinks to the LDD
 
         # ----------------------------------------------------------------------
         # ----------------------------------------------------------------------
 
         self.routing_module.initialSecond()
         # CHANNEL INITIAL SPLIT UP IN SECOND CHANNEL
+
         self.surface_routing_module.initialSecond()
 
+        # MCT confluence must be in the LddKinematic at this point
+        self.routing_module.initialKinematicWave()
+
         if option.get('MCTRouting'):
             self.routing_module.initialMCT()
             # initialising Muskingum-Cunge-Todini routing for channel
+            if option.get('MCTRoutingInterface'):
+                self.mctconfluence_module.dynamic_init()
 
         self.evapowater_module.initial()
         self.riceirrigation_module.initial()

src/lisflood/global_modules/add1.py

@@ -183,7 +183,7 @@ def loadsetclone(name):
         # settings x is first
         # setclone row col cellsize xupleft yupleft
         try:
-            setclone(int(coord[1]), int(coord[0]), float(coord[2]), float(coord[3]), float(coord[4]))   # CM: pcraster
+            setclone(int(coord[1]), int(coord[0]), float(coord[2]), float(coord[3]), float(coord[4]))   # pcraster
         except:
             rem = "["+str(coord[0])+" "+ str(coord[1])+" "+ str(coord[2])+" "+ str(coord[3])+" "+str(coord[4])+"]"
             msg = "Maskmap: " + rem + \
@@ -357,7 +357,7 @@ def loadmap_base(name, pcr=False, lddflag=False, timestampflag='exact', averagey
 
     :param name: name of key in Settings.xml input file containing path and name of the map file (as string)
     :param pcr: flag for output maps in pcraster format 
-    :param lddflag: flag for local drain direction map (CM??)
+    :param lddflag: flag for local drain direction map
     :param timestampflag: look for exact time stamp in netcdf file ('exact') or for the closest (left) time stamp available ('closest')
     :param averageyearflag: if True, use "average year" netcdf file over the entire model simulation period
     :param force_load_with_nans: if True, loads the map without checking for nan values inside area Map. 
@@ -709,11 +709,11 @@ def readnetcdf(name, time, timestampflag='exact', averageyearflag=False):
     t_steps = nf1.variables['time'][:]    # get values for timesteps ([  0.,  24.,  48.,  72.,  96.])
     t_unit = nf1.variables['time'].units  # get unit (u'hours since 2015-01-01 06:00:00')
     t_cal = get_calendar_type(nf1)
-    # CM: get year from time unit in case average year is used
+    # get year from time unit in case average year is used
     if averageyearflag:
-        # CM: get date of the first step in netCDF file containing average year values
+        # get date of the first step in netCDF file containing average year values
         first_date = num2date(t_steps[0], t_unit, t_cal)
-        # CM: get year of the first step in netCDF file containing average year values
+        # get year of the first step in netCDF file containing average year values
         t_ref_year = first_date.year
     settings = LisSettings.instance()
     binding = settings.binding
@@ -733,7 +733,7 @@ def readnetcdf(name, time, timestampflag='exact', averageyearflag=False):
         try:
             currentDate = currentDate.replace(year=t_ref_year)
         except:
-            # CM: if simulation year is leap and average year is not, switch 29/2 with 28/2
+            # if simulation year is leap and average year is not, switch 29/2 with 28/2
             currentDate = currentDate.replace(day=28)
             currentDate = currentDate.replace(year=t_ref_year)
 
@@ -747,9 +747,9 @@ def readnetcdf(name, time, timestampflag='exact', averageyearflag=False):
             msg = "Date " + str(currentDate) + " not stored in " + filename
             raise LisfloodError(msg)
         elif (timestampflag == 'closest'):
-            # CM: get the closest value
+            # get the closest value
             current_ncdf_step_new = takeClosest(t_steps, current_ncdf_step)
-            # CM: set current_ncdf_step to the closest available time step in netCDF file
+            # set current_ncdf_step to the closest available time step in netCDF file
             current_ncdf_step = current_ncdf_step_new
 
     # get index of timestep in netCDF file corresponding to current simulation date

src/lisflood/global_modules/checkers.py

@@ -26,7 +26,7 @@
 from ..hydrological_modules import (surface_routing, evapowater, snow, routing, leafarea, inflow, waterlevel,
                                     waterbalance, wateruse, waterabstraction, lakes, riceirrigation, indicatorcalc,
                                     landusechange, frost, groundwater, miscInitial, soilloop, soil,
-                                    reservoir, transmission)
+                                    reservoir, transmission, mctconfluence)
 
 
 class ModulesInputs:

src/lisflood/global_modules/default_options.py

@@ -10,6 +10,7 @@
                    'MonteCarlo': False,
                    'SplitRouting': False,
                    'MCTRouting': False,
+                   'MCTRoutingInterface': False,
                    'TemperatureInKelvin': False,
                    'TransLoss': False,
                    'TransientLandUseChange': False,
@@ -1288,6 +1289,7 @@
                    'repwateruseGauges': False,
                    'repwateruseSites': False,
                    'riceIrrigation': False,
+                   'simulateCalibrationPoints': False,
                    'simulateLakes': False,
                    'simulatePF': False,
                    'simulatePolders': False,

src/lisflood/global_modules/netcdf.py

@@ -530,8 +530,8 @@ def write_netcdf_header(settings,
     # time coordinates and associated values
     if frequency is not None:  # output file with "time" dimension
         n_steps = len(rep_steps)
-        #Get initial and final dates for data to be stored in nerCDF file
-        # CM: Create time stamps for each step stored in netCDF file
+        # Get initial and final dates for data to be stored in nerCDF file
+        # Create time stamps for each step stored in netCDF file
         all_dates = np.array([start_date + datetime.timedelta(days=(int(d)-1)*DtDay) for d in rep_steps])
         all_steps = np.array(rep_steps)
         if frequency == "all":
@@ -556,16 +556,16 @@ def write_netcdf_header(settings,
         time = nf1.createVariable('time', float, ('time'))
         time.standard_name = 'time'
         time.calendar = binding["calendar_type"]
-        # CM: select the time unit according to model time step
+        # select the time unit according to model time step
         DtDay_in_sec = DtDay * 86400
         if DtDay_in_sec >= 86400:
             # Daily model time steps or larger
             time.units = 'days since %s' % start_date.strftime("%Y-%m-%d %H:%M:%S.0")
         elif DtDay_in_sec >= 3600 and DtDay_in_sec < 86400:
-            # CM: hours to days model time steps
+            # hours to days model time steps
             time.units = 'hours since %s' % start_date.strftime("%Y-%m-%d %H:%M:%S.0")
         elif DtDay_in_sec >= 60 and DtDay_in_sec <3600:
-            # CM: minutes to hours model time step
+            # minutes to hours model time step
             time.units = 'minutes since %s' % start_date.strftime("%Y-%m-%d %H:%M:%S.0")
         nf1.variables["time"][:] = date2num(time_stamps, time.units, time.calendar)
 

src/lisflood/global_modules/settings.py

@@ -507,7 +507,7 @@ def _out_dir(user_settings):
             else:
                 pathout = pathout.replace(pathout[a1:a2 + 1], s2)
 
-        # CM: output folder
+        # output folder
         return pathout
 
     @staticmethod
@@ -790,13 +790,13 @@ def inttodate(int_in, ref_date, binding=None):
         settings = LisSettings.instance()
         binding = settings.binding
 
-    # CM: get model time step as float form 'DtSec' in Settings.xml file
+    # get model time step as float form 'DtSec' in Settings.xml file
     DtSec = float(binding['DtSec'])
-    # CM: compute fraction of day corresponding to model time step as float
+    # compute fraction of day corresponding to model time step as float
     DtDay = DtSec / 86400.
     # Time step, expressed as fraction of day (same as self.var.DtSec and self.var.DtDay)
 
-    # CM: compute date corresponding to intIn steps from reference date refDate
+    # compute date corresponding to intIn steps from reference date refDate
     stepDate = ref_date + datetime.timedelta(days=(int_in * DtDay))
 
     return stepDate

src/lisflood/global_modules/stateVar.py

@@ -20,7 +20,7 @@
 from .add1 import *
 
 
-# CM: new-style class in Python 2.x
+# new-style class in Python 2.x
 class stateVar(object):
 
     """

src/lisflood/hydrological_modules/inflow.py

@@ -54,6 +54,7 @@ def initial(self):
         # ************************************************************
         settings = LisSettings.instance()
         option = settings.options
+
         if option['inflow']:
             self.var.InflowPoints = loadmap('InflowPoints') #1D array size is pixels belonging to basin mask
 
@@ -110,7 +111,7 @@ def dynamic_init(self):
         settings = LisSettings.instance()
         option = settings.options
         if option['inflow']:
-            self.var.QDelta = (self.var.QInM3 - self.var.QInM3Old) * self.var.InvNoRoutSteps
+            self.var.QDeltaM3 = (self.var.QInM3 - self.var.QInM3Old) * self.var.InvNoRoutSteps
             # difference between old and new inlet flow  per sub step
             # in order to calculate the amount of inlet flow in the routing loop
 
@@ -146,6 +147,6 @@ def dynamic_inloop(self, NoRoutingExecuted):
 
         if option['inflow']:
 
-            self.var.QInDt = (self.var.QInM3Old + (NoRoutingExecuted + 1) * self.var.QDelta) * self.var.InvNoRoutSteps
+            self.var.QInM3Dt = (self.var.QInM3Old + (NoRoutingExecuted + 1) * self.var.QDeltaM3) * self.var.InvNoRoutSteps
             # flow from inlets per sub step
-            self.var.QinADDEDM3 += self.var.QInDt 
+            self.var.QinADDEDM3 += self.var.QInM3Dt

src/lisflood/hydrological_modules/kinematic_wave_parallel.py

@@ -74,6 +74,7 @@ def streamLookups(flow_dir, land_mask):
     '''
     Compute the downstream lookup vector for a D8 water flow channel network,
     i.e. the adjecency list of the directed graph describing flow direction from each pixel.
+    Each land_mask pixel is assigned a unique id from 0 to numpix-1, numbered by row
     Arguments:
         flow_dir (numpy.ndarray): LISFLOOD flow matrix values (FLOW_CODE).
         land_mask (numpy.ndarray): land mask on coordinate mesh.
@@ -82,12 +83,20 @@ def streamLookups(flow_dir, land_mask):
         upstream lookup (numpy.ndarray): each row gives the immediately upstream pixels (-1 = fill value); size = num_pixels, max_ups_pixs <= 8
     '''
     flow_dir[~land_mask] = 8 # exceeds number of rows of IX_ADDS
+    # assign 8 to all pixels not belonging to land_mask
     num_pixs = land_mask.sum()
+    # count number of pixels in land_mask
+
     # Create 2D array of indices of land pixels (each index is unique)
     # The IDs are numbered from 0 to num_pixs-1
     land_points = -np.ones(land_mask.shape, int)
     land_points[land_mask] = np.arange(num_pixs, dtype=int)
+    # every land pixel has now a unique id (pixel id)
+
     downstream_lookup, upstream_lookup = kwpt.upDownLookups(flow_dir, np.ascontiguousarray(land_mask).astype(np.uint8), land_points, num_pixs, IX_ADDS)
+    # downstream_lookup[i] = index of the cell receiving flow from node i
+    # upstream_lookup[i, :] = indices of all cells draining into node i
+
     max_num_ups_pixs = max(1, np.any(upstream_lookup != -1, 0).sum()) # maximum number of upstreams pixels
     return downstream_lookup, np.ascontiguousarray(upstream_lookup[:,:max_num_ups_pixs]).astype(int)