From d1082f4b480d940415ef27ee4192b16b2e2cb620 Mon Sep 17 00:00:00 2001 From: Marco Dal Molin Date: Tue, 19 Nov 2019 10:09:51 +0100 Subject: [PATCH] Demo added --- demos/ZHydro_GR4J.ipynb | 501 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 501 insertions(+) create mode 100644 demos/ZHydro_GR4J.ipynb diff --git a/demos/ZHydro_GR4J.ipynb b/demos/ZHydro_GR4J.ipynb new file mode 100644 index 0000000..dda7596 --- /dev/null +++ b/demos/ZHydro_GR4J.ipynb @@ -0,0 +1,501 @@ +{ + "nbformat": 4, + "nbformat_minor": 0, + "metadata": { + "colab": { + "name": "SuperflexPy_demo_ZHydro.ipynb", + "provenance": [], + "collapsed_sections": [] + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + } + }, + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "FLE3y_HYMquG", + "colab_type": "text" + }, + "source": [ + "# SuperflexPy: A new open source framework for building conceptual hydrological models\n", + "\n", + "Authors: Marco Dal Molin, Dmitri Kavetski, Mario Schirmer, and Fabrizio Fenicia\n", + "\n", + "This is a tutorial for the implementation of the GR4J model using the SuperflexPy framework.\n", + "\n", + "For reference, look at the [official guide](https://https://superflexpy.readthedocs.io/)\n", + "\n", + "Framework code available on [Github](https://github.com/dalmo1991/superflexPy)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "hmPg1xOiyAut", + "colab_type": "text" + }, + "source": [ + "![alt text](https://github.com/dalmo1991/superflexPy/raw/master/doc/pics/GR4J.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Bnqev2ztNzwc", + "colab_type": "text" + }, + "source": [ + "## 01 - Installation" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "fopPfb-533oo", + "colab_type": "code", + "outputId": "799f02e5-11a3-461f-8316-1b2a04cd709a", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 34 + } + }, + "source": [ + "!pip install superflexpy" + ], + "execution_count": 0, + "outputs": [ + { + "output_type": "stream", + "text": [ + "Requirement already satisfied: superflexpy in /usr/local/lib/python3.6/dist-packages (0.2.2)\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "IkX3vXruN9Ws", + "colab_type": "text" + }, + "source": [ + "## 02 - Import the elements" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "mukW7GgJ-55L", + "colab_type": "code", + "colab": {} + }, + "source": [ + "# Elements of the framework\n", + "from superflexpy.framework.unit import Unit\n", + "\n", + "# Structural elements\n", + "from superflexpy.implementation.elements.structure_elements import Splitter, Junction, Transparent\n", + "\n", + "# GR4J elements\n", + "from superflexpy.implementation.elements.gr4j import InterceptionFilter, ProductionStore, UnitHydrograph1, UnitHydrograph2, RoutingStore, FluxAggregator\n", + "\n", + "# Root finder\n", + "from superflexpy.implementation.computation.pegasus_root_finding import PegasusPython" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "eo-f5aSYPrlx", + "colab_type": "text" + }, + "source": [ + "## 03 - Initialize the elements" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "4cpeZE3fPu9h", + "colab_type": "code", + "colab": {} + }, + "source": [ + "# Parameters of GR4J -> the others are set to the default value\n", + "x1, x2, x3, x4 = (50.0, 0.1, 20.0, 3.5)\n", + "\n", + "# Root finder (used to solve the differentail equations)\n", + "solver_pegasus = PegasusPython() # Use the default parameters\n", + "\n", + "# Interception Filter\n", + "\u000binterception_filter = InterceptionFilter(id='ir')\n", + "\n", + "# Production Store\n", + "production_store = ProductionStore(parameters={'x1': x1, 'alpha': 2.0, \n", + " 'beta': 5.0, 'ni': 4/9},\n", + " states={'S0': 10.0},\n", + " solver=solver_pegasus,\n", + " id='ps')\n", + "\n", + "# Splitter\n", + "splitter = Splitter(weight=[[0.9], [0.1]], \n", + " direction=[[0], [0]],\n", + " id='spl')\n", + "\n", + "# Unit Hydrograph 1\n", + "unit_hydrograph_1 = UnitHydrograph1(parameters={'lag-time': x4},\n", + " states={'lag': None},\n", + " id='uh1')\n", + "\n", + "# Unit Hydrograph 2\n", + "unit_hydrograph_2 = UnitHydrograph2(parameters={'lag-time': 2*x4},\n", + " states={'lag': None},\n", + " id='uh2')\n", + "\n", + "# Routing Store\n", + "routing_store = RoutingStore(parameters={'x2': x2, 'x3': x3,\n", + " 'gamma': 5.0, 'omega': 3.5},\n", + " states={'S0': 10.0},\n", + " solver=solver_pegasus,\n", + " id='rs')\n", + "\n", + "# Transparent Element\n", + "transparent = Transparent(id='tr') \n", + "\n", + "# Junction\n", + "junction = Junction(direction=[[0, None], # First output\n", + " [1, None], # Second output\n", + " [None, 0]], # Third output\n", + " id='jun')\n", + "\n", + "# Flux Aggregator\n", + "flux_aggregator = FluxAggregator(id='fa')\n", + "\n", + "# Build the unit\n", + "\u000bstructure = Unit(layers=[[interception_filter],\n", + " [production_store],\n", + " [splitter],\n", + " [unit_hydrograph_1, unit_hydrograph_2],\n", + " [routing_store, transparent],\n", + " [junction],\n", + " [flux_aggregator]],\n", + " id='structure')" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "IahQgEuqS12e", + "colab_type": "text" + }, + "source": [ + "## 04 - Play with the parameters and the states\n", + "\n", + "Once the elements have been initialized, the parameters and the states can be updated. Note how the identifiers of the element and of the unit are added to the original name of the parameter or state." + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "FcDnOgL1SnGO", + "colab_type": "code", + "outputId": "61f640fa-6e3a-4dcb-e9c9-92b41583b5e7", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 54 + } + }, + "source": [ + "# Get all the parameters\n", + "print(structure.get_parameters())" + ], + "execution_count": 0, + "outputs": [ + { + "output_type": "stream", + "text": [ + "{'structure_ps_x1': 50.0, 'structure_ps_alpha': 2.0, 'structure_ps_beta': 5.0, 'structure_ps_ni': 0.4444444444444444, 'structure_uh1_lag-time': 3.5, 'structure_uh2_lag-time': 7.0, 'structure_rs_x2': 0.1, 'structure_rs_x3': 20.0, 'structure_rs_gamma': 5.0, 'structure_rs_omega': 3.5}\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "axlKjZx6Tbfi", + "colab_type": "code", + "outputId": "f5506376-8187-4ad5-ad0e-94433be369b0", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 72 + } + }, + "source": [ + "# Change the parameter 'structure_ps_x1'\n", + "\n", + "structure.set_parameters(parameters={'structure_ps_x1' : 45.0})\n", + "\n", + "# Notice how the elements are copied into the structure and not linked: the\n", + "# parameter has changed only in the structure and not in the element.\n", + "\n", + "print(production_store.get_parameters())\n", + "print(structure.get_parameters())" + ], + "execution_count": 0, + "outputs": [ + { + "output_type": "stream", + "text": [ + "{'ps_x1': 50.0, 'ps_alpha': 2.0, 'ps_beta': 5.0, 'ps_ni': 0.4444444444444444}\n", + "{'structure_ps_x1': 45.0, 'structure_ps_alpha': 2.0, 'structure_ps_beta': 5.0, 'structure_ps_ni': 0.4444444444444444, 'structure_uh1_lag-time': 3.5, 'structure_uh2_lag-time': 7.0, 'structure_rs_x2': 0.1, 'structure_rs_x3': 20.0, 'structure_rs_gamma': 5.0, 'structure_rs_omega': 3.5}\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "RibvfPG4UTMy", + "colab_type": "text" + }, + "source": [ + "## 05 - Run the structure" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "v1CBzl2OUbav", + "colab_type": "code", + "outputId": "3ef5dae8-5ee7-4ddc-d5f7-5882cd631082", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 72 + } + }, + "source": [ + "# Generate inputs\n", + "import numpy as np\n", + "\n", + "NUM_TIMESTEPS = 100\n", + "\n", + "# Create input fluxes\n", + "precipitation = np.random.rand(NUM_TIMESTEPS)*10 # Precipitation input\n", + "precipitation[10:20] = 0\n", + "precipitation[55:60] = 0\n", + "precipitation[70:] = 0\n", + "pet = np.array([2.0]*NUM_TIMESTEPS)\n", + "\n", + "# Set the timestep\n", + "structure.set_timestep(dt=1.0)\n", + "\n", + "# Set the inputs\n", + "structure.set_input([pet, precipitation])\n", + "\n", + "# Get the output\n", + "output = structure.get_output(solve=True)\n", + "\n", + "# Notice that now the states of the storages are the ones at the ones at the\n", + "# end of the simulation. If you run the model again, it will start from those\n", + "# states.\n", + "\n", + "print(structure.get_states())" + ], + "execution_count": 0, + "outputs": [ + { + "output_type": "stream", + "text": [ + "{'structure_ps_S0': 4.944337435652198, 'structure_uh1_lag': [array([9.91303081e-06, 5.76975894e-06, 2.02315437e-06, 0.00000000e+00])], 'structure_uh2_lag': [array([1.77692002e-06, 1.14351078e-06, 6.63589307e-07, 2.93493533e-07,\n", + " 9.43696994e-08, 1.53355828e-08, 0.00000000e+00])], 'structure_rs_S0': 8.658425007791461}\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "2Q8r00fLejc_", + "colab_type": "text" + }, + "source": [ + "## 06 - Inspect after running\n", + "\n", + "Note that when calling the method `get_output` we need to pass the argument `solve=False`" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "Z2izkAAYec7S", + "colab_type": "code", + "colab": {} + }, + "source": [ + "# Interception\n", + "output_interception = structure.call_internal(id='ir', method='get_output', solve=False)\n", + "\n", + "# Production\n", + "output_production = structure.call_internal(id='ps', method='get_output', solve=False)\n", + "aet_production = structure.call_internal(id='ps', method='get_aet')\n", + "state_production= structure.get_internal(id='ps', attribute='state_array')\n", + "\n", + "# Unit hydrographs\n", + "output_uh1 = structure.call_internal(id='uh1', method='get_output', solve=False)\n", + "output_uh2 = structure.call_internal(id='uh2', method='get_output', solve=False)\n", + "\n", + "# Routing\n", + "output_routing = structure.call_internal(id='rs', method='get_output', solve=False)\n", + "state_routing= structure.get_internal(id='rs', attribute='state_array')" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "XacV9fvnfpkz", + "colab_type": "text" + }, + "source": [ + "## 07 - Plot\n", + "\n", + "Click on the legend to turn on/off the traces" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "EDgjW70oXX0t", + "colab_type": "code", + "outputId": "e5844b31-6ae3-4f9e-b4fb-a2a4a8c9ddbd", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 767 + } + }, + "source": [ + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", + "fig = make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.2)\n", + "\n", + "fig.add_trace(trace=go.Bar(x=list(range(NUM_TIMESTEPS)), y=precipitation, name='precipitation'),\n", + " row=1, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output[0], name='Final Output'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output_interception[0], name='Interception_PET'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output_interception[1], name='Interception_P'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output_production[0], name='Production_Output'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=aet_production, name='Production_AET'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=state_production[:, 0], name='Production_Storage'),\n", + " row=3, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output_uh1[0], name='UH1_Output'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output_uh2[0], name='UH2_Output'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output_routing[0], name='Routing_Output'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=output_routing[1], name='Loss'),\n", + " row=2, col=1)\n", + "\n", + "fig.add_trace(trace=go.Scatter(x=list(range(NUM_TIMESTEPS)), y=state_routing[:, 0], name='Routing_Storage'),\n", + " row=3, col=1)\n", + "\n", + "fig.update_xaxes(title_text=\"Timestep\", row=1, col=1)\n", + "fig.update_yaxes(title_text=\"Precipitation [mm/timestep]\", row=1, col=1)\n", + "fig.update_xaxes(title_text=\"Timestep\", row=2, col=1)\n", + "fig.update_yaxes(title_text=\"Fluxes [mm/timestep]\", row=2, col=1)\n", + "fig.update_xaxes(title_text=\"Timestep\", row=3, col=1)\n", + "fig.update_yaxes(title_text=\"Storage [mm]\", row=3, col=1)\n", + "fig.update_layout(height=750)" + ], + "execution_count": 0, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "
\n", + " \n", + " \n", + " \n", + "
\n", + " \n", + "
\n", + "\n", + "" + ] + }, + "metadata": { + "tags": [] + } + } + ] + } + ] +} \ No newline at end of file