WSIMOD model demonstration - Oxford (.py)¶
Note - this script can also be opened in interactive Python if you wanted to play around. On the GitHub it is in docs/demo/scripts
-
3.1 Freshwater Treatment Works
3.2 Land
3.3 Demand
3.4 Reservoir
3.5 Distribution
3.6 Wastewater Treatment Works
3.7 Sewers
3.8 Groundwater
3.9 Node list
-
4.1 Arc parameters
4.2 Create the arcs
-
7.1 Validation
We will cover a demo WSIMOD case study¶
The glamorous town of Oxford will be our demo case study. Below, we will create these nodes and arcs, orchestrate them into a model, and run simulations.
Although GIS is pretty, the schematic below is a more accurate representation of what will be created. WSIMOD treats everything as a node or an arc.
Imports and forcing data¶
Import packages
import os
import pandas as pd
from wsimod.arcs.arcs import Arc
from wsimod.core import constants
from wsimod.nodes.catchment import Catchment
from wsimod.nodes.demand import ResidentialDemand
from wsimod.nodes.land import Land
from wsimod.nodes.nodes import Node
from wsimod.nodes.sewer import Sewer
from wsimod.nodes.storage import Groundwater, Reservoir
from wsimod.nodes.waste import Waste
from wsimod.nodes.wtw import FWTW, WWTW
from wsimod.orchestration.model import Model
os.environ["USE_PYGEOS"] = "0"
import geopandas as gpd
from matplotlib import pyplot as plt
from shapely.geometry import LineString
Load input data
# Select the root path for the data folder. Use the appropriate value for your case.
data_folder = os.path.join(os.path.abspath(""), "docs", "demo", "data")
input_fid = os.path.join(data_folder, "processed", "timeseries_data.csv")
input_data = pd.read_csv(input_fid)
input_data.loc[input_data.variable == "flow", "value"] *= constants.M3_S_TO_M3_DT
input_data.loc[input_data.variable == "precipitation", "value"] *= constants.MM_TO_M
input_data.date = pd.to_datetime(input_data.date)
data_input_dict = input_data.set_index(["variable", "date"]).value.to_dict()
data_input_dict = (
input_data.groupby("site")
.apply(lambda x: x.set_index(["variable", "date"]).value.to_dict())
.to_dict()
)
print(input_data.sample(10))
site date variable value 65444 cherwell 2012-12-12 boron 0.000045 60568 evenlode 2011-07-23 magnesium 0.004543 55338 ray 2009-03-13 calcium 0.144914 3270 ray 2010-02-24 temperature 4.457143 8395 evenlode 2012-03-22 silicon 0.001302 77719 oxford_land 2010-09-05 et0 0.002000 24561 cherwell 2012-08-19 nitrate 0.004608 32492 ray 2010-06-06 sulphate 0.100093 7673 evenlode 2010-03-31 silicon 0.004629 51733 thames 2011-04-15 potassium 0.003850
Input data is stored in dicts
print(data_input_dict["cherwell"][("boron", pd.to_datetime("2010-11-20"))])
6.985714285714286e-05
We select dates that are available in the input data
dates = input_data.date.unique()
dates = dates[dates.argsort()]
dates = [pd.Timestamp(x) for x in dates]
print(dates[0:10])
[Timestamp('2009-03-03 00:00:00'), Timestamp('2009-03-04 00:00:00'), Timestamp('2009-03-05 00:00:00'), Timestamp('2009-03-06 00:00:00'), Timestamp('2009-03-07 00:00:00'), Timestamp('2009-03-08 00:00:00'), Timestamp('2009-03-09 00:00:00'), Timestamp('2009-03-10 00:00:00'), Timestamp('2009-03-11 00:00:00'), Timestamp('2009-03-12 00:00:00')]
We can specify the pollutants. In this example we choose based on what pollutants we have input data for.
constants.POLLUTANTS = input_data.variable.unique().tolist()
constants.POLLUTANTS.remove("flow")
constants.POLLUTANTS.remove("precipitation")
constants.POLLUTANTS.remove("et0")
constants.NON_ADDITIVE_POLLUTANTS = ["temperature"]
constants.ADDITIVE_POLLUTANTS = list(
set(constants.POLLUTANTS).difference(constants.NON_ADDITIVE_POLLUTANTS)
)
constants.FLOAT_ACCURACY = 1e-11
print(constants.POLLUTANTS)
['temperature', 'silicon', 'fluoride', 'chloride', 'nitrate', 'sulphate', 'nitrogen', 'sodium', 'potassium', 'calcium', 'magnesium', 'boron']
Create nodes¶
For waste nodes, no parameters are needed, they are just the model outlet
thames_above_abingdon = Waste(name="thames_above_abingdon")
For junctions and abstraction locations, we can simply use the default nodes
farmoor_abstraction = Node(name="farmoor_abstraction")
evenlode_thames = Node(name="evenlode_thames")
cherwell_ray = Node(name="cherwell_ray")
cherwell_thames = Node(name="cherwell_thames")
thames_mixer = Node(name="thames_mixer")
For catchment nodes, we only need to specify the input data (as a dictionary format).
evenlode = Catchment(name="evenlode", data_input_dict=data_input_dict["evenlode"])
thames = Catchment(name="thames", data_input_dict=data_input_dict["thames"])
ray = Catchment(name="ray", data_input_dict=data_input_dict["ray"])
cherwell = Catchment(name="cherwell", data_input_dict=data_input_dict["cherwell"])
We can see that, even though we provided mimimal information (name and input data) each node comes with many predefined functions.
print(dir(evenlode))
['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'blend_vqip', 'check_basic', 'compare_vqip', 'concentration_to_total', 'copy_vqip', 'data_input_dict', 'ds_vqip', 'ds_vqip_c', 'empty_vqip', 'empty_vqip_predefined', 'end_timestep', 'end_timestep_', 'extract_vqip', 'extract_vqip_c', 'get_avail', 'get_connected', 'get_data_input', 'get_direction_arcs', 'get_flow', 'in_arcs', 'in_arcs_type', 'mass_balance', 'mass_balance_ds', 'mass_balance_in', 'mass_balance_out', 'name', 'node_mass_balance', 'out_arcs', 'out_arcs_type', 'pull_check', 'pull_check_abstraction', 'pull_check_basic', 'pull_check_deny', 'pull_check_handler', 'pull_distributed', 'pull_set', 'pull_set_abstraction', 'pull_set_deny', 'pull_set_handler', 'push_check', 'push_check_accept', 'push_check_basic', 'push_check_deny', 'push_check_handler', 'push_distributed', 'push_set', 'push_set_deny', 'push_set_handler', 'query_handler', 'reinit', 'route', 'sum_vqip', 't', 'total_in', 'total_out', 'total_to_concentration', 'unrouted_water', 'v_change_vqip', 'v_change_vqip_c', 'v_distill_vqip', 'v_distill_vqip_c']
Freshwater treatment works¶
Each type of node uses different parameters (see API reference). Below we create a freshwater treatment works (FWTW)
oxford_fwtw = FWTW(
service_reservoir_storage_capacity=1e5,
service_reservoir_storage_area=2e4,
treatment_throughput_capacity=4.5e4,
name="oxford_fwtw",
)
Each node type has different types of functionality available
print(dir(oxford_fwtw))
['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'blend_vqip', 'check_basic', 'compare_vqip', 'concentration_to_total', 'copy_vqip', 'current_input', 'data_input_dict', 'ds_vqip', 'ds_vqip_c', 'empty_vqip', 'empty_vqip_predefined', 'end_timestep', 'end_timestep_', 'extract_vqip', 'extract_vqip_c', 'get_connected', 'get_data_input', 'get_direction_arcs', 'get_excess_throughput', 'in_arcs', 'in_arcs_type', 'liquor', 'liquor_multiplier', 'mass_balance', 'mass_balance_ds', 'mass_balance_in', 'mass_balance_out', 'name', 'node_mass_balance', 'out_arcs', 'out_arcs_type', 'percent_solids', 'previous_pulled', 'process_parameters', 'pull_check', 'pull_check_basic', 'pull_check_deny', 'pull_check_fwtw', 'pull_check_handler', 'pull_distributed', 'pull_set', 'pull_set_deny', 'pull_set_fwtw', 'pull_set_handler', 'push_check', 'push_check_accept', 'push_check_basic', 'push_check_deny', 'push_check_handler', 'push_distributed', 'push_set', 'push_set_deny', 'push_set_handler', 'query_handler', 'reinit', 'service_reservoir_initial_storage', 'service_reservoir_storage_area', 'service_reservoir_storage_capacity', 'service_reservoir_storage_elevation', 'service_reservoir_tank', 'solids', 'sum_vqip', 't', 'total_deficit', 'total_in', 'total_out', 'total_pulled', 'total_to_concentration', 'treat_current_input', 'treat_water', 'treated', 'treatment_throughput_capacity', 'unpushed_sludge', 'v_change_vqip', 'v_change_vqip_c', 'v_distill_vqip', 'v_distill_vqip_c']
The FWTW node has a tank representing the service reservoirs, we can see that it has been initialised empty.
print(oxford_fwtw.service_reservoir_tank.storage)
{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}
If we try to pull water from the FWTW, it responds that there is no water to pull.
print(oxford_fwtw.pull_check({"volume": 10}))
{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}
If we add in some water, we see the pull check responds that water is available.
oxford_fwtw.service_reservoir_tank.storage["volume"] += 25
print(oxford_fwtw.pull_check({"volume": 10}))
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 10.0}
When we set a pull request, we see that we successfully receive the water and the tank is updated.
reply = oxford_fwtw.pull_set({"volume": 10})
print(reply)
print(oxford_fwtw.service_reservoir_tank.storage)
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 10.0} {'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 15.0}
Land¶
We will now create a land node, it is a bit involved so you might want to skip ahead to demand, or check out the land node tutorial
Data inputs are a single dictionary
land_inputs = data_input_dict["oxford_land"]
print(land_inputs[("precipitation", pd.to_datetime("2010-11-20"))])
print(land_inputs[("et0", pd.to_datetime("2010-11-20"))])
print(land_inputs[("temperature", pd.to_datetime("2009-10-15"))])
0.0003 0.002 13.107142857142858
Assign some default pollutant deposition values (kg/m2/d)
pollutant_deposition = {
"boron": 100e-10,
"calcium": 70e-7,
"chloride": 60e-10,
"fluoride": 0.2e-7,
"magnesium": 6e-7,
"nitrate": 2e-9,
"nitrogen": 4e-7,
"potassium": 7e-7,
"silicon": 7e-9,
"sodium": 30e-9,
"sulphate": 70e-7,
}
Create two surfaces as a list of dicts
surface = [
{
"type_": "PerviousSurface",
"area": 2e7,
"pollutant_load": pollutant_deposition,
"surface": "rural",
"field_capacity": 0.3,
"depth": 0.5,
"initial_storage": 2e7 * 0.4 * 0.5,
},
{
"type_": "ImperviousSurface",
"area": 1e7,
"pollutant_load": pollutant_deposition,
"surface": "urban",
"initial_storage": 5e6,
},
]
Create the land node from these surfaces and the input data
oxford_land = Land(surfaces=surface, name="oxford_land", data_input_dict=land_inputs)
We can see the land node has various tanks that have been initialised empty
print(oxford_land.surface_runoff.storage)
print(oxford_land.subsurface_runoff.storage)
print(oxford_land.percolation.storage)
{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0} {'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0} {'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}
We can see the surfaces have also been initialised, although they are not empty because we provided 'initial_storage' parameters.
rural_surface = oxford_land.get_surface("rural")
urban_surface = oxford_land.get_surface("urban")
print("{0}-{1}".format("rural", rural_surface.storage))
print("{0}-{1}".format("urban", urban_surface.storage))
rural-{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 4000000.0} urban-{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 5000000.0}
We can run a timestep of the land node with the 'run' command
oxford_land.t = pd.to_datetime("2012-12-22")
oxford_land.run()
We can see that the land and surface tanks have been updated
print(oxford_land.surface_runoff.storage)
print(oxford_land.subsurface_runoff.storage)
print(oxford_land.percolation.storage)
print("{0}-{1}".format("rural", rural_surface.storage))
print("{0}-{1}".format("urban", urban_surface.storage))
{'temperature': 3.521534008631931, 'silicon': 0.00020062235262918442, 'fluoride': 0.0005732067217976697, 'chloride': 0.0001719620165393009, 'nitrate': 5.732067217976698e-05, 'sulphate': 0.20062235262918443, 'nitrogen': 0.011464134435953396, 'sodium': 0.0008598100826965045, 'potassium': 0.020062235262918445, 'calcium': 0.20062235262918443, 'magnesium': 0.017196201653930095, 'boron': 0.00028660336089883487, 'volume': 12400.000000000002} {'temperature': 0.9722493642718524, 'silicon': 0.0017554455855053634, 'fluoride': 0.0050155588157296096, 'chloride': 0.0015046676447188828, 'nitrate': 0.000501555881572961, 'sulphate': 1.7554455855053634, 'nitrogen': 0.1003111763145922, 'sodium': 0.007523338223594414, 'potassium': 0.17554455855053633, 'calcium': 1.7554455855053634, 'magnesium': 0.1504667644718883, 'boron': 0.0025077794078648048, 'volume': 58900.0} {'temperature': 0.3349282031818326, 'silicon': 0.005868203814403644, 'fluoride': 0.016766296612581843, 'chloride': 0.005029888983774552, 'nitrate': 0.001676629661258184, 'sulphate': 5.868203814403645, 'nitrogen': 0.3353259322516368, 'sodium': 0.02514944491887276, 'potassium': 0.5868203814403645, 'calcium': 5.868203814403645, 'magnesium': 0.5029888983774553, 'boron': 0.008383148306290921, 'volume': 176700.00000000003} rural-{'temperature': 6.203473168254872, 'silicon': 0.1321757282474618, 'fluoride': 0.3776449378498909, 'chloride': 0.11329348135496725, 'nitrate': 0.037764493784989084, 'sulphate': 132.17572824746182, 'nitrogen': 7.552898756997817, 'sodium': 0.5664674067748363, 'potassium': 13.21757282474618, 'calcium': 132.17572824746182, 'magnesium': 11.329348135496726, 'boron': 0.18882246892494545, 'volume': 3980000.0} urban-{'temperature': 0.13343969249141666, 'silicon': 0.06999999999999999, 'fluoride': 0.2, 'chloride': 0.06, 'nitrate': 0.02, 'sulphate': 70.0, 'nitrogen': 4.0, 'sodium': 0.3, 'potassium': 7.0, 'calcium': 70.0, 'magnesium': 6.0, 'boron': 0.1, 'volume': 5104000.0}
Residential demand¶
The residential demand node requires population, per capita demand and a pollutant_load dictionary that defines how much (weight in kg) pollution is generated per person per day.
oxford = ResidentialDemand(
name="oxford",
population=2e5,
per_capita=0.15,
pollutant_load={
"boron": 500 * constants.UG_L_TO_KG_M3 * 0.15,
"calcium": 150 * constants.MG_L_TO_KG_M3 * 0.15,
"chloride": 180 * constants.MG_L_TO_KG_M3 * 0.15,
"fluoride": 0.4 * constants.MG_L_TO_KG_M3 * 0.15,
"magnesium": 30 * constants.MG_L_TO_KG_M3 * 0.15,
"nitrate": 60 * constants.MG_L_TO_KG_M3 * 0.15,
"nitrogen": 50 * constants.MG_L_TO_KG_M3 * 0.15,
"potassium": 30 * constants.MG_L_TO_KG_M3 * 0.15,
"silicon": 20 * constants.MG_L_TO_KG_M3 * 0.15,
"sodium": 200 * constants.MG_L_TO_KG_M3 * 0.15,
"sulphate": 250 * constants.MG_L_TO_KG_M3 * 0.15,
"temperature": 14,
},
data_input_dict=land_inputs,
)
# pollutant_load calculated based on expected effluent at WWTW
Reservoir¶
A reservoir node is used to make abstractions from rivers and supply FWTWs
farmoor = Reservoir(
name="farmoor", capacity=1e7, initial_storage=1e7, area=1.5e6, datum=62
)
distribution = Node(name="oxford_distribution")
Wastewater treatment works¶
Wastewater treatment works (WWTW) are nodes that can store sewage water temporarily in storm tanks, and reduce the pollution amounts in water before releasing them onwards to rivers.
oxford_wwtw = WWTW(
stormwater_storage_capacity=2e4,
stormwater_storage_area=2e4,
treatment_throughput_capacity=5e4,
name="oxford_wwtw",
)
Sewers¶
Sewer nodes enable water to transition between households and WWTWs, and between impervious surfaces and rivers or WWTWs. They use a timearea diagram to represent travel time, which assigns a specified percentage of water to take a specified duration to pass through the sewer node.
combined_sewer = Sewer(
capacity=4e6, pipe_timearea={0: 0.8, 1: 0.15, 2: 0.05}, name="combined_sewer"
)
Groundwater¶
Groundwater nodes implement a simple residence time to determine baseflow
gw = Groundwater(capacity=3.2e9, area=3.2e8, name="gw", residence_time=20)
Create a nodelist¶
To keep all the nodes in one place, we put them into a list
nodelist = [
thames_above_abingdon,
evenlode,
thames,
ray,
cherwell,
oxford,
distribution,
farmoor,
oxford_fwtw,
oxford_wwtw,
combined_sewer,
oxford_land,
gw,
farmoor_abstraction,
evenlode_thames,
cherwell_ray,
cherwell_thames,
thames_mixer,
]
print(nodelist)
[<wsimod.nodes.waste.Waste object at 0x7f3b184dde90>, <wsimod.nodes.catchment.Catchment object at 0x7f3abe0dc610>, <wsimod.nodes.catchment.Catchment object at 0x7f3abe0c7610>, <wsimod.nodes.catchment.Catchment object at 0x7f3abe281e50>, <wsimod.nodes.catchment.Catchment object at 0x7f3abe280f50>, <wsimod.nodes.demand.ResidentialDemand object at 0x7f3abe3b5850>, <wsimod.nodes.nodes.Node object at 0x7f3abe3b7290>, <wsimod.nodes.storage.Reservoir object at 0x7f3abe3b7690>, <wsimod.nodes.wtw.FWTW object at 0x7f3abe0dcb10>, <wsimod.nodes.wtw.WWTW object at 0x7f3abe36f390>, <wsimod.nodes.sewer.Sewer object at 0x7f3abe36da90>, <wsimod.nodes.land.Land object at 0x7f3abe1c99d0>, <wsimod.nodes.storage.Groundwater object at 0x7f3abe36c190>, <wsimod.nodes.nodes.Node object at 0x7f3abe14de10>, <wsimod.nodes.nodes.Node object at 0x7f3abe14d690>, <wsimod.nodes.nodes.Node object at 0x7f3abe14c950>, <wsimod.nodes.nodes.Node object at 0x7f3abe123010>, <wsimod.nodes.nodes.Node object at 0x7f3abe122290>]
# Standard simple arcs
fwtw_to_distribution = Arc(
in_port=oxford_fwtw, out_port=distribution, name="fwtw_to_distribution"
)
print(fwtw_to_distribution)
<wsimod.arcs.arcs.Arc object at 0x7f3abe35b410>
As with nodes, even though we only gave it a few parameters, the arc comes with a lot built in
print(dir(fwtw_to_distribution))
['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'arc_mass_balance', 'blend_vqip', 'capacity', 'compare_vqip', 'concentration_to_total', 'copy_vqip', 'ds_vqip', 'ds_vqip_c', 'empty_vqip', 'empty_vqip_predefined', 'end_timestep', 'extract_vqip', 'extract_vqip_c', 'flow_in', 'flow_out', 'get_excess', 'in_port', 'mass_balance', 'mass_balance_ds', 'mass_balance_in', 'mass_balance_out', 'name', 'out_port', 'preference', 'reinit', 'send_pull_check', 'send_pull_request', 'send_push_check', 'send_push_request', 'sum_vqip', 'total_to_concentration', 'v_change_vqip', 'v_change_vqip_c', 'v_distill_vqip', 'v_distill_vqip_c', 'vqip_in', 'vqip_out']
We can see that the arc links the two nodes
print(fwtw_to_distribution.in_port)
<wsimod.nodes.wtw.FWTW object at 0x7f3abe0dcb10>
print(fwtw_to_distribution.out_port)
<wsimod.nodes.nodes.Node object at 0x7f3abe3b7290>
And that it has updated the nodes that it is connecting.
print(oxford_fwtw.out_arcs)
{'fwtw_to_distribution': <wsimod.arcs.arcs.Arc object at 0x7f3abe35b410>}
print(distribution.in_arcs)
{'fwtw_to_distribution': <wsimod.arcs.arcs.Arc object at 0x7f3abe35b410>}
We use arcs to send checks and requests..
print(fwtw_to_distribution.send_pull_check({"volume": 20}))
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 15.0}
reply = fwtw_to_distribution.send_pull_request({"volume": 20})
print(reply)
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 15.0}
They convey this information to the nodes that they connect to, which update their state variables
print(oxford_fwtw.service_reservoir_tank.storage)
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 0.0}
In turn, the arcs update their own state variables.
print(fwtw_to_distribution.flow_in)
print(fwtw_to_distribution.flow_out)
15.0 15.0
Arc parameters¶
Besides the in/out ports and names, arcs can have a parameter for their capacity, to limit the flow that may pass through it each timestep. A typical example would be on river abstractions to a reservoir
abstraction_to_farmoor = Arc(
in_port=farmoor_abstraction,
out_port=farmoor,
name="abstraction_to_farmoor",
capacity=5e4,
)
A bit more sophisticated is the 'preference' parameter. We use preference to express where we would prefer the model to send water. In this example, the sewer can send water to both the treatment plant and directly into the river. Of course we would always to prefer to send water to the plant, so we give it a very high preference. Discharging into the river should only be done if there is no capacity left at the WWTW, so we give the arc a very low preference.
sewer_to_wwtw = Arc(
in_port=combined_sewer, out_port=oxford_wwtw, preference=1e10, name="sewer_to_wwtw"
)
sewer_overflow = Arc(
in_port=combined_sewer,
out_port=thames_mixer,
preference=1e-10,
name="sewer_overflow",
)
Create arcs¶
Arcs are a bit less interesting than nodes because they generally don't capture complicated physical behaviours. So we just create all of them below.
evenlode_to_thames = Arc(
in_port=evenlode, out_port=evenlode_thames, name="evenlode_to_thames"
)
thames_to_thames = Arc(
in_port=thames, out_port=evenlode_thames, name="thames_to_thames"
)
ray_to_cherwell = Arc(in_port=ray, out_port=cherwell_ray, name="ray_to_cherwell")
cherwell_to_cherwell = Arc(
in_port=cherwell, out_port=cherwell_ray, name="cherwell_to_cherwell"
)
thames_to_farmoor = Arc(
in_port=evenlode_thames, out_port=farmoor_abstraction, name="thames_to_farmoor"
)
farmoor_to_mixer = Arc(
in_port=farmoor_abstraction, out_port=thames_mixer, name="farmoor_to_mixer"
)
cherwell_to_mixer = Arc(
in_port=cherwell_ray, out_port=thames_mixer, name="cherwell_to_mixer"
)
wwtw_to_mixer = Arc(in_port=oxford_wwtw, out_port=thames_mixer, name="wwtw_to_mixer")
mixer_to_waste = Arc(
in_port=thames_mixer, out_port=thames_above_abingdon, name="mixer_to_waste"
)
distribution_to_demand = Arc(
in_port=distribution, out_port=oxford, name="distribution_to_demand"
)
reservoir_to_fwtw = Arc(in_port=farmoor, out_port=oxford_fwtw, name="reservoir_to_fwtw")
fwtw_to_sewer = Arc(in_port=oxford_fwtw, out_port=combined_sewer, name="fwtw_to_sewer")
demand_to_sewer = Arc(in_port=oxford, out_port=combined_sewer, name="demand_to_sewer")
land_to_sewer = Arc(in_port=oxford_land, out_port=combined_sewer, name="land_to_sewer")
land_to_gw = Arc(in_port=oxford_land, out_port=gw, name="land_to_gw")
garden_to_gw = Arc(in_port=oxford, out_port=gw, name="garden_to_gw")
gw_to_mixer = Arc(in_port=gw, out_port=thames_mixer, name="gw_to_mixer")
Again, we keep all the arcs in a tidy list together.
arclist = [
evenlode_to_thames,
thames_to_thames,
ray_to_cherwell,
cherwell_to_cherwell,
thames_to_farmoor,
farmoor_to_mixer,
cherwell_to_mixer,
wwtw_to_mixer,
sewer_overflow,
mixer_to_waste,
abstraction_to_farmoor,
distribution_to_demand,
demand_to_sewer,
land_to_sewer,
sewer_to_wwtw,
fwtw_to_sewer,
fwtw_to_distribution,
reservoir_to_fwtw,
land_to_gw,
garden_to_gw,
gw_to_mixer,
]
Mapping¶
Remember, WSIMOD is an integrated model. Because it covers so many different things, it is very easy to make mistakes. Thus it is always good practice to plot your data!
Below we load the node location data and create arcs from the information in the arclist.
location_fn = os.path.join(data_folder, "raw", "points_locations.geojson")
nodes_gdf = gpd.read_file(location_fn).set_index("name")
arcs_gdf = []
for arc in arclist:
arcs_gdf.append(
{
"name": arc.name,
"geometry": LineString(
[
nodes_gdf.loc[arc.in_port.name, "geometry"],
nodes_gdf.loc[arc.out_port.name, "geometry"],
]
),
}
)
arcs_gdf = gpd.GeoDataFrame(arcs_gdf, crs=nodes_gdf.crs)
Because we converted the information as GeoDataFrames, we can simply plot them below
f, ax = plt.subplots()
arcs_gdf.plot(ax=ax)
nodes_gdf.plot(color="r", ax=ax, zorder=10)
<Axes: >
Orchestration¶
Orchestration is making the simulation happen by calling functions in the nodes. These functions simulate physical behaviour within the node, and cause pulls/pushes to happen which in turn triggers physical behaviour in other nodes.
Orchestrating an individual timestep¶
We will start below by manually orchestrating a single timestep.
We start by setting the date, so that every node knows what forcing data to read for this timestep.
date = dates[0]
for node in nodelist:
node.t = date
print(date)
print(oxford_fwtw.t)
2009-03-03 00:00:00 2009-03-03 00:00:00
We can see the service reservoirs are empty but the supply reservoir is not!
print(oxford_fwtw.service_reservoir_tank.storage)
print(farmoor.tank.storage)
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 0.0} {'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 10000000.0}
If we call the FWTW's treat_water function it will pull water from the supply reservoir and update its service reservoirs
oxford_fwtw.treat_water()
print(oxford_fwtw.service_reservoir_tank.storage)
print(farmoor.tank.storage)
{'temperature': 0.0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 43641.0} {'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 9955000.0}
This information is tracked in the arcs that enter the FWTW
print(oxford_fwtw.in_arcs)
{'reservoir_to_fwtw': <wsimod.arcs.arcs.Arc object at 0x7f3abe45d810>}
print(reservoir_to_fwtw.flow_in)
print(reservoir_to_fwtw.flow_out)
45000.0 45000.0
Although none of that water has yet entered the distribution network (only some small flow from the earlier demonstration)
print(fwtw_to_distribution.flow_in)
15.0
That is because no water consumption demand had yet been generated.
If we call the demand node's create_demand function we see that the distribution arc becomes utilised.
oxford.create_demand()
print(fwtw_to_distribution.flow_in)
30015.0
We also see that this gets pushed onwards into the sewer system
print(demand_to_sewer.flow_in)
30000.0
Many nodes have functions intended to be called during orchestration. These functions are described in the documentation. For example, we see in the Land node API reference that the 'run' function is intended to be called from orchestration.
oxford_land.run()
Below we call the functions for other nodes
# Discharge GW
gw.distribute()
# Discharge sewers (pushed to other sewers or WWTW)
combined_sewer.make_discharge()
# Run WWTW model
oxford_wwtw.calculate_discharge()
# Make abstractions
farmoor.make_abstractions()
# Discharge WW
oxford_wwtw.make_discharge()
# Route catchments
evenlode.route()
thames.route()
ray.route()
cherwell.route()
Ending a timestep¶
Because mistakes happen, it is essential to carry out mass balance testing. Each node has a mass balance function that can be called. We see a mass balance violation resulting from the demonstration with the FWTW earlier.
for node in nodelist:
in_, ds_, out_ = node.node_mass_balance()
mass balance error for volume of 15.0 in oxford_distribution mass balance error for volume of -15.0 in oxford_fwtw
We should also call the end_timestep function in nodes and arcs. This is important for mass balance testing and capturing the behaviour of some dynamic processes in nodes.
for node in nodelist:
node.end_timestep()
for arc in arclist:
arc.end_timestep()
Model object¶
Of course it would be a massive pain to manually orchestrate every timestep. So instead we store node and arc information in a model object that will do the orchestration for us.
Because we have already created the nodes/arcs above, we simply need to add the instantiated lists above.
my_model = Model()
my_model.add_instantiated_nodes(nodelist)
my_model.add_instantiated_arcs(arclist)
my_model.dates = dates
The model object lets us reinitialise the nodes/arcs, and run all of the orchestration with a 'run' function.
my_model.reinit()
flows, _, _, _ = my_model.run()
0%| | 0/1456 [00:00<?, ?it/s]
3%|█████▊ | 39/1456 [00:00<00:03, 386.91it/s]
5%|███████████▋ | 78/1456 [00:00<00:04, 334.74it/s]
8%|████████████████▌ | 112/1456 [00:00<00:04, 327.72it/s]
10%|██████████████████████ | 149/1456 [00:00<00:03, 339.81it/s]
13%|███████████████████████████▌ | 186/1456 [00:00<00:03, 348.65it/s]
16%|█████████████████████████████████▌ | 226/1456 [00:00<00:03, 363.43it/s]
18%|███████████████████████████████████████▏ | 264/1456 [00:00<00:03, 368.54it/s]
21%|████████████████████████████████████████████▉ | 303/1456 [00:00<00:03, 374.45it/s]
23%|██████████████████████████████████████████████████▌ | 341/1456 [00:00<00:03, 365.58it/s]
26%|████████████████████████████████████████████████████████▋ | 382/1456 [00:01<00:02, 376.28it/s]
29%|██████████████████████████████████████████████████████████████▎ | 420/1456 [00:01<00:02, 365.23it/s]
31%|███████████████████████████████████████████████████████████████████▊ | 457/1456 [00:01<00:02, 336.47it/s]
34%|████████████████████████████████████████████████████████████████████████▉ | 492/1456 [00:01<00:03, 308.10it/s]
36%|█████████████████████████████████████████████████████████████████████████████▋ | 524/1456 [00:01<00:03, 293.34it/s]
38%|██████████████████████████████████████████████████████████████████████████████████▏ | 554/1456 [00:01<00:03, 288.56it/s]
40%|██████████████████████████████████████████████████████████████████████████████████████▋ | 584/1456 [00:01<00:03, 284.74it/s]
42%|██████████████████████████████████████████████████████████████████████████████████████████▉ | 613/1456 [00:01<00:02, 282.63it/s]
44%|███████████████████████████████████████████████████████████████████████████████████████████████▋ | 645/1456 [00:01<00:02, 291.91it/s]
47%|█████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 685/1456 [00:02<00:02, 322.46it/s]
50%|███████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 726/1456 [00:02<00:02, 347.38it/s]
53%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 767/1456 [00:02<00:01, 363.46it/s]
55%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 804/1456 [00:02<00:02, 282.70it/s]
58%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 845/1456 [00:02<00:01, 313.77it/s]
61%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 886/1456 [00:02<00:01, 338.04it/s]
63%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 923/1456 [00:02<00:01, 340.63it/s]
66%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 961/1456 [00:02<00:01, 350.23it/s]
69%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 1003/1456 [00:02<00:01, 367.95it/s]
71%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 1041/1456 [00:03<00:01, 301.01it/s]
74%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 1081/1456 [00:03<00:01, 325.18it/s]
77%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 1123/1456 [00:03<00:00, 348.70it/s]
80%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 1164/1456 [00:03<00:00, 365.15it/s]
83%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 1206/1456 [00:03<00:00, 379.54it/s]
86%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 1249/1456 [00:03<00:00, 391.29it/s]
89%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 1289/1456 [00:03<00:00, 336.06it/s]
91%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 1325/1456 [00:03<00:00, 303.27it/s]
93%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 1358/1456 [00:04<00:00, 302.68it/s]
96%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 1395/1456 [00:04<00:00, 318.03it/s]
98%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 1428/1456 [00:04<00:00, 312.30it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1456/1456 [00:04<00:00, 330.97it/s]
The model outputs flows as a dictionary which can be converted to a dataframe
flows = pd.DataFrame(flows)
print(flows.sample(10))
arc flow time temperature \ 13882 thames_to_thames 4.233600e+05 2010-12-24 0.914286 13029 mixer_to_waste 1.796833e+06 2010-11-13 7.089646 26078 reservoir_to_fwtw 3.093421e+04 2012-07-26 11.338005 14037 mixer_to_waste 3.759096e+06 2010-12-31 2.729300 7066 abstraction_to_farmoor 3.093421e+04 2010-02-02 5.419641 14855 sewer_overflow 4.419173e-06 2011-02-08 11.528663 18863 farmoor_to_mixer 1.153410e+05 2011-08-18 17.134427 23335 thames_to_farmoor 6.998400e+05 2012-03-18 9.073898 4585 wwtw_to_mixer 4.998918e+04 2009-10-07 15.648625 29392 land_to_sewer 3.400000e+04 2012-12-31 7.716667 silicon fluoride chloride nitrate sulphate \ 13882 3.314631e+03 7.197120e+01 2.145619e+04 2.880855e+03 24371.928000 13029 1.493185e+04 1.969042e+02 7.191189e+04 1.049136e+04 106483.517749 26078 1.787156e+02 4.219087e+00 1.069621e+03 1.689888e+02 1591.591106 14037 2.907590e+04 7.016430e+02 1.631560e+05 2.839105e+04 225507.406041 7066 1.895541e+02 3.712106e+00 7.110377e+02 2.268185e+02 1240.169284 14855 1.114294e-07 2.228038e-09 8.992449e-07 2.809677e-07 0.000001 18863 7.875501e+02 1.483179e+01 5.235773e+03 4.830787e+02 7326.739583 23335 1.692908e+03 9.808869e+01 2.531403e+04 4.251599e+03 39288.931200 4585 1.938479e+01 4.237405e-01 1.664853e+02 5.293940e+01 242.989656 29392 1.400000e-01 4.000000e-01 1.200000e-01 4.000000e-02 140.000000 nitrogen sodium potassium calcium magnesium \ 13882 2.977430e+03 1.227139e+04 2.162160e+03 45731.952000 2.183328e+03 13029 1.254961e+04 4.720654e+04 1.033646e+04 178708.039241 9.899578e+03 26078 1.826501e+02 6.969689e+02 1.401725e+02 3094.635873 1.470924e+02 14037 2.949883e+04 9.504294e+04 1.748362e+04 404921.039208 2.026835e+04 7066 2.397946e+02 4.180010e+02 8.230724e+01 3374.696467 1.343434e+02 14855 2.394971e-07 9.366397e-07 1.447647e-07 0.000001 1.469513e-07 18863 5.346622e+02 3.937609e+03 7.091529e+02 11419.821462 6.270838e+02 23335 4.749766e+03 1.589797e+04 2.995611e+03 76657.412571 3.591648e+03 4585 4.525273e+01 1.782782e+02 2.835228e+01 180.602250 2.840135e+01 29392 8.000000e+00 6.000000e-01 1.400000e+01 140.000000 1.200000e+01 boron 13882 1.977696e+01 13029 1.099404e+02 26078 1.509684e+00 14037 1.890550e+02 7066 1.063915e+00 14855 2.328990e-09 18863 6.690212e+00 23335 3.221757e+01 4585 4.620017e-01 29392 2.000000e-01
Validation plots¶
Of course we shouldn't trust a model without proof, and this is doubly true for an integrated model whose errors may easily propagate.
Thankfully we have rich in-river water quality sampling in Oxford that we can use for validation
We load and format this data for dates and pollutants that overlap with what we have simulated
mixer_val_df = pd.read_csv(os.path.join(data_folder, "raw", "mixer_wims_val.csv"))
mixer_val_df.date = pd.to_datetime(mixer_val_df.date)
mixer_val_df = mixer_val_df.loc[mixer_val_df.date.isin(dates)]
val_pollutants = mixer_val_df.variable.unique()
mixer_val_df = mixer_val_df.pivot(index="date", columns="variable", values="value")
We get rid of the first day of flows because our tanks were initialised empty and will not be informative
flows = flows.loc[flows.time != dates[0]]
We convert our flows, which are simulated in kg/d into a concentration value (kg/m3/d).
flows_plot = flows.copy()
# Convert to KG/M3
for pol in set(val_pollutants):
if pol != "temperature":
flows_plot[pol] /= flows_plot.flow
We pick the model outlet as a validation location and make a pretty plot - fantastic!
plot_arc = "mixer_to_waste"
f, axs = plt.subplots(val_pollutants.size, 1)
for pol, ax in zip(val_pollutants, axs):
ax.plot(
flows_plot.loc[flows_plot.arc == plot_arc, [pol, "time"]].set_index("time"),
color="b",
label="simulation",
)
ax.plot(mixer_val_df[pol], ls="", marker="o", color="r", label="spot sample")
ax.set_ylabel(pol)
plt.legend()
<matplotlib.legend.Legend at 0x7f3abe4478d0>