diff --git a/.nojekyll b/.nojekyll index 94e2c2ef4..422552408 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -4831267a \ No newline at end of file +3c5d5384 \ No newline at end of file diff --git a/python/examples.html b/python/examples.html index 6419b96aa..d41161d5c 100644 --- a/python/examples.html +++ b/python/examples.html @@ -554,7 +554,7 @@

2 Update the basi ax = df_flow.pivot_table(index="time", columns="edge", values="flow_m3d").plot() ax.legend(bbox_to_anchor=(1.3, 1), title="Edge")
-
<matplotlib.legend.Legend at 0x7f33deef3e50>
+
<matplotlib.legend.Legend at 0x7fc375cb2f10>

diff --git a/search.json b/search.json index f53bb7422..6a5179a31 100644 --- a/search.json +++ b/search.json @@ -4,7 +4,7 @@ "href": "python/examples.html", "title": "Examples", "section": "", - "text": "1 Basic model with static forcing\n\nimport geopandas as gpd\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\n\nimport ribasim\n\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin,\n (1.0, 0.0), # 2: ManningResistance\n (2.0, 0.0), # 3: Basin\n (3.0, 0.0), # 4: TabulatedRatingCurve\n (3.0, 1.0), # 5: FractionalFlow\n (3.0, 2.0), # 6: Basin\n (4.0, 1.0), # 7: Pump\n (4.0, 0.0), # 8: FractionalFlow\n (5.0, 0.0), # 9: Basin\n (6.0, 0.0), # 10: LinearResistance\n (2.0, 2.0), # 11: LevelBoundary\n (2.0, 1.0), # 12: LinearResistance\n (3.0, -1.0), # 13: FractionalFlow\n (3.0, -2.0), # 14: Terminal\n (3.0, 3.0), # 15: FlowBoundary\n (0.0, 1.0), # 16: FlowBoundary\n (6.0, 1.0), # 17: LevelBoundary\n ]\n)\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"Basin\",\n \"ManningResistance\",\n \"Basin\",\n \"TabulatedRatingCurve\",\n \"FractionalFlow\",\n \"Basin\",\n \"Pump\",\n \"FractionalFlow\",\n \"Basin\",\n \"LinearResistance\",\n \"LevelBoundary\",\n \"LinearResistance\",\n \"FractionalFlow\",\n \"Terminal\",\n \"FlowBoundary\",\n \"FlowBoundary\",\n \"LevelBoundary\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n static=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array(\n [1, 2, 3, 4, 4, 5, 6, 8, 7, 9, 11, 12, 4, 13, 15, 16, 10], dtype=np.int64\n)\nto_id = np.array(\n [2, 3, 4, 5, 8, 6, 7, 9, 9, 10, 12, 3, 13, 14, 6, 1, 17], dtype=np.int64\n)\nlines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id)\nedge = ribasim.Edge(\n static=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": len(from_id) * [\"flow\"],\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1, 3, 3, 6, 6, 9, 9],\n \"area\": [0.0, 1000.0] * 4,\n \"level\": [0.0, 1.0] * 4,\n }\n)\n\n# Convert steady forcing to m/s\n# 2 mm/d precipitation, 1 mm/d evaporation\nseconds_in_day = 24 * 3600\nprecipitation = 0.002 / seconds_in_day\nevaporation = 0.001 / seconds_in_day\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [0],\n \"drainage\": [0.0],\n \"potential_evaporation\": [evaporation],\n \"infiltration\": [0.0],\n \"precipitation\": [precipitation],\n \"urban_runoff\": [0.0],\n }\n)\nstatic = static.iloc[[0, 0, 0, 0]]\nstatic[\"node_id\"] = [1, 3, 6, 9]\n\nbasin = ribasim.Basin(profile=profile, static=static)\n\nSetup linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\"node_id\": [10, 12], \"resistance\": [5e3, (3600.0 * 24) / 100.0]}\n )\n)\n\nSetup Manning resistance:\n\nmanning_resistance = ribasim.ManningResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [2],\n \"length\": [900.0],\n \"manning_n\": [0.04],\n \"profile_width\": [6.0],\n \"profile_slope\": [3.0],\n }\n )\n)\n\nSet up a rating curve node:\n\n# Discharge: lose 1% of storage volume per day at storage = 1000.0.\nq1000 = 1000.0 * 0.01 / seconds_in_day\n\nrating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\n \"node_id\": [4, 4],\n \"level\": [0.0, 1.0],\n \"discharge\": [0.0, q1000],\n }\n )\n)\n\nSetup fractional flows:\n\nfractional_flow = ribasim.FractionalFlow(\n static=pd.DataFrame(\n data={\n \"node_id\": [5, 8, 13],\n \"fraction\": [0.3, 0.6, 0.1],\n }\n )\n)\n\nSetup pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": [7],\n \"flow_rate\": [0.5 / 3600],\n }\n )\n)\n\nSetup level boundary:\n\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [11, 17],\n \"level\": [0.5, 1.5],\n }\n )\n)\n\nSetup flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [15, 16],\n \"flow_rate\": [-1e-4, 1e-4],\n }\n )\n)\n\nSetup terminal:\n\nterminal = ribasim.Terminal(\n static=pd.DataFrame(\n data={\n \"node_id\": [14],\n }\n )\n)\n\nSetup a model:\n\nmodel = ribasim.Model(\n modelname=\"basic\",\n node=node,\n edge=edge,\n basin=basin,\n level_boundary=level_boundary,\n flow_boundary=flow_boundary,\n pump=pump,\n linear_resistance=linear_resistance,\n manning_resistance=manning_resistance,\n tabulated_rating_curve=rating_curve,\n fractional_flow=fractional_flow,\n terminal=terminal,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"basic\")\n\n\n\n2 Update the basic model with transient forcing\nThis assumes you have already created the basic model with static forcing.\n\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\n\nimport ribasim\n\n\nmodel = ribasim.Model.from_toml(datadir / \"basic/basic.toml\")\n\n\ntime = pd.date_range(model.starttime, model.endtime)\nday_of_year = time.day_of_year.to_numpy()\nseconds_per_day = 24 * 60 * 60\nevaporation = (\n (-1.0 * np.cos(day_of_year / 365.0 * 2 * np.pi) + 1.0) * 0.0025 / seconds_per_day\n)\nrng = np.random.default_rng(seed=0)\nprecipitation = (\n rng.lognormal(mean=-1.0, sigma=1.7, size=time.size) * 0.001 / seconds_per_day\n)\n\nWe’ll use xarray to easily broadcast the values.\n\ntimeseries = (\n pd.DataFrame(\n data={\n \"node_id\": 1,\n \"time\": time,\n \"drainage\": 0.0,\n \"potential_evaporation\": evaporation,\n \"infiltration\": 0.0,\n \"precipitation\": precipitation,\n \"urban_runoff\": 0.0,\n }\n )\n .set_index(\"time\")\n .to_xarray()\n)\n\nbasin_ids = model.basin.static[\"node_id\"].to_numpy()\nbasin_nodes = xr.DataArray(\n np.ones(len(basin_ids)), coords={\"node_id\": basin_ids}, dims=[\"node_id\"]\n)\nforcing = (timeseries * basin_nodes).to_dataframe().reset_index()\n\n\nstate = pd.DataFrame(\n data={\n \"node_id\": basin_ids,\n \"storage\": 1000.0,\n \"concentration\": 0.0,\n }\n)\n\n\nmodel.basin.forcing = forcing\nmodel.basin.state = state\n\n\nmodel.modelname = \"basic-transient\"\nmodel.write(datadir / \"basic-transient\")\n\nNow run the model with ribasim basic-transient/basic.toml. After running the model, read back the output:\n\ndf_basin = pd.read_feather(datadir / \"basic-transient/output/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\ndf_basin_wide[\"level\"].plot()\n\n<Axes: xlabel='time'>\n\n\n\n\n\n\ndf_flow = pd.read_feather(datadir / \"basic-transient/output/flow.arrow\")\ndf_flow[\"edge\"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))\ndf_flow[\"flow_m3d\"] = df_flow.flow * 86400\nax = df_flow.pivot_table(index=\"time\", columns=\"edge\", values=\"flow_m3d\").plot()\nax.legend(bbox_to_anchor=(1.3, 1), title=\"Edge\")\n\n<matplotlib.legend.Legend at 0x7f33deef3e50>\n\n\n\n\n\n\ntype(df_flow)\n\npandas.core.frame.DataFrame\n\n\n\n\n3 Model with discrete control\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin\n (1.0, -1.0), # 2: LinearResistance\n (2.0, 0.0), # 3: Basin\n (1.0, 0.0), # 4: Pump\n (1.0, 1.0), # 5: Control\n ]\n)\n\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"Basin\",\n \"LinearResistance\",\n \"Basin\",\n \"Pump\",\n \"DiscreteControl\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n static=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array([1, 2, 1, 4, 5], dtype=np.int64)\nto_id = np.array([2, 3, 4, 3, 4], dtype=np.int64)\n\n# Note: currently only controlling pumps is supported in the Julia core\nedge_type = 4 * [\"flow\"] + [\"control\"]\n\nlines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id)\nedge = ribasim.Edge(\n static=gpd.GeoDataFrame(\n data={\"from_node_id\": from_id, \"to_node_id\": to_id, \"edge_type\": edge_type},\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1, 1, 3, 3, 3],\n \"area\": [0.0, 100.0, 100.0] * 2,\n \"level\": [0.0, 0.001, 1.0] * 2,\n }\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [1, 3],\n \"drainage\": [0.0] * 2,\n \"potential_evaporation\": [0.0] * 2,\n \"infiltration\": [0.0] * 2,\n \"precipitation\": [0.0] * 2,\n \"urban_runoff\": [0.0] * 2,\n }\n)\n\nstate = pd.DataFrame(data={\"node_id\": [1, 3], \"storage\": [100.0, 0.0]})\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the control:\n\ncondition = pd.DataFrame(\n data={\n \"node_id\": [5, 5],\n \"listen_feature_id\": [1, 3],\n \"variable\": [\"level\", \"level\"],\n \"greater_than\": [0.8, 0.4],\n }\n)\n\n# False, False -> \"on\"\n# True, False -> \"off\"\n# False, True -> \"off\"\n# True, True -> \"On\"\n\n# Truth state as subset of the conditions above and in that order\n\nlogic = pd.DataFrame(\n data={\n \"node_id\": [5, 5, 5, 5],\n \"truth_state\": [\"FF\", \"TF\", \"FT\", \"TT\"],\n \"control_state\": [\"on\", \"off\", \"off\", \"on\"],\n }\n)\n\ndiscrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)\n\nSetup the pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"control_state\": [\"off\", \"on\"],\n \"node_id\": [4, 4],\n \"flow_rate\": [0.0, 1e-5],\n }\n )\n)\n\nSetup the linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [2],\n \"resistance\": [1e5],\n }\n )\n)\n\nSetup a model:\n\nmodel = ribasim.Model(\n modelname=\"control\",\n node=node,\n edge=edge,\n basin=basin,\n linear_resistance=linear_resistance,\n pump=pump,\n discrete_control=discrete_control,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\n\n\n\n\n\nListen edges are plotted with a dashed line since they are not present in the “Edge / static” schema but only in the “Control / condition” schema.\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"control\")\n\n\nfrom matplotlib.dates import date2num\n\ndf_basin = pd.read_feather(datadir / \"control/output/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\n\nax = df_basin_wide[\"level\"].plot()\n\nax.hlines(\n model.discrete_control.condition.greater_than,\n df_basin.time[0],\n df_basin.time.max(),\n lw=1,\n ls=\":\",\n color=[\"C0\", \"C1\"],\n)\n\ndf_control = pd.read_feather(datadir / \"control/output/control.arrow\")\n\nax.vlines(df_control.time, 0, 1, lw=1, ls=\":\", color=\"k\")\nax.set_xticks(\n list(ax.get_xticks()) + date2num(df_control.time).tolist(),\n list(ax.get_xticklabels()) + df_control.control_state.tolist(),\n rotation=50,\n);\n\n\n\n\nLet’s print an overview of what happened with control:\n\nmodel.print_discrete_control_record(datadir / \"control/output/control.arrow\")\n\n0. At 2020-01-01 00:00:00 the control node with ID 5 reached truth state TF:\n For node ID 1 (Basin): level > 0.8\n For node ID 3 (Basin): level < 0.4\n\n This yielded control state \"off\":\n For node ID 4 (Pump): flow_rate = 0.0\n\n1. At 2020-01-30 14:26:08.048000 the control node with ID 5 reached truth state FF:\n For node ID 1 (Basin): level < 0.8\n For node ID 3 (Basin): level < 0.4\n\n This yielded control state \"on\":\n For node ID 4 (Pump): flow_rate = 1e-05\n\n2. At 2020-02-16 05:05:05.564000 the control node with ID 5 reached truth state FT:\n For node ID 1 (Basin): level < 0.8\n For node ID 3 (Basin): level > 0.4\n\n This yielded control state \"off\":\n For node ID 4 (Pump): flow_rate = 0.0" + "text": "1 Basic model with static forcing\n\nimport geopandas as gpd\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\n\nimport ribasim\n\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin,\n (1.0, 0.0), # 2: ManningResistance\n (2.0, 0.0), # 3: Basin\n (3.0, 0.0), # 4: TabulatedRatingCurve\n (3.0, 1.0), # 5: FractionalFlow\n (3.0, 2.0), # 6: Basin\n (4.0, 1.0), # 7: Pump\n (4.0, 0.0), # 8: FractionalFlow\n (5.0, 0.0), # 9: Basin\n (6.0, 0.0), # 10: LinearResistance\n (2.0, 2.0), # 11: LevelBoundary\n (2.0, 1.0), # 12: LinearResistance\n (3.0, -1.0), # 13: FractionalFlow\n (3.0, -2.0), # 14: Terminal\n (3.0, 3.0), # 15: FlowBoundary\n (0.0, 1.0), # 16: FlowBoundary\n (6.0, 1.0), # 17: LevelBoundary\n ]\n)\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"Basin\",\n \"ManningResistance\",\n \"Basin\",\n \"TabulatedRatingCurve\",\n \"FractionalFlow\",\n \"Basin\",\n \"Pump\",\n \"FractionalFlow\",\n \"Basin\",\n \"LinearResistance\",\n \"LevelBoundary\",\n \"LinearResistance\",\n \"FractionalFlow\",\n \"Terminal\",\n \"FlowBoundary\",\n \"FlowBoundary\",\n \"LevelBoundary\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n static=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array(\n [1, 2, 3, 4, 4, 5, 6, 8, 7, 9, 11, 12, 4, 13, 15, 16, 10], dtype=np.int64\n)\nto_id = np.array(\n [2, 3, 4, 5, 8, 6, 7, 9, 9, 10, 12, 3, 13, 14, 6, 1, 17], dtype=np.int64\n)\nlines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id)\nedge = ribasim.Edge(\n static=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": len(from_id) * [\"flow\"],\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1, 3, 3, 6, 6, 9, 9],\n \"area\": [0.0, 1000.0] * 4,\n \"level\": [0.0, 1.0] * 4,\n }\n)\n\n# Convert steady forcing to m/s\n# 2 mm/d precipitation, 1 mm/d evaporation\nseconds_in_day = 24 * 3600\nprecipitation = 0.002 / seconds_in_day\nevaporation = 0.001 / seconds_in_day\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [0],\n \"drainage\": [0.0],\n \"potential_evaporation\": [evaporation],\n \"infiltration\": [0.0],\n \"precipitation\": [precipitation],\n \"urban_runoff\": [0.0],\n }\n)\nstatic = static.iloc[[0, 0, 0, 0]]\nstatic[\"node_id\"] = [1, 3, 6, 9]\n\nbasin = ribasim.Basin(profile=profile, static=static)\n\nSetup linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\"node_id\": [10, 12], \"resistance\": [5e3, (3600.0 * 24) / 100.0]}\n )\n)\n\nSetup Manning resistance:\n\nmanning_resistance = ribasim.ManningResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [2],\n \"length\": [900.0],\n \"manning_n\": [0.04],\n \"profile_width\": [6.0],\n \"profile_slope\": [3.0],\n }\n )\n)\n\nSet up a rating curve node:\n\n# Discharge: lose 1% of storage volume per day at storage = 1000.0.\nq1000 = 1000.0 * 0.01 / seconds_in_day\n\nrating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\n \"node_id\": [4, 4],\n \"level\": [0.0, 1.0],\n \"discharge\": [0.0, q1000],\n }\n )\n)\n\nSetup fractional flows:\n\nfractional_flow = ribasim.FractionalFlow(\n static=pd.DataFrame(\n data={\n \"node_id\": [5, 8, 13],\n \"fraction\": [0.3, 0.6, 0.1],\n }\n )\n)\n\nSetup pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": [7],\n \"flow_rate\": [0.5 / 3600],\n }\n )\n)\n\nSetup level boundary:\n\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [11, 17],\n \"level\": [0.5, 1.5],\n }\n )\n)\n\nSetup flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [15, 16],\n \"flow_rate\": [-1e-4, 1e-4],\n }\n )\n)\n\nSetup terminal:\n\nterminal = ribasim.Terminal(\n static=pd.DataFrame(\n data={\n \"node_id\": [14],\n }\n )\n)\n\nSetup a model:\n\nmodel = ribasim.Model(\n modelname=\"basic\",\n node=node,\n edge=edge,\n basin=basin,\n level_boundary=level_boundary,\n flow_boundary=flow_boundary,\n pump=pump,\n linear_resistance=linear_resistance,\n manning_resistance=manning_resistance,\n tabulated_rating_curve=rating_curve,\n fractional_flow=fractional_flow,\n terminal=terminal,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"basic\")\n\n\n\n2 Update the basic model with transient forcing\nThis assumes you have already created the basic model with static forcing.\n\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\n\nimport ribasim\n\n\nmodel = ribasim.Model.from_toml(datadir / \"basic/basic.toml\")\n\n\ntime = pd.date_range(model.starttime, model.endtime)\nday_of_year = time.day_of_year.to_numpy()\nseconds_per_day = 24 * 60 * 60\nevaporation = (\n (-1.0 * np.cos(day_of_year / 365.0 * 2 * np.pi) + 1.0) * 0.0025 / seconds_per_day\n)\nrng = np.random.default_rng(seed=0)\nprecipitation = (\n rng.lognormal(mean=-1.0, sigma=1.7, size=time.size) * 0.001 / seconds_per_day\n)\n\nWe’ll use xarray to easily broadcast the values.\n\ntimeseries = (\n pd.DataFrame(\n data={\n \"node_id\": 1,\n \"time\": time,\n \"drainage\": 0.0,\n \"potential_evaporation\": evaporation,\n \"infiltration\": 0.0,\n \"precipitation\": precipitation,\n \"urban_runoff\": 0.0,\n }\n )\n .set_index(\"time\")\n .to_xarray()\n)\n\nbasin_ids = model.basin.static[\"node_id\"].to_numpy()\nbasin_nodes = xr.DataArray(\n np.ones(len(basin_ids)), coords={\"node_id\": basin_ids}, dims=[\"node_id\"]\n)\nforcing = (timeseries * basin_nodes).to_dataframe().reset_index()\n\n\nstate = pd.DataFrame(\n data={\n \"node_id\": basin_ids,\n \"storage\": 1000.0,\n \"concentration\": 0.0,\n }\n)\n\n\nmodel.basin.forcing = forcing\nmodel.basin.state = state\n\n\nmodel.modelname = \"basic-transient\"\nmodel.write(datadir / \"basic-transient\")\n\nNow run the model with ribasim basic-transient/basic.toml. After running the model, read back the output:\n\ndf_basin = pd.read_feather(datadir / \"basic-transient/output/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\ndf_basin_wide[\"level\"].plot()\n\n<Axes: xlabel='time'>\n\n\n\n\n\n\ndf_flow = pd.read_feather(datadir / \"basic-transient/output/flow.arrow\")\ndf_flow[\"edge\"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))\ndf_flow[\"flow_m3d\"] = df_flow.flow * 86400\nax = df_flow.pivot_table(index=\"time\", columns=\"edge\", values=\"flow_m3d\").plot()\nax.legend(bbox_to_anchor=(1.3, 1), title=\"Edge\")\n\n<matplotlib.legend.Legend at 0x7fc375cb2f10>\n\n\n\n\n\n\ntype(df_flow)\n\npandas.core.frame.DataFrame\n\n\n\n\n3 Model with discrete control\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin\n (1.0, -1.0), # 2: LinearResistance\n (2.0, 0.0), # 3: Basin\n (1.0, 0.0), # 4: Pump\n (1.0, 1.0), # 5: Control\n ]\n)\n\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"Basin\",\n \"LinearResistance\",\n \"Basin\",\n \"Pump\",\n \"DiscreteControl\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n static=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array([1, 2, 1, 4, 5], dtype=np.int64)\nto_id = np.array([2, 3, 4, 3, 4], dtype=np.int64)\n\n# Note: currently only controlling pumps is supported in the Julia core\nedge_type = 4 * [\"flow\"] + [\"control\"]\n\nlines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id)\nedge = ribasim.Edge(\n static=gpd.GeoDataFrame(\n data={\"from_node_id\": from_id, \"to_node_id\": to_id, \"edge_type\": edge_type},\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1, 1, 3, 3, 3],\n \"area\": [0.0, 100.0, 100.0] * 2,\n \"level\": [0.0, 0.001, 1.0] * 2,\n }\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [1, 3],\n \"drainage\": [0.0] * 2,\n \"potential_evaporation\": [0.0] * 2,\n \"infiltration\": [0.0] * 2,\n \"precipitation\": [0.0] * 2,\n \"urban_runoff\": [0.0] * 2,\n }\n)\n\nstate = pd.DataFrame(data={\"node_id\": [1, 3], \"storage\": [100.0, 0.0]})\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the control:\n\ncondition = pd.DataFrame(\n data={\n \"node_id\": [5, 5],\n \"listen_feature_id\": [1, 3],\n \"variable\": [\"level\", \"level\"],\n \"greater_than\": [0.8, 0.4],\n }\n)\n\n# False, False -> \"on\"\n# True, False -> \"off\"\n# False, True -> \"off\"\n# True, True -> \"On\"\n\n# Truth state as subset of the conditions above and in that order\n\nlogic = pd.DataFrame(\n data={\n \"node_id\": [5, 5, 5, 5],\n \"truth_state\": [\"FF\", \"TF\", \"FT\", \"TT\"],\n \"control_state\": [\"on\", \"off\", \"off\", \"on\"],\n }\n)\n\ndiscrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)\n\nSetup the pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"control_state\": [\"off\", \"on\"],\n \"node_id\": [4, 4],\n \"flow_rate\": [0.0, 1e-5],\n }\n )\n)\n\nSetup the linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [2],\n \"resistance\": [1e5],\n }\n )\n)\n\nSetup a model:\n\nmodel = ribasim.Model(\n modelname=\"control\",\n node=node,\n edge=edge,\n basin=basin,\n linear_resistance=linear_resistance,\n pump=pump,\n discrete_control=discrete_control,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\n\n\n\n\n\nListen edges are plotted with a dashed line since they are not present in the “Edge / static” schema but only in the “Control / condition” schema.\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"control\")\n\n\nfrom matplotlib.dates import date2num\n\ndf_basin = pd.read_feather(datadir / \"control/output/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\n\nax = df_basin_wide[\"level\"].plot()\n\nax.hlines(\n model.discrete_control.condition.greater_than,\n df_basin.time[0],\n df_basin.time.max(),\n lw=1,\n ls=\":\",\n color=[\"C0\", \"C1\"],\n)\n\ndf_control = pd.read_feather(datadir / \"control/output/control.arrow\")\n\nax.vlines(df_control.time, 0, 1, lw=1, ls=\":\", color=\"k\")\nax.set_xticks(\n list(ax.get_xticks()) + date2num(df_control.time).tolist(),\n list(ax.get_xticklabels()) + df_control.control_state.tolist(),\n rotation=50,\n);\n\n\n\n\nLet’s print an overview of what happened with control:\n\nmodel.print_discrete_control_record(datadir / \"control/output/control.arrow\")\n\n0. At 2020-01-01 00:00:00 the control node with ID 5 reached truth state TF:\n For node ID 1 (Basin): level > 0.8\n For node ID 3 (Basin): level < 0.4\n\n This yielded control state \"off\":\n For node ID 4 (Pump): flow_rate = 0.0\n\n1. At 2020-01-30 14:26:08.048000 the control node with ID 5 reached truth state FF:\n For node ID 1 (Basin): level < 0.8\n For node ID 3 (Basin): level < 0.4\n\n This yielded control state \"on\":\n For node ID 4 (Pump): flow_rate = 1e-05\n\n2. At 2020-02-16 05:05:05.564000 the control node with ID 5 reached truth state FT:\n For node ID 1 (Basin): level < 0.8\n For node ID 3 (Basin): level > 0.4\n\n This yielded control state \"off\":\n For node ID 4 (Pump): flow_rate = 0.0" }, { "objectID": "python/reference/Node.html",