Muskingum routing in the IJssel - verification scriptΒΆ

Code adjusted from https://github.com/quaquel/epa1361_open. The output of this model is equal to my own model, as can be seen in this script.

Source: https://onlinelibrary.wiley.com/doi/abs/10.1111/jfr3.12532

Ciullo, A., de Bruijn, K. M., Kwakkel, J. H., & Klijn, F. (2019). Accounting for the uncertain effects of hydraulic interactions in optimising embankments heights: Proof of principle for the IJssel River. Journal of Flood Risk Management, e12532.

[1]:
import sys
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
[2]:
import generate_network
from functions_ijssel_muskingum import Muskingum
[3]:
G, dike_list = generate_network.get_network()
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-3-022e56b79809> in <module>
----> 1 G, dike_list = generate_network.get_network()

~/checkouts/readthedocs.org/user_builds/rna/checkouts/stable/docs/ijssel/generate_network.py in get_network()
      9
     10     # load data
---> 11     df = pd.read_excel('../../data/network-structure.xlsx')
     12     df = df.set_index('node')
     13     nodes = df.to_dict('index')

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/util/_decorators.py in wrapper(*args, **kwargs)
    206                 else:
    207                     kwargs[new_arg_name] = new_arg_value
--> 208             return func(*args, **kwargs)
    209
    210         return wrapper

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_base.py in read_excel(io, sheet_name, header, names, index_col, usecols, squeeze, dtype, engine, converters, true_values, false_values, skiprows, nrows, na_values, keep_default_na, verbose, parse_dates, date_parser, thousands, comment, skip_footer, skipfooter, convert_float, mangle_dupe_cols, **kwds)
    308
    309     if not isinstance(io, ExcelFile):
--> 310         io = ExcelFile(io, engine=engine)
    311     elif engine and engine != io.engine:
    312         raise ValueError(

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_base.py in __init__(self, io, engine)
    817         self._io = _stringify_path(io)
    818
--> 819         self._reader = self._engines[engine](self._io)
    820
    821     def __fspath__(self):

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_xlrd.py in __init__(self, filepath_or_buffer)
     19         err_msg = "Install xlrd >= 1.0.0 for Excel support"
     20         import_optional_dependency("xlrd", extra=err_msg)
---> 21         super().__init__(filepath_or_buffer)
     22
     23     @property

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_base.py in __init__(self, filepath_or_buffer)
    357             self.book = self.load_workbook(filepath_or_buffer)
    358         elif isinstance(filepath_or_buffer, str):
--> 359             self.book = self.load_workbook(filepath_or_buffer)
    360         else:
    361             raise ValueError(

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_xlrd.py in load_workbook(self, filepath_or_buffer)
     34             return open_workbook(file_contents=data)
     35         else:
---> 36             return open_workbook(filepath_or_buffer)
     37
     38     @property

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/xlrd/__init__.py in open_workbook(filename, logfile, verbosity, use_mmap, file_contents, encoding_override, formatting_info, on_demand, ragged_rows)
    109     else:
    110         filename = os.path.expanduser(filename)
--> 111         with open(filename, "rb") as f:
    112             peek = f.read(peeksz)
    113     if peek == b"PK\x03\x04": # a ZIP file

FileNotFoundError: [Errno 2] No such file or directory: '../../data/network-structure.xlsx'
[4]:
class DikeNetwork(object):
    def __init__(self):
        # planning steps
        self.num_events = 30

        # load network
        G, dike_list = generate_network.get_network()

        self.Qpeaks = 2000 #np.random.uniform(1000,16000,100)

        self.G = G
        self.dikelist = dike_list

    def printG(self):
        print(G.nodes.data())

    def getG(self):
        return G


    def init_node(self,value, time):
        init = np.repeat(value, len(time)).tolist()
        return init

    def _initialize_hydroloads(self, node, time, Q_0):
        #node['cumVol'], node['wl'], node['Qpol'], node['hbas'] = (
        #    self.init_node(0, time) for _ in range(4))
        node['Qin'], node['Qout'] = (self.init_node(Q_0, time) for _ in range(2))
        #node['status'] = self.init_node(False, time)
        #node['tbreach'] = np.nan
        return node

    def calc_wave(self,timestep=1):
        startnode = G.node['A.0']
        waveshape_id = 0
        Qpeak = self.Qpeaks#[0]
        dikelist = self.dikelist
        time = np.arange(0, startnode['Qevents_shape'].loc[waveshape_id].shape[0],
                             timestep)
        startnode['Qout'] = Qpeak * startnode['Qevents_shape'].loc[waveshape_id]

        # Initialize hydrological event:
        for key in dikelist:
            node = G.node[key]
            #Q_0 = int(G.node['A.0']['Qout'][0])
            Q_0 = G.node['A.0']['Qout'][0]
            self._initialize_hydroloads(node, time, Q_0)

        # Run the simulation:
        # Run over the discharge wave:
        for t in range(1, len(time)):
            # Run over each node of the branch:
            for n in range(0, len(dikelist)):
                # Select current node:
                node = G.node[dikelist[n]]
                if node['type'] == 'dike':
                    # Muskingum parameters:
                    C1 = node['C1']
                    C2 = node['C2']
                    C3 = node['C3']

                    prev_node = G.node[node['pnode']]
                    # Evaluate Q coming in a given node at time t:
                    node['Qin'][t] = Muskingum(C1, C2, C3,
                                                   prev_node['Qout'][t],
                                                   prev_node['Qout'][t - 1],
                                                   node['Qin'][t - 1])

                    node['Qout'][t] = node['Qin'][t]

    def __call__(self, timestep=1, **kwargs):
        G = self.G
        Qpeaks = self.Qpeaks
        dikelist = self.dikelist
[5]:
dikeNetwork = DikeNetwork()
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-5-e47df4b58fc1> in <module>
----> 1 dikeNetwork = DikeNetwork()

<ipython-input-4-f3ab11cb7e64> in __init__(self)
      5
      6         # load network
----> 7         G, dike_list = generate_network.get_network()
      8
      9         self.Qpeaks = 2000 #np.random.uniform(1000,16000,100)

~/checkouts/readthedocs.org/user_builds/rna/checkouts/stable/docs/ijssel/generate_network.py in get_network()
      9
     10     # load data
---> 11     df = pd.read_excel('../../data/network-structure.xlsx')
     12     df = df.set_index('node')
     13     nodes = df.to_dict('index')

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/util/_decorators.py in wrapper(*args, **kwargs)
    206                 else:
    207                     kwargs[new_arg_name] = new_arg_value
--> 208             return func(*args, **kwargs)
    209
    210         return wrapper

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_base.py in read_excel(io, sheet_name, header, names, index_col, usecols, squeeze, dtype, engine, converters, true_values, false_values, skiprows, nrows, na_values, keep_default_na, verbose, parse_dates, date_parser, thousands, comment, skip_footer, skipfooter, convert_float, mangle_dupe_cols, **kwds)
    308
    309     if not isinstance(io, ExcelFile):
--> 310         io = ExcelFile(io, engine=engine)
    311     elif engine and engine != io.engine:
    312         raise ValueError(

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_base.py in __init__(self, io, engine)
    817         self._io = _stringify_path(io)
    818
--> 819         self._reader = self._engines[engine](self._io)
    820
    821     def __fspath__(self):

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_xlrd.py in __init__(self, filepath_or_buffer)
     19         err_msg = "Install xlrd >= 1.0.0 for Excel support"
     20         import_optional_dependency("xlrd", extra=err_msg)
---> 21         super().__init__(filepath_or_buffer)
     22
     23     @property

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_base.py in __init__(self, filepath_or_buffer)
    357             self.book = self.load_workbook(filepath_or_buffer)
    358         elif isinstance(filepath_or_buffer, str):
--> 359             self.book = self.load_workbook(filepath_or_buffer)
    360         else:
    361             raise ValueError(

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/pandas/io/excel/_xlrd.py in load_workbook(self, filepath_or_buffer)
     34             return open_workbook(file_contents=data)
     35         else:
---> 36             return open_workbook(filepath_or_buffer)
     37
     38     @property

~/checkouts/readthedocs.org/user_builds/rna/envs/stable/lib/python3.7/site-packages/xlrd/__init__.py in open_workbook(filename, logfile, verbosity, use_mmap, file_contents, encoding_override, formatting_info, on_demand, ragged_rows)
    109     else:
    110         filename = os.path.expanduser(filename)
--> 111         with open(filename, "rb") as f:
    112             peek = f.read(peeksz)
    113     if peek == b"PK\x03\x04": # a ZIP file

FileNotFoundError: [Errno 2] No such file or directory: '../../data/network-structure.xlsx'
[6]:
dikeNetwork.calc_wave()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-81275308749e> in <module>
----> 1 dikeNetwork.calc_wave()

NameError: name 'dikeNetwork' is not defined
[7]:
G = dikeNetwork.getG()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-22b4b66c74d7> in <module>
----> 1 G = dikeNetwork.getG()

NameError: name 'dikeNetwork' is not defined
[8]:
#G.nodes['A.1']
[9]:
plt.figure(figsize=(20,10))
plt.plot(G.node['A.0']['Qout'])
df = pd.DataFrame({'Qin':G.node['A.0']['Qout']})

dikelist = dikeNetwork.dikelist
for n in range(0, len(dikelist)):
    node = G.node[dikelist[n]]
    plt.plot(node['Qin'])
    df[dikelist[n]] = node['Qin']
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-b770146659e7> in <module>
      1 plt.figure(figsize=(20,10))
----> 2 plt.plot(G.node['A.0']['Qout'])
      3 df = pd.DataFrame({'Qin':G.node['A.0']['Qout']})
      4
      5 dikelist = dikeNetwork.dikelist

NameError: name 'G' is not defined
<Figure size 1920x960 with 0 Axes>
[10]:
df
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-00cf07b74dcd> in <module>
----> 1 df

NameError: name 'df' is not defined