Coverage for wsimod/demo/create_oxford.py: 0%
181 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-01-11 16:39 +0000
1# -*- coding: utf-8 -*-
2"""Created on Tue Nov 1 09:03:17 2022.
4@author: Barney
5"""
6import os
7import re
9import pandas as pd
11from wsimod.arcs.arcs import Arc
12from wsimod.core import constants
13from wsimod.nodes.catchment import Catchment
14from wsimod.nodes.demand import ResidentialDemand
15from wsimod.nodes.land import Land
16from wsimod.nodes.nodes import Node
17from wsimod.nodes.sewer import Sewer
18from wsimod.nodes.storage import Groundwater, Reservoir
19from wsimod.nodes.waste import Waste
20from wsimod.nodes.wtw import FWTW, WWTW
21from wsimod.orchestration.model import Model
23# NOTE - these are only supplements to the demos in the documentation
24# if you want to see demos, go to docs/demo/scripts!
27def preprocess_oxford(data_dir=None):
28 """
30 Args:
31 data_dir:
32 """
33 # This function is used to create the timeseries data for the oxford demo
34 if data_dir is None:
35 data_dir = os.path.join(
36 os.path.dirname(
37 os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
38 ),
39 "demo",
40 "data",
41 )
43 UG_L_TO_KG_M3 = 1e-6
44 MG_L_TO_KG_M3 = 1e-3
46 q_lab = [
47 ("39008_gdf.csv", "thames"),
48 ("39034_gdf.csv", "evenlode"),
49 ("39021_gdf.csv", "cherwell"),
50 ("39140_gdf.csv", "ray"),
51 ]
53 flows = []
54 for fn, river in q_lab:
55 df = pd.read_csv(
56 os.path.join(data_dir, "raw", fn), skiprows=19, error_bad_lines=False
57 )
58 if df.shape[1] == 2:
59 df.columns = ["date", "flow"]
60 else:
61 df.columns = ["date", "flow", "missing"]
62 df["site"] = river
63 flows.append(df)
65 flows = pd.concat(flows)
66 flows = flows.pivot(columns="site", index="date", values="flow")
67 flows.index = pd.to_datetime(flows.index)
69 wq_data = pd.read_csv(
70 os.path.join(
71 data_dir, "raw", "CEHThamesInitiative_WaterQualityData_2009-2013.csv"
72 )
73 )
75 rain = pd.read_csv(
76 os.path.join(data_dir, "raw", "39008_cdr.csv"),
77 skiprows=19,
78 error_bad_lines=False,
79 )
80 rain.columns = ["date", "value", "misc"]
81 rain["site"] = "oxford_land"
82 rain["variable"] = "precipitation"
83 rain = rain.drop("misc", axis=1)
84 rain.date = pd.to_datetime(rain.date)
86 sites = {
87 "River Thames at Swinford": "thames",
88 "River Evenlode at Cassington Mill": "evenlode",
89 "River Ray at Islip": "ray",
90 "River Cherwell at Hampton Poyle": "cherwell",
91 }
92 wq = wq_data.loc[wq_data.Site.isin(sites.keys())]
93 wq["Sampling date (dd/mm/yyyy)"] = pd.to_datetime(
94 wq["Sampling date (dd/mm/yyyy)"], format="%d/%m/%Y"
95 )
96 wq = wq.rename(columns={"Sampling date (dd/mm/yyyy)": "date", "Site": "site"})
98 wq.site = [sites[x] for x in wq.site]
100 wq = wq.set_index("date").drop("Sampling time (hh:mm)", axis=1)
101 wq[wq.columns.drop("site")] = wq[wq.columns.drop("site")].apply(
102 lambda x: pd.to_numeric(x, errors="coerce")
103 )
104 wq = wq.dropna(axis=1, how="any")
105 wq = wq.drop("Mean daily river discharge (m3 s-1)", axis=1)
106 wq = wq.groupby("site").resample("d").interpolate().drop("site", axis=1)
107 wq = wq.reset_index()
108 wq.loc[:, wq.columns.str.contains("ug")] *= UG_L_TO_KG_M3
109 wq.loc[:, wq.columns.str.contains("mg")] *= MG_L_TO_KG_M3
111 columns = []
112 for pol in wq.columns.unique():
113 text = pol.lower()
114 for sub in ["water", "dissolved", "total", " ", r"\(.*\)"]:
115 text = re.sub(sub, "", text)
116 columns.append(text)
117 wq.columns = columns
119 # Convert to nitrate as N
120 wq["nitrate"] /= 4.43
122 # Convert to Silica as SiO2
123 wq["silicon"] *= 2.14
125 wq = wq.melt(id_vars=["site", "date"])
127 # wq.date = pd.to_datetime(wq.date.dt.date)
129 flows = flows.loc[flows.index.isin(wq.date)]
130 flows = flows.unstack().rename("value").reset_index()
131 flows["variable"] = "flow"
133 rain = rain.loc[rain.date.isin(wq.date)]
134 evaporation = create_timeseries(2 / 1000, rain.date, "et0")
135 evaporation["site"] = "oxford_land"
137 temperature_ = (
138 wq.loc[wq.variable == "temperature"].groupby("date").mean().reset_index()
139 )
140 temperature_["site"] = "oxford_land"
141 temperature_["variable"] = "temperature"
143 input_data = pd.concat([wq, flows, rain, evaporation, temperature_], axis=0)
144 input_data.to_csv(
145 os.path.join(data_dir, "processed", "timeseries_data.csv"), index=False
146 )
149def create_oxford_model(data_folder):
150 """
152 Args:
153 data_folder:
155 Returns:
157 """
158 input_fid = os.path.join(data_folder, "processed", "timeseries_data.csv")
159 input_data = pd.read_csv(input_fid)
160 input_data.loc[input_data.variable == "flow", "value"] *= constants.M3_S_TO_M3_DT
161 input_data.loc[input_data.variable == "precipitation", "value"] *= constants.MM_TO_M
162 input_data.date = pd.to_datetime(input_data.date)
163 data_input_dict = input_data.set_index(["variable", "date"]).value.to_dict()
164 data_input_dict = (
165 input_data.groupby("site")
166 .apply(lambda x: x.set_index(["variable", "date"]).value.to_dict())
167 .to_dict()
168 )
169 dates = input_data.date.unique()
170 dates = dates[dates.argsort()]
171 dates = [pd.Timestamp(x) for x in dates]
172 constants.POLLUTANTS = input_data.variable.unique().tolist()
173 constants.POLLUTANTS.remove("flow")
174 constants.POLLUTANTS.remove("precipitation")
175 constants.POLLUTANTS.remove("et0")
176 constants.NON_ADDITIVE_POLLUTANTS = ["temperature"]
177 constants.ADDITIVE_POLLUTANTS = list(
178 set(constants.POLLUTANTS).difference(constants.NON_ADDITIVE_POLLUTANTS)
179 )
180 constants.FLOAT_ACCURACY = 1e-8
181 thames_above_abingdon = Waste(name="thames_above_abingdon")
182 farmoor_abstraction = Node(name="farmoor_abstraction")
183 evenlode_thames = Node(name="evenlode_thames")
184 cherwell_ray = Node(name="cherwell_ray")
185 cherwell_thames = Node(name="cherwell_thames")
186 thames_mixer = Node(name="thames_mixer")
187 evenlode = Catchment(name="evenlode", data_input_dict=data_input_dict["evenlode"])
188 thames = Catchment(name="thames", data_input_dict=data_input_dict["thames"])
189 ray = Catchment(name="ray", data_input_dict=data_input_dict["ray"])
190 cherwell = Catchment(name="cherwell", data_input_dict=data_input_dict["cherwell"])
191 oxford_fwtw = FWTW(
192 service_reservoir_storage_capacity=1e5,
193 service_reservoir_storage_area=2e4,
194 service_reservoir_initial_storage=0.9e5,
195 treatment_throughput_capacity=4.5e4,
196 name="oxford_fwtw",
197 )
198 land_inputs = data_input_dict["oxford_land"]
199 pollutant_deposition = {
200 "boron": 100e-10,
201 "calcium": 70e-7,
202 "chloride": 60e-10,
203 "fluoride": 0.2e-7,
204 "magnesium": 6e-7,
205 "nitrate": 2e-9,
206 "nitrogen": 4e-7,
207 "potassium": 7e-7,
208 "silicon": 7e-9,
209 "sodium": 30e-9,
210 "sulphate": 70e-7,
211 }
212 surface = [
213 {
214 "type_": "PerviousSurface",
215 "area": 2e7,
216 "pollutant_load": pollutant_deposition,
217 "surface": "rural",
218 "field_capacity": 0.3,
219 "depth": 0.5,
220 "initial_storage": 2e7 * 0.4 * 0.5,
221 },
222 {
223 "type_": "ImperviousSurface",
224 "area": 1e7,
225 "pollutant_load": pollutant_deposition,
226 "surface": "urban",
227 "initial_storage": 5e6,
228 },
229 ]
230 oxford_land = Land(
231 surfaces=surface, name="oxford_land", data_input_dict=land_inputs
232 )
233 oxford = ResidentialDemand(
234 name="oxford",
235 population=2e5,
236 per_capita=0.15,
237 pollutant_load={
238 "boron": 500 * constants.UG_L_TO_KG_M3 * 0.15,
239 "calcium": 150 * constants.MG_L_TO_KG_M3 * 0.15,
240 "chloride": 180 * constants.MG_L_TO_KG_M3 * 0.15,
241 "fluoride": 0.4 * constants.MG_L_TO_KG_M3 * 0.15,
242 "magnesium": 30 * constants.MG_L_TO_KG_M3 * 0.15,
243 "nitrate": 60 * constants.MG_L_TO_KG_M3 * 0.15,
244 "nitrogen": 50 * constants.MG_L_TO_KG_M3 * 0.15,
245 "potassium": 30 * constants.MG_L_TO_KG_M3 * 0.15,
246 "silicon": 20 * constants.MG_L_TO_KG_M3 * 0.15,
247 "sodium": 200 * constants.MG_L_TO_KG_M3 * 0.15,
248 "sulphate": 250 * constants.MG_L_TO_KG_M3 * 0.15,
249 "temperature": 14,
250 },
251 data_input_dict=land_inputs,
252 )
253 farmoor = Reservoir(
254 name="farmoor", capacity=1e7, initial_storage=1e7, area=1.5e6, datum=62
255 )
256 distribution = Node(name="oxford_distribution")
257 oxford_wwtw = WWTW(
258 stormwater_storage_capacity=2e4,
259 stormwater_storage_area=2e4,
260 treatment_throughput_capacity=5e4,
261 name="oxford_wwtw",
262 )
263 combined_sewer = Sewer(
264 capacity=4e6, pipe_timearea={0: 0.8, 1: 0.15, 2: 0.05}, name="combined_sewer"
265 )
266 gw = Groundwater(capacity=3.2e9, area=3.2e8, name="gw", residence_time=20)
267 nodelist = [
268 thames_above_abingdon,
269 evenlode,
270 thames,
271 ray,
272 cherwell,
273 oxford,
274 distribution,
275 farmoor,
276 oxford_fwtw,
277 oxford_wwtw,
278 combined_sewer,
279 oxford_land,
280 gw,
281 farmoor_abstraction,
282 evenlode_thames,
283 cherwell_ray,
284 cherwell_thames,
285 thames_mixer,
286 ]
287 fwtw_to_distribution = Arc(
288 in_port=oxford_fwtw, out_port=distribution, name="fwtw_to_distribution"
289 )
290 abstraction_to_farmoor = Arc(
291 in_port=farmoor_abstraction,
292 out_port=farmoor,
293 name="abstraction_to_farmoor",
294 capacity=5e4,
295 )
296 sewer_to_wwtw = Arc(
297 in_port=combined_sewer,
298 out_port=oxford_wwtw,
299 preference=1e10,
300 name="sewer_to_wwtw",
301 )
302 sewer_overflow = Arc(
303 in_port=combined_sewer,
304 out_port=thames_mixer,
305 preference=1e-10,
306 name="sewer_overflow",
307 )
308 evenlode_to_thames = Arc(
309 in_port=evenlode, out_port=evenlode_thames, name="evenlode_to_thames"
310 )
312 thames_to_thames = Arc(
313 in_port=thames, out_port=evenlode_thames, name="thames_to_thames"
314 )
316 ray_to_cherwell = Arc(in_port=ray, out_port=cherwell_ray, name="ray_to_cherwell")
318 cherwell_to_cherwell = Arc(
319 in_port=cherwell, out_port=cherwell_ray, name="cherwell_to_cherwell"
320 )
322 thames_to_farmoor = Arc(
323 in_port=evenlode_thames, out_port=farmoor_abstraction, name="thames_to_farmoor"
324 )
326 farmoor_to_mixer = Arc(
327 in_port=farmoor_abstraction, out_port=thames_mixer, name="farmoor_to_mixer"
328 )
330 cherwell_to_mixer = Arc(
331 in_port=cherwell_ray, out_port=thames_mixer, name="cherwell_to_mixer"
332 )
334 wwtw_to_mixer = Arc(
335 in_port=oxford_wwtw, out_port=thames_mixer, name="wwtw_to_mixer"
336 )
338 mixer_to_waste = Arc(
339 in_port=thames_mixer, out_port=thames_above_abingdon, name="mixer_to_waste"
340 )
342 distribution_to_demand = Arc(
343 in_port=distribution, out_port=oxford, name="distribution_to_demand"
344 )
346 reservoir_to_fwtw = Arc(
347 in_port=farmoor, out_port=oxford_fwtw, name="reservoir_to_fwtw"
348 )
350 fwtw_to_sewer = Arc(
351 in_port=oxford_fwtw, out_port=combined_sewer, name="fwtw_to_sewer"
352 )
354 demand_to_sewer = Arc(
355 in_port=oxford, out_port=combined_sewer, name="demand_to_sewer"
356 )
358 land_to_sewer = Arc(
359 in_port=oxford_land, out_port=combined_sewer, name="land_to_sewer"
360 )
362 land_to_gw = Arc(in_port=oxford_land, out_port=gw, name="land_to_gw")
364 garden_to_gw = Arc(in_port=oxford, out_port=gw, name="garden_to_gw")
366 gw_to_mixer = Arc(in_port=gw, out_port=thames_mixer, name="gw_to_mixer")
367 arclist = [
368 evenlode_to_thames,
369 thames_to_thames,
370 ray_to_cherwell,
371 cherwell_to_cherwell,
372 thames_to_farmoor,
373 farmoor_to_mixer,
374 cherwell_to_mixer,
375 wwtw_to_mixer,
376 sewer_overflow,
377 mixer_to_waste,
378 abstraction_to_farmoor,
379 distribution_to_demand,
380 demand_to_sewer,
381 land_to_sewer,
382 sewer_to_wwtw,
383 fwtw_to_sewer,
384 fwtw_to_distribution,
385 reservoir_to_fwtw,
386 land_to_gw,
387 garden_to_gw,
388 gw_to_mixer,
389 ]
390 my_model = Model()
391 my_model.add_instantiated_nodes(nodelist)
392 my_model.add_instantiated_arcs(arclist)
393 my_model.dates = dates
394 return my_model
397def pull_check_mrf(node, vqip=None):
398 """
400 Args:
401 node:
402 vqip:
404 Returns:
406 """
407 reply = node.pull_check_basic()
408 reply["volume"] = max(
409 reply["volume"] - (node.mrf - node.mrf_satisfied_this_timestep), 0
410 )
411 if vqip is not None:
412 reply["volume"] = min(reply["volume"], vqip["volume"])
413 return reply
416def pull_set_mrf(node, vqip):
417 """
419 Args:
420 node:
421 vqip:
423 Returns:
425 """
426 request = node.copy_vqip(vqip)
427 request["volume"] += node.mrf - node.mrf_satisfied_this_timestep
428 reply = node.pull_distributed(request)
429 reply_to_mrf = min((node.mrf - node.mrf_satisfied_this_timestep), reply["volume"])
430 node.mrf_satisfied_this_timestep += reply_to_mrf
431 reply_to_mrf = node.v_change_vqip(reply, reply_to_mrf)
432 reply = node.extract_vqip(reply, reply_to_mrf)
433 mrf_route_reply = node.push_distributed(reply_to_mrf, of_type=["Node"])
434 if mrf_route_reply["volume"] > constants.FLOAT_ACCURACY:
435 print("warning MRF not able to push")
436 return reply
439def push_set_mrf(node, vqip):
440 """
442 Args:
443 node:
444 vqip:
446 Returns:
448 """
449 reply = node.push_distributed(vqip)
450 total_pushed_downstream = vqip["volume"] - reply["volume"]
451 node.mrf_satisfied_this_timestep = min(
452 node.mrf, node.mrf_satisfied_this_timestep + total_pushed_downstream
453 )
454 return reply
457def end_timestep(node):
458 """
460 Args:
461 node:
462 """
463 node.mrf_satisfied_this_timestep = 0
466def convert_to_mrf(node, mrf=5):
467 """
469 Args:
470 node:
471 mrf:
472 """
473 node.mrf = mrf
474 node.mrf_satisfied_this_timestep = 0
475 node.pull_set_handler["default"] = lambda x: pull_set_mrf(node, x)
476 node.pull_check_handler["default"] = lambda x: pull_check_mrf(node, x)
477 node.end_timestep = lambda: end_timestep(node)
480def create_oxford_model_mrf(data_folder):
481 """
483 Args:
484 data_folder:
486 Returns:
488 """
489 my_model = create_oxford_model(data_folder)
490 convert_to_mrf(my_model.nodes["farmoor_abstraction"], mrf=3 * 86400)
491 return my_model
494def create_timeseries(amount, dates, variable):
495 """Create a timeseries with a constant value, formatted as a dataframe.
497 Args:
498 amount (float): Constant value to be applied over the timeseries
499 dates (list): list or iterable of dates to be generated
500 variable (str): String to store in the 'variable' column of the dataframe
502 Returns:
503 (DataFrame): the formatted dataframe
504 """
505 df = pd.DataFrame(columns=["date", "variable", "value"])
506 df["date"] = dates
507 df["variable"] = variable
508 df["value"] = amount
509 return df