Interpolation documentation

core/src/validation.jl

@@ -389,8 +389,8 @@ function valid_tabulated_curve_level(
             basin_bottom_level = basin_bottom(basin, id_in)[2]
             # the second level is the bottom, the first is added to control extrapolation
             if table.t[1] + 1.0 < basin_bottom_level
-                @error "Lowest levels of $id is lower than bottom of upstream $id_in" table.t[1] +
-                                                                                      1.0 basin_bottom_level
+                @error "Lowest level of $id is lower than bottom of upstream $id_in" table.t[1] +
+                                                                                     1.0 basin_bottom_level
                 errors = true
             end
         end

core/test/validation_test.jl

@@ -394,7 +394,7 @@ end
     @test length(logger.logs) == 1
     @test logger.logs[1].level == Error
     @test logger.logs[1].message ==
-          "Lowest levels of TabulatedRatingCurve #5 is lower than bottom of upstream Basin #1"
+          "Lowest level of TabulatedRatingCurve #5 is lower than bottom of upstream Basin #1"
 end
 
 @testitem "Outlet upstream level validation" begin

docs/_quarto.yml

@@ -66,7 +66,6 @@ website:
         - reference/index.qmd
         - reference/usage.qmd
         - reference/validation.qmd
-        - reference/allocation.qmd
         - reference/python/index.qmd
         - reference/test-models.qmd
         - section: "Nodes"

docs/reference/node/basin.qmd

@@ -32,8 +32,60 @@ to have a reasonable and safe default, a value must be provided in the static ta
 This table is the transient form of the `Basin` table.
 The only difference is that a time column is added.
 The table must by sorted by time, and per time it must be sorted by `node_id`.
+
+### Interpolation
+
 At the given timestamps the values are set in the simulation, such that the timeseries can be seen as forward filled.
 
+```{python}
+# | code-fold: true
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from IPython.display import display, Markdown
+
+np.random.seed(1)
+fig, ax = plt.subplots()
+fontsize = 15
+
+N = 5
+y = np.random.rand(N)
+x = np.cumsum(np.random.rand(N))
+
+def forward_fill(x_):
+    i = min(max(0, np.searchsorted(x, x_)-1), len(x)-1)
+    return y[i]
+
+def plot_forward_fill(i):
+    ax.plot([x[i], x[i+1]], [y[i], y[i]], color = "C0", label = "interpolation" if i == 0 else None)
+
+ax.scatter(x[:-1],y[:-1], label = "forcing at data points")
+for i in range(N-1):
+    plot_forward_fill(i)
+
+x_missing_data = np.sort(x[0] + (x[-1] - x[0]) * np.random.rand(5))
+y_missing_data = [forward_fill(x_) for x_ in x_missing_data]
+ax.scatter(x_missing_data, y_missing_data, color = "C0", marker = "x", label = "missing data")
+ax.set_xticks([])
+ax.set_yticks([])
+ax.set_xlabel("time", fontsize = fontsize)
+ax.set_ylabel("forcing", fontsize = fontsize)
+xlim = ax.get_xlim()
+ax.set_xlim(xlim[0], (x[-2] + x[-1])/2)
+ax.legend()
+
+markdown_table = pd.DataFrame(
+        data = {
+            "time" : x,
+            "forcing" : y
+        }
+    ).to_markdown(index = False)
+
+display(Markdown(markdown_table))
+```
+
+As shown this interpolation type supports missing data, and just maintains the last available value. Because of this for instance precipitation can be updated while evaporation stays the same.
+
 ## State {#sec-state}
 
 The state table gives the initial water levels of all Basins.
@@ -83,6 +135,220 @@ Internally this get converted to two functions, $A(S)$ and $h(S)$, by integratin
 The minimum area cannot be zero to avoid numerical issues.
 The maximum area is used to convert the precipitation flux into an inflow.
 
+### Interpolation
+
+#### Level to area
+
+The level to area relationship is defined with the `Basin / profile` data using linear interpolation. An example of such a relationship is shown below.
+
+```{python}
+# | code-fold: true
+fig, ax = plt.subplots()
+
+# Data
+N = 3
+area = 25 * np.cumsum(np.random.rand(N))
+level = np.cumsum(np.random.rand(N))
+
+# Interpolation
+ax.scatter(level,area, label = "data")
+ax.plot(level,area, label = "interpolation")
+ax.set_xticks([level[0], level[-1]])
+ax.set_xticklabels(["bottom", "last supplied level"])
+ax.set_xlabel("level", fontsize = fontsize)
+ax.set_ylabel("area", fontsize = fontsize)
+ax.set_yticks([0])
+
+# Extrapolation
+level_extrap = 2 * level[-1] - level[-2]
+area_extrap = 2 * area[-1] - area[-2]
+ax.plot([level[-1], level_extrap], [area[-1], area_extrap], color = "C0", ls = "dashed", label = "extrapolation")
+xlim = ax.get_xlim()
+ax.set_xlim(xlim[0], (level[-1] + level_extrap)/2)
+
+ax.legend()
+fig.tight_layout()
+
+markdown_table = pd.DataFrame(
+        data = {
+            "level" : level,
+            "area" : area
+        }
+    ).to_markdown(index = False)
+
+display(Markdown(markdown_table))
+```
+
+For this interpolation it is validated that:
+
+- The areas are positive and are non-decreasing;
+- There are at least 2 data points.
+
+This interpolation is used in each evaluation of the right hand side function of the ODE.
+
+#### Level to storage
+
+The level to storage relationship gives the volume of water in the basin at a given level, which is given by the integral over the level to area relationship from the basin bottom to the given level:
+
+$$
+    S(h) = \int_{h_0}^h A(h')\text{d}h'.
+$$
+
+```{python}
+# | code-fold: true
+storage = np.diff(level) * area[:-1] + 0.5 * np.diff(area) * np.diff(level)
+storage = np.cumsum(storage)
+storage = np.insert(storage, 0, 0.0)
+def S(h):
+    i = min(max(0, np.searchsorted(level, h)-1), len(level)-2)
+    return storage[i] + area[i] * (h - level[i]) + 0.5 * (area[i+1] - area[i]) / (level[i+1] - level[i]) * (h - level[i])**2
+
+S = np.vectorize(S)
+
+# Interpolation
+fig, ax = plt.subplots()
+level_eval = np.linspace(level[0], level[-1], 100)
+storage_eval = S(np.linspace(level[0], level[-1], 100))
+ax.scatter(level, storage, label = "storage at datapoints")
+ax.plot(level_eval, storage_eval, label = "interpolation")
+ax.set_xticks([level[0], level[-1]])
+ax.set_xticklabels(["bottom", "last supplied level"])
+ax.set_yticks([0])
+ax.set_xlabel("level", fontsize = fontsize)
+ax.set_ylabel("storage", fontsize = fontsize)
+
+# Extrapolation
+level_eval_extrap = np.linspace(level[-1], level_extrap, 35)
+storage_eval_extrap = S(level_eval_extrap)
+ax.plot(level_eval_extrap, storage_eval_extrap, color = "C0", linestyle = "dashed", label = "extrapolation")
+xlim = ax.get_xlim()
+ax.set_xlim(xlim[0], (level[-1] + level_extrap)/2)
+ax.legend()
+```
+
+for converting the initial state in terms of levels to an initial state in terms of storages used in the core.
+
+#### Interactive basin example
+
+The profile data is not detailed enough to create a full 3D picture of the basin. However, if we assume the profile data is for a stretch of canal of given length, the following plot shows a cross section of the basin.
+```{python}
+# | code-fold: true
+import plotly.graph_objects as go
+import numpy as np
+
+def linear_interpolation(X,Y,x):
+    i = min(max(0, np.searchsorted(X, x)-1), len(X)-2)
+    return Y[i] + (Y[i+1] - Y[i])/(X[i+1] - X[i]) * (x - X[i])
+
+def A(h):
+    return linear_interpolation(level, area, h)
+
+fig = go.Figure()
+
+x = area/2
+x = np.concat([-x[::-1], x])
+y = np.concat([level[::-1], level])
+
+# Basin profile
+fig.add_trace(
+    go.Scatter(
+        x = x,
+        y = y,
+        line = dict(color = "green"),
+        name = "Basin profile"
+    )
+)
+
+# Basin profile extrapolation
+y_extrap = np.array([level[-1], level_extrap])
+x_extrap = np.array([area[-1]/2, area_extrap/2])
+fig.add_trace(
+    go.Scatter(
+        x = x_extrap,
+        y = y_extrap,
+        line = dict(color = "green", dash = "dash"),
+        name = "Basin extrapolation"
+    )
+)
+fig.add_trace(
+    go.Scatter(
+        x = -x_extrap,
+        y = y_extrap,
+        line = dict(color = "green", dash = "dash"),
+        showlegend = False
+    )
+)
+
+# Water level
+fig.add_trace(
+    go.Scatter(x = [-area[0]/2, area[0]/2],
+               y = [level[0], level[0]],
+               line = dict(color = "blue"),
+               name= "Water level")
+)
+
+# Fill area
+fig.add_trace(
+    go.Scatter(
+        x = [],
+        y = [],
+        fill = 'tonexty',
+        fillcolor = 'rgba(0, 0, 255, 0.2)',
+        line = dict(color = 'rgba(255, 255, 255, 0)'),
+        name = "Filled area"
+    )
+)
+
+# Create slider steps
+steps = []
+for h in np.linspace(level[0], level_extrap, 100):
+    a = A(h)
+    s = S(h).item()
+
+
+    i = min(max(0, np.searchsorted(level, h)-1), len(level)-2)
+    fill_area = np.append(area[:i+1], a)
+    fill_level = np.append(level[:i+1], h)
+    fill_x = np.concat([-fill_area[::-1]/2, fill_area/2])
+    fill_y = np.concat([fill_level[::-1], fill_level])
+
+    step = dict(
+        method = "update",
+        args=[
+            {
+                "x": [x, x_extrap, -x_extrap, [-a/2, a/2], fill_x],
+                "y": [y, y_extrap, y_extrap, [h, h], fill_y]
+            },
+            {"title": f"Interactive water level <br> Area: {a:.2f}, Storage: {s:.2f}"}
+        ],
+        label=str(round(h, 2))
+    )
+    steps.append(step)
+
+# Create slider
+sliders = [dict(
+    active=0,
+    currentvalue={"prefix": "Level: "},
+    pad={"t": 25},
+    steps=steps
+)]
+
+fig.update_layout(
+    title = {
+        "text": f"Interactive water level <br> Area: {area[0]:.2f}, Storage: 0.0",
+    },
+    yaxis_title = "level",
+    sliders = sliders,
+    margin = {"t": 100, "b": 100}
+)
+
+fig.show()
+```
+
+#### Storage to level
+
+The level is computed from the storage by inverting the level to storage relationship shown above. See [here](https://docs.sciml.ai/DataInterpolations/stable/inverting_integrals/) for more details.
+
 ## Area
 
 The optional area table is not used during computation, but provides a place to associate areas in the form of polygons to Basins.

docs/reference/node/tabulated-rating-curve.qmd

@@ -31,6 +31,51 @@ Below the lowest given level of -0.10, the flow rate is kept at 0.
 Between given levels the flow rate is interpolated linearly.
 Above the maximum given level of 20.09, the flow rate keeps increases linearly according to the slope of the last segment.
 
+### Interpolation
+
+The $Q(h)$ relationship of a tabulated rating curve is defined as a linear interpolation.
+
+```{python}
+# | code-fold: true
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from IPython.display import display, Markdown
+
+np.random.seed(0)
+fontsize = 15
+
+N = 10
+level = np.cumsum(np.random.rand(N))
+flow = level**2 + np.random.rand(N)
+flow[0] = 0.0
+fig, ax = plt.subplots()
+ax.set_xticks([])
+ax.set_yticks([0])
+ax.scatter(level, flow, label = "data")
+ax.plot(level, flow, label = "interpolation", color = "C0")
+ax.plot([level[0] - 1, level[0]], [0, 0], label = "extrapolation", linestyle = "dashed")
+ax.legend()
+ax.set_xlabel("level", fontsize = fontsize)
+ax.set_ylabel("flow", fontsize = fontsize)
+
+level_extrap = 2 * level[-1] - level[-2]
+flow_extrap = 2 * flow[-1] - flow[-2]
+ax.plot([level[-1], level_extrap], [flow[-1], flow_extrap], color = "C0", linestyle = "dashed")
+ax.set_xlim(level[0] - 0.5, (level[-1] + level_extrap)/2)
+
+markdown_table = pd.DataFrame(
+        data = {
+            "level" : level,
+            "flow" : flow
+        }
+    ).to_markdown(index = False)
+
+display(Markdown(markdown_table))
+```
+
+Here it is validated that the flow starts at $0$ and is non-decreasing. The flow is extrapolated as $0$ backward and linearly forward.
+
 ## Time
 
 This table is the transient form of the `TabulatedRatingCurve` table.