diff --git a/dev-requirements.txt b/dev-requirements.txt index 218d771f..687ea873 100644 --- a/dev-requirements.txt +++ b/dev-requirements.txt @@ -4,6 +4,10 @@ # # pip-compile --extra=dev --output-file=dev-requirements.txt # +aenum==3.1.11 + # via + # pyswmm + # swmm-toolkit affine==2.4.0 # via # pysheds @@ -102,6 +106,8 @@ imageio==2.33.1 # via scikit-image iniconfig==2.0.0 # via pytest +julian==0.14 + # via pyswmm kiwisolver==1.4.5 # via matplotlib lazy-loader==0.3 @@ -162,6 +168,7 @@ packaging==23.2 # fastparquet # geopandas # matplotlib + # pyswmm # pytest # rioxarray # scikit-image @@ -203,6 +210,8 @@ pyproject-hooks==1.0.0 # via build pysheds==0.3.5 # via swmmanywhere (pyproject.toml) +pyswmm==1.5.1 + # via swmmanywhere (pyproject.toml) pytest==7.4.4 # via # pytest-cov @@ -264,6 +273,8 @@ snkit==1.9.0 # via swmmanywhere (pyproject.toml) snuggs==1.4.7 # via rasterio +swmm-toolkit==0.15.3 + # via pyswmm tifffile==2023.12.9 # via scikit-image tqdm==4.66.1 diff --git a/pyproject.toml b/pyproject.toml index e317ac9c..3024c537 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ dependencies = [ # TODO definitely don't need all of these "pandas", "pyarrow", "pysheds", + "pyswmm", "PyYAML", "rasterio", "rioxarray", diff --git a/requirements.txt b/requirements.txt index 120ed9d3..a65d2cac 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,6 +4,10 @@ # # pip-compile # +aenum==3.1.11 + # via + # pyswmm + # swmm-toolkit affine==2.4.0 # via # pysheds @@ -80,6 +84,8 @@ idna==3.6 # via requests imageio==2.33.1 # via scikit-image +julian==0.14 + # via pyswmm kiwisolver==1.4.5 # via matplotlib lazy-loader==0.3 @@ -131,6 +137,7 @@ packaging==23.2 # fastparquet # geopandas # matplotlib + # pyswmm # rioxarray # scikit-image # xarray @@ -161,6 +168,8 @@ pyproj==3.6.1 # rioxarray pysheds==0.3.5 # via swmmanywhere (pyproject.toml) +pyswmm==1.5.1 + # via swmmanywhere (pyproject.toml) python-dateutil==2.8.2 # via # matplotlib @@ -206,6 +215,8 @@ snkit==1.9.0 # via swmmanywhere (pyproject.toml) snuggs==1.4.7 # via rasterio +swmm-toolkit==0.15.3 + # via pyswmm tifffile==2023.12.9 # via scikit-image tqdm==4.66.1 diff --git a/swmmanywhere/defs/basic_drainage_all_bits.inp b/swmmanywhere/defs/basic_drainage_all_bits.inp new file mode 100644 index 00000000..0d8a4f5b --- /dev/null +++ b/swmmanywhere/defs/basic_drainage_all_bits.inp @@ -0,0 +1,173 @@ +[TITLE] +;;Project Title/Notes + +[OPTIONS] +;;Option Value +FLOW_UNITS LPS +INFILTRATION HORTON +FLOW_ROUTING KINWAVE +LINK_OFFSETS DEPTH +MIN_SLOPE 0 +ALLOW_PONDING NO +SKIP_STEADY_STATE NO + +START_DATE 10/03/2020 +START_TIME 00:00:00 +REPORT_START_DATE 10/03/2020 +REPORT_START_TIME 00:00:00 +END_DATE 10/08/2020 +END_TIME 00:00:00 +SWEEP_START 1/1 +SWEEP_END 12/31 +DRY_DAYS 0 +REPORT_STEP 00:15:00 +WET_STEP 00:05:00 +DRY_STEP 01:00:00 +ROUTING_STEP 0:00:05 +RULE_STEP 00:00:00 + +INERTIAL_DAMPING PARTIAL +NORMAL_FLOW_LIMITED BOTH +FORCE_MAIN_EQUATION H-W +VARIABLE_STEP 0.75 +LENGTHENING_STEP 0 +MIN_SURFAREA 0 +MAX_TRIALS 0 +HEAD_TOLERANCE 0 +SYS_FLOW_TOL 5 +LAT_FLOW_TOL 5 +MINIMUM_STEP 0.5 +THREADS 1 + +[EVAPORATION] +;;Data Source Parameters +;;-------------- ---------------- +CONSTANT 0.0 +DRY_ONLY NO + +[RAINGAGES] +;;Name Format Interval SCF Source +;;-------------- --------- ------ ------ ---------- +1 INTENSITY 00:05 1.0 FILE "january.dat" 1 MM + +[SUBCATCHMENTS] +;;Name Rain Gage Outlet Area %Imperv Width %Slope CurbLen SnowPack +;;-------------- ---------------- ---------------- -------- -------- -------- -------- -------- ---------------- +1 1 5 1.166 100 500 0.5 0 empty + +[SUBAREAS] +;;Subcatchment N-Imperv N-Perv S-Imperv S-Perv PctZero RouteTo PctRouted +;;-------------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- +1 0.01 0.1 0.05 0.05 25 OUTLET + +[INFILTRATION] +;;Subcatchment MaxRate MinRate Decay DryTime MaxInfil +;;-------------- ---------- ---------- ---------- ---------- ---------- +1 3.0 0.5 4 7 0 + +[SNOWPACKS] +;;Name Surface Parameters +;;-------------- ---------- ---------- +empty PLOWABLE 0.001 0.001 32.0 0.10 0.00 0.00 0.0 +empty IMPERVIOUS 0.001 0.001 32.0 0.10 0.00 0.00 0.00 +empty PERVIOUS 0.001 0.001 32.0 0.10 0.00 0.00 0.00 +empty REMOVAL 1.0 0.0 0.0 0.0 0.0 0.0 + +[JUNCTIONS] +;;Name Elevation MaxDepth InitDepth SurDepth Aponded +;;-------------- ---------- ---------- ---------- ---------- ---------- +4 30.85 0.5 0 0 0 +6 30.85 0 0 0 0 +7 30.85 0 0 0 0 + +[OUTFALLS] +;;Name Elevation Type Stage Data Gated Route To +;;-------------- ---------- ---------- ---------------- -------- ---------------- +3 29.759 FREE NO + +[STORAGE] +;;Name Elev. MaxDepth InitDepth Shape Curve Name/Params N/A Fevap Psi Ksat IMD +;;-------------- -------- ---------- ----------- ---------- ---------------------------- -------- -------- -------- -------- +5 30.85 0.5 0 FUNCTIONAL 1000 0 0 0 0 +8 30.85 1 0 TABULAR t1 0 0 + +[CONDUITS] +;;Name From Node To Node Length Roughness InOffset OutOffset InitFlow MaxFlow +;;-------------- ---------------- ---------------- ---------- ---------- ---------- ---------- ---------- ---------- +1 4 3 37.9 0.014344444444 0 0 0 0 +6 6 3 400 0.01 0 0 0 0 +7 7 3 400 0.01 0 0 0 0 +8 5 8 400 0.01 0 0 0 0 +9 8 7 400 0.01 0 0 0 0 + +[PUMPS] +;;Name From Node To Node Pump Curve Status Sartup Shutoff +;;-------------- ---------------- ---------------- ---------------- ------ -------- -------- +4 5 7 * ON 0 0 + +[ORIFICES] +;;Name From Node To Node Type Offset Qcoeff Gated CloseTime +;;-------------- ---------------- ---------------- ------------ ---------- ---------- -------- ---------- +3 5 6 SIDE 0 0.65 NO 0 + +[WEIRS] +;;Name From Node To Node Type CrestHt Qcoeff Gated EndCon EndCoeff Surcharge RoadWidth RoadSurf Coeff. Curve +;;-------------- ---------------- ---------------- ------------ ---------- ---------- -------- -------- ---------- ---------- ---------- ---------- ---------------- +2 5 4 TRANSVERSE 0.25 3.33 NO 0 0 YES + +[XSECTIONS] +;;Link Shape Geom1 Geom2 Geom3 Geom4 Barrels Culvert +;;-------------- ------------ ---------------- ---------- ---------- ---------- ---------- ---------- +1 TRIANGULAR 10.7 43.0 0.0 0.0 1 +6 CIRCULAR 1 0 0 0 1 +7 CIRCULAR 1 0 0 0 1 +8 RECT_CLOSED 1.5 0.5 0 0 1 +9 CIRCULAR 1 0 0 0 1 +3 CIRCULAR 1 0 0 0 +2 RECT_OPEN 1 2 0 0 + +[CURVES] +;;Name Type X-Value Y-Value +;;-------------- ---------- ---------- ---------- +t1 Storage 0 1 +t1 1 2 + +[REPORT] +;;Reporting Options +SUBCATCHMENTS ALL +NODES ALL +LINKS ALL + +[TAGS] + +[MAP] +DIMENSIONS 0.000 0.000 10000.000 10000.000 +Units None + +[COORDINATES] +;;Node X-Coord Y-Coord +;;-------------- ------------------ ------------------ +4 4743.590 5396.270 +6 4615.807 4829.857 +7 4484.083 4204.171 +3 6072.261 5477.855 +5 3298.368 5069.930 +8 3092.287 3801.653 + +[VERTICES] +;;Link X-Coord Y-Coord +;;-------------- ------------------ ------------------ + +[Polygons] +;;Subcatchment X-Coord Y-Coord +;;-------------- ------------------ ------------------ +1 -594.002 2399.077 +1 2750.865 5801.615 +1 1516.724 7347.174 +1 -836.217 5963.091 + +[SYMBOLS] +;;Gage X-Coord Y-Coord +;;-------------- ------------------ ------------------ +1 268.065 7715.618 + diff --git a/swmmanywhere/defs/storm.dat b/swmmanywhere/defs/storm.dat new file mode 100644 index 00000000..4b450941 --- /dev/null +++ b/swmmanywhere/defs/storm.dat @@ -0,0 +1,144 @@ +1 2020 10 03 00 00 0.06383868499999999 +1 2020 10 03 01 00 0.041554296 +1 2020 10 03 02 00 0.05283641 +1 2020 10 03 03 00 0.06305583499999999 +1 2020 10 03 04 00 0.0631444 +1 2020 10 03 05 00 0.025762885 +1 2020 10 03 06 00 0.015543459999999999 +1 2020 10 03 07 00 0.002143075 +1 2020 10 03 08 00 0.0005809343000000001 +1 2020 10 03 09 00 0.00029047078 +1 2020 10 03 10 00 7.275957999999999e-09 +1 2020 10 03 11 00 7.275957999999999e-09 +1 2020 10 03 12 00 7.275957999999999e-09 +1 2020 10 03 13 00 7.275957999999999e-09 +1 2020 10 03 14 00 7.275957999999999e-09 +1 2020 10 03 15 00 7.275957999999999e-09 +1 2020 10 03 16 00 7.275957999999999e-09 +1 2020 10 03 17 00 7.275957999999999e-09 +1 2020 10 03 18 00 0.013442892 +1 2020 10 03 19 00 0.05766097 +1 2020 10 03 20 00 0.040084262 +1 2020 10 03 21 00 0.0140769625 +1 2020 10 03 22 00 0.012433353 +1 2020 10 03 23 00 0.0035387275 +1 2020 10 04 00 00 0.006939299 +1 2020 10 04 01 00 0.004810397500000001 +1 2020 10 04 02 00 0.0047572685 +1 2020 10 04 03 00 0.007219139699999999 +1 2020 10 04 04 00 0.0034714249 +1 2020 10 04 05 00 0.0043428226 +1 2020 10 04 06 00 0.006737391200000001 +1 2020 10 04 07 00 7.275957999999999e-09 +1 2020 10 04 08 00 7.275957999999999e-09 +1 2020 10 04 09 00 7.275957999999999e-09 +1 2020 10 04 10 00 0.0031667878 +1 2020 10 04 11 00 0.011328164 +1 2020 10 04 12 00 0.024292849999999998 +1 2020 10 04 13 00 0.020297179 +1 2020 10 04 14 00 0.022971587999999998 +1 2020 10 04 15 00 0.023254965000000002 +1 2020 10 04 16 00 0.026450084999999998 +1 2020 10 04 17 00 0.02760841 +1 2020 10 04 18 00 0.03214959 +1 2020 10 04 19 00 0.13477959 +1 2020 10 04 20 00 0.08938194000000001 +1 2020 10 04 21 00 0.09367871 +1 2020 10 04 22 00 0.07757912 +1 2020 10 04 23 00 0.038086422 +1 2020 10 05 00 00 0.0710649 +1 2020 10 05 01 00 0.025242174 +1 2020 10 05 02 00 0.019326595 +1 2020 10 05 03 00 0.055535613 +1 2020 10 05 04 00 0.025847905 +1 2020 10 05 05 00 0.012900928 +1 2020 10 05 06 00 0.0010201766 +1 2020 10 05 07 00 0.00010981603 +1 2020 10 05 08 00 7.275957999999999e-09 +1 2020 10 05 09 00 7.275957999999999e-09 +1 2020 10 05 10 00 7.275957999999999e-09 +1 2020 10 05 11 00 0.0005809343000000001 +1 2020 10 05 12 00 0.0005809343000000001 +1 2020 10 05 13 00 0.0026885828 +1 2020 10 05 14 00 0.009766023 +1 2020 10 05 15 00 0.03863193 +1 2020 10 05 16 00 0.018175357 +1 2020 10 05 17 00 0.0071270406 +1 2020 10 05 18 00 0.00036840356 +1 2020 10 05 19 00 0.06358718399999999 +1 2020 10 05 20 00 0.035277407999999996 +1 2020 10 05 21 00 0.057979774 +1 2020 10 05 22 00 0.06011222 +1 2020 10 05 23 00 0.06900684 +1 2020 10 06 00 00 0.07822027 +1 2020 10 06 01 00 0.15135737999999999 +1 2020 10 06 02 00 0.107465195 +1 2020 10 06 03 00 0.08602388 +1 2020 10 06 04 00 0.06740928 +1 2020 10 06 05 00 0.0111085465 +1 2020 10 06 06 00 0.0130497065 +1 2020 10 06 07 00 0.0067621877 +1 2020 10 06 08 00 0.0020120133 +1 2020 10 06 09 00 0.0049308364 +1 2020 10 06 10 00 0.008625414999999999 +1 2020 10 06 11 00 0.010024611000000001 +1 2020 10 06 12 00 0.019184903000000003 +1 2020 10 06 13 00 0.015600133 +1 2020 10 06 14 00 0.014664976 +1 2020 10 06 15 00 0.028685259999999997 +1 2020 10 06 16 00 0.099735975 +1 2020 10 06 17 00 0.05485196 +1 2020 10 06 18 00 0.06388118999999999 +1 2020 10 06 19 00 0.08494349 +1 2020 10 06 20 00 0.23213516 +1 2020 10 06 21 00 0.09611933 +1 2020 10 06 22 00 0.14322790000000002 +1 2020 10 06 23 00 0.03817852 +1 2020 10 07 00 00 0.014767706 +1 2020 10 07 01 00 0.017080797999999998 +1 2020 10 07 02 00 0.012560871 +1 2020 10 07 03 00 0.020651401 +1 2020 10 07 04 00 0.02854711 +1 2020 10 07 05 00 0.028164540000000002 +1 2020 10 07 06 00 0.020105894 +1 2020 10 07 07 00 0.00052780524 +1 2020 10 07 08 00 0.0014239921999999999 +1 2020 10 07 09 00 0.0008253555 +1 2020 10 07 10 00 0.007109331 +1 2020 10 07 11 00 0.06388118999999999 +1 2020 10 07 12 00 0.039464365 +1 2020 10 07 13 00 0.025624737 +1 2020 10 07 14 00 0.06524496 +1 2020 10 07 15 00 0.08188297 +1 2020 10 07 16 00 0.0747276 +1 2020 10 07 17 00 0.04946063 +1 2020 10 07 18 00 0.08798275 +1 2020 10 07 19 00 0.019656029 +1 2020 10 07 20 00 0.005929752700000001 +1 2020 10 07 21 00 0.0019057406 +1 2020 10 07 22 00 0.0037654318000000003 +1 2020 10 07 23 00 0.0055011405 +1 2020 10 08 00 00 0.010166303 +1 2020 10 08 01 00 0.010956223999999999 +1 2020 10 08 02 00 0.0071553804 +1 2020 10 08 03 00 0.0034112018 +1 2020 10 08 04 00 0.00032943353 +1 2020 10 08 05 00 7.275957999999999e-09 +1 2020 10 08 06 00 7.275957999999999e-09 +1 2020 10 08 07 00 7.275957999999999e-09 +1 2020 10 08 08 00 7.275957999999999e-09 +1 2020 10 08 09 00 7.275957999999999e-09 +1 2020 10 08 10 00 7.275957999999999e-09 +1 2020 10 08 11 00 7.275957999999999e-09 +1 2020 10 08 12 00 7.275957999999999e-09 +1 2020 10 08 13 00 7.275957999999999e-09 +1 2020 10 08 14 00 0.0019942962999999997 +1 2020 10 08 15 00 0.017066632 +1 2020 10 08 16 00 0.1024281 +1 2020 10 08 17 00 0.09838992 +1 2020 10 08 18 00 0.04274804 +1 2020 10 08 19 00 0.013393306 +1 2020 10 08 20 00 0.005975802 +1 2020 10 08 21 00 0.0024016626 +1 2020 10 08 22 00 0.01663802 +1 2020 10 08 23 00 0.048688416 diff --git a/swmmanywhere/defs/swmm_conversion.yml b/swmmanywhere/defs/swmm_conversion.yml new file mode 100644 index 00000000..eddd4fbc --- /dev/null +++ b/swmmanywhere/defs/swmm_conversion.yml @@ -0,0 +1,54 @@ +# Adapted from SWMMIO: https://github.com/pyswmm/swmmio/tree/master +# 'columns' gives the names and order of input of data that SWMM expects +# 'iwcolumns' gives that names of columns in our dataframes that match with the +# corresponding values in 'columns' +# / in 'iwcolumns' sets a default value for that column +CONDUITS: + columns: [Name, InletNode, OutletNode, Length, Roughness, InOffset, OutOffset, InitFlow, MaxFlow] + iwcolumns: [id, u, v, length, roughness, /0, /0, /0, capacity] +INFILTRATION: + columns: [Subcatchment, Suction, HydCon, IMDmax, ] + iwcolumns: [subcatchment, /0, /0, /0, /0, /0] +JUNCTIONS: + columns: [Name, InvertElev, MaxDepth, InitDepth, SurchargeDepth, PondedArea] + iwcolumns: [id, chamber_floor_elevation, max_depth, /0, surcharge_depth, flooded_area] +OUTFALLS: + columns: [Name, InvertElev, OutfallType, StageOrTimeseries, TideGate, RouteTo] + iwcolumns: [id, chamber_floor_elevation, /FREE, / , /NO, /*] +STORAGE: + columns: [Name, InvertElev, MaxD, InitDepth, StorageCurve, Coefficient, Exponent, + Constant, EvapFrac, SuctionHead, Conductivity, InitialDeficit] + iwcolumns: [id, chamber_floor_elevation, max_depth, /0, /FUNCTIONAL, /0, /0, + manhole_area, /0, /0, /0, /0] +SUBCATCHMENTS: + columns: [Name, Raingage, Outlet, Area, PercImperv, Width, PercSlope, + CurbLength, SnowPack] + iwcolumns: [subcatchment, rain_gage, id, area, rc, width, slope, + /0, /empty] +SUBAREAS: + columns: [Name, N-Imperv, N-Perv, S-Imperv, S-Perv, PctZero, RouteTo, PctRouted] + iwcolumns: [subcatchment, /1, /0, /0, /0, /0, /OUTLET, /0] +XSECTIONS: + columns: [Link, Shape, Geom1, Geom2, Geom3, Geom4, Barrels, XX] + iwcolumns: [id, /CIRCULAR, diameter, /0, /0, /0, /1, /*] +COORDINATES: + columns: [Name, X, Y] + iwcolumns: [id, x, y] +VERTICES: + columns: [Name, X, Y] + columns: [id, x, y] +Polygons: + columns: [Name, X, Y] + iwcolumns: [subcatchment, x, y] +POLYGONS: + columns: [Name, X, Y] + iwcolumns: [subcatchment, x, y] +MAP: + columns: [Param, x1, y1, x2, y2] + iwcolumns: [/DIMENSIONS, x1, y1, x2, y2] +RAINGAGES: + columns: [Name,Format,Interval,SCF,Source,Filename,StationID,Unit] + iwcolumns: [name,/INTENSITY,interval,/1,/FILE,fid,/1,unit] +SYMBOLS: + columns: [Gage, X, Y] + iwcolumns: [name, x, y] \ No newline at end of file diff --git a/swmmanywhere/post_processing.py b/swmmanywhere/post_processing.py new file mode 100644 index 00000000..51fb633b --- /dev/null +++ b/swmmanywhere/post_processing.py @@ -0,0 +1,454 @@ +# -*- coding: utf-8 -*- +"""Created 2024-01-22. + +A module containing functions to format and write processed data into SWMM .inp +files. + +@author: Barnaby Dobson +""" +import os +import re +import shutil +from pathlib import Path +from typing import Literal + +import geopandas as gpd +import numpy as np +import pandas as pd +import yaml + + +def synthetic_write(model_dir: str): + """Load synthetic data and write to SWMM input file. + + Loads nodes, edges and subcatchments from synthetic data, assumes that + these are all located in `model_dir`. Fills in appropriate default values + for many SWMM parameters. More parameters are available to edit (see + defs/swmm_conversion.yml). Identifies outfalls and automatically ensures + that they have only one link to them (as is required by SWMM). Formats + (with format_to_swmm_dict) and writes (with data_dict_to_inp) the data to + a SWMM input (.inp) file. + + Args: + model_dir (str): model directory address. Assumes a format along the + lines of 'model_2', where the number is the model number. + """ + # TODO these node/edge names are probably not good or extendible defulats + # revisit once overall software architecture is more clear. + nodes = gpd.read_file(os.path.join(model_dir, + 'pipe_by_pipe_nodes.geojson')) + edges = gpd.read_file(os.path.join(model_dir, + 'pipe_by_pipe_edges.geojson')) + subs = gpd.read_file(os.path.join(model_dir, + 'subcatchments.geojson')) + + # Extract SWMM relevant data + edges = edges[['u','v','diameter','length']] + nodes = nodes[['id', + 'x', + 'y', + 'chamber_floor_elevation', + 'surface_elevation']] + subs = subs[['id', + 'geometry', + 'area', + 'slope', + 'width', + 'rc']] + + # Nodes + nodes['id'] = nodes['id'].astype(str) + nodes['max_depth'] = nodes.surface_elevation - nodes.chamber_floor_elevation + nodes['surcharge_depth'] = 0 + nodes['flooded_area'] = 100 # TODO arbitrary... not sure how to calc this + nodes['manhole_area'] = 0.5 + + # Subs + subs['id'] = subs['id'].astype(str) + subs['subcatchment'] = subs['id'] + '-sub' + subs['rain_gage'] = 1 # TODO revise when revising storms + + # Edges + edges['u'] = edges['u'].astype(str) + edges['v'] = edges['v'].astype(str) + edges['roughness'] = 0.01 + edges['capacity'] = 1E10 # capacity in swmm is a hard limit + + # Outfalls (create new nodes that link to the stores connected tothem + outfalls = nodes.loc[~nodes.id.isin(edges.u)].copy() + outfalls['id'] = outfalls['id'] + '_outfall' + + # Reduce elevation to ensure flow + outfalls['chamber_floor_elevation'] -= 1 + outfalls['x'] -= 1 + outfalls['y'] -= 1 + + # Link stores to outfalls + new_edges = edges.iloc[0:outfalls.shape[0]].copy() + new_edges['u'] = outfalls['id'].str.replace('_outfall','').values + new_edges['v'] = outfalls['id'].values + new_edges['diameter'] = 15 # TODO .. big pipe to enable all outfall... + new_edges['length'] = 1 + + # Append new edges + edges = pd.concat([edges, new_edges], ignore_index = True) + + # Name all edges + edges['id'] = edges.u.astype(str) + '-' + edges.v.astype(str) + + # Create event + # TODO will need some updating if multiple rain gages + # TODO automatically match units to storm.csv? + event = {'name' : '1', + 'unit' : 'mm', + 'interval' : '01:00', + 'fid' : 'storm.dat' # overwritten at runtime + } + + # Locate raingage(s) on the map + symbol = {'x' : nodes.x.min(), + 'y' : nodes.y.min(), + 'name' : '1' # matches event name(s) + } + + # Template SWMM input file + existing_input_file = os.path.join(os.path.dirname(__file__), + 'defs', + 'basic_drainage_all_bits.inp') + + # New input file + model_number = model_dir.split('_')[-1] + new_input_file = os.path.join(model_dir, + 'model_{0}.inp'.format(model_number)) + + # Format to dict + data_dict = format_to_swmm_dict(nodes, + outfalls, + edges, + subs, + event, + symbol) + + # Write new input file + data_dict_to_inp(data_dict, existing_input_file, new_input_file) + + +def overwrite_section(data: np.ndarray, + section: str, + fid: str): + """Overwrite a section of a SWMM .inp file with new data. + + Args: + data (np.ndarray): Data array to be written to the SWMM .inp file. + section (str): Section of the SWMM .inp file to be overwritten. + fid (str): File path to the SWMM .inp file. + + Example: + data = np.array([ + ['1', '1', '1', '1.166', '100', '500', '0.5', '0', 'empty'], + ['2', '1', '1', '1.1', '100', '500', '0.5', '0', 'empty'], + ['3', '1', '1', '2', '100', '400', '0.5', '0', 'empty']]) + fid = 'my_pre_existing_swmm_input_file.inp' + section = '[SUBCATCHMENTS]' + overwrite_section(data, section, fid) + """ + # Read the existing SWMM .inp file + with open(fid, 'r') as infile: + lines = infile.readlines() + + # Create a flag to indicate whether we are within the target section + within_target_section = False + + # Iterate through the lines and make modifications as needed + with open(fid, 'w') as outfile: + for ix, line in enumerate(lines): + if line.strip() != section and re.search(r'\[.*?\]', line): + within_target_section = False + + if line.strip() != section and not within_target_section: + outfile.write(line) # Write lines outside the target section + + if line.strip() != section: + continue + + within_target_section = True + outfile.write(line) # Write the start section header + + # Write headers + i = 1 + while lines[ix + i][0] == ';': + outfile.write(lines[ix + i]) # Write column headers + i += 1 + + example_line = lines[ix + i] + print('example_line {1}: {0}'.format( + example_line.replace('\n', ''), section)) + print('note - this line must have at least as many column') + print('entries as all other rows in this section\n') + pattern = r'(\s+)' + + # Find all matches of the pattern in the input line + matches = re.findall(pattern, example_line) + + # Calculate the space counts by taking the length of each match + space_counts = [len(x) + len(y) + for x, y in zip(matches, example_line.split())] + if len(space_counts) == 0: + if data.shape[0] != 0: + print('no template for data?') + continue + + space_counts[-1] -= 1 + new_text = '' + if data is None: + continue + + for i, row in enumerate(data): + if section == '[CONTROLS]': + new_text = row + continue + + formatted_row = [] + for x, y in zip(row, space_counts): + formatted_value = '{0:<{1}}'.format(x, max(y, len(str(x)) + 1)) + formatted_row.append(formatted_value) + new_line = '{0}\n'.format(''.join(formatted_row)) + new_text += new_line + + outfile.write(new_text) # Write the new content + outfile.write('\n') + +def change_flow_routing(routing_method: Literal["KINWAVE", "DYNWAVE", "STEADY"], + file_path: str | Path)-> None: + """Replace the flow routing method in a SWMM inp file with a new method, in-place. + + Args: + file_path : str or Path + Path to the SWMM inp file to be modified. + routing_method : {"KINWAVE", "DYNWAVE", "STEADY"} + The new flow routing method to be used. Available options are: + ``KINWAVE``, ``DYNWAVE``, or ``STEADY``. + """ + if routing_method.upper() not in ["KINWAVE", "DYNWAVE", "STEADY"]: + raise ValueError( + "routing_method must be one of 'KINWAVE', 'DYNWAVE', or 'STEADY'." + ) + file_path = Path(file_path) + updated_contents = re.sub( + r'^FLOW_ROUTING\s+.*$', + f'FLOW_ROUTING {routing_method.upper()}', + file_path.read_text(), + flags=re.MULTILINE, + ) + file_path.write_text(updated_contents) + +def data_dict_to_inp(data_dict: dict[str, np.ndarray], + base_input_file: str, + new_input_file: str, + routing: Literal["KINWAVE", "DYNWAVE", "STEADY"] = "DYNWAVE"): + """Write a SWMM .inp file from a dictionary of data arrays. + + Args: + data_dict (dict[str, np.ndarray]): Dictionary of data arrays. Where + each key is a SWMM section and each value is a numpy array of + data to be written to that section. The existing section is + overwritten + base_input_file (str): File path to the example/template .inp file. + new_input_file (str): File path to the new SWMM .inp file. + routing (str, optional): Flow routing method (KINWAVE, DYNWAVE, + STEADY). Defaults to "DYNWAVE". + """ + shutil.copy2(base_input_file, new_input_file) + + # Write the inp file + for key, data in data_dict.items(): + print(key) + start_section = '[{0}]'.format(key) + + overwrite_section(data, start_section, new_input_file) + + # Set the flow routing + change_flow_routing(routing, new_input_file) + +def explode_polygon(row): + """Explode a polygon into a DataFrame of coordinates. + + Args: + row (pd.Series): A row of a GeoDataFrame containing a polygon. + + Example: + >>> import geopandas as gpd + >>> from shapely.geometry import Polygon + >>> df = pd.Series({'subcatchment' : '1', + ... 'geometry' : Polygon([(0,0), (1,0), + ... (1,1), (0,1)])}) + >>> explode_polygon(df) + x y subcatchment + 0 0 0 1 + 1 1 0 1 + 2 1 1 1 + 3 0 1 1 + 4 0 0 1 + """ + # Get the vertices of the polygon + vertices = list(row['geometry'].exterior.coords) + + # Create a new DataFrame for this row + df = pd.DataFrame(columns = ['x','y'], + data =vertices) + df['subcatchment'] = row['subcatchment'] + return df + +def format_to_swmm_dict(nodes, + outfalls, + conduits, + subs, + event, + symbol): + """Format data to a dictionary of data arrays with columns matching SWMM. + + These data are the parameters of all assets that are written to the SWMM + input file. More parameters are available to edit (see + defs/swmm_conversion.yml). + + Args: + nodes (pd.DataFrame): GeoDataFrame of nodes. With at least columns: + 'id', 'x', 'y', 'max_depth', 'chamber_floor_elevation', + 'manhole_area'. + outfalls (pd.DataFrame): GeoDataFrame of outfalls. With at least + columns: 'id', 'chamber_floor_elevation'. + conduits (pd.DataFrame): GeoDataFrame of conduits. With at least + columns: 'id', 'u', 'v', 'length', 'roughness', 'shape_swmm', + 'diameter'. + subs (gpd.GeoDataFrame): GeoDataFrame of subcatchments. With at least + columns: 'subcatchment', 'rain_gage', 'id', 'area', 'rc', 'width', + 'geometry', + event (dict): Dict describing storm event. With at least + keys: 'name', 'unit', 'interval', 'fid'. + symbol (dict): Dict with coordinates of rain gage. With at least keys: + 'x', 'y', 'name'. + + Example: + >>> import geopandas as gpd + >>> from shapely.geometry import Point + >>> nodes = gpd.GeoDataFrame({'id' : ['node1', 'node2'], + ... 'x' : [0, 1], + ... 'y' : [0, 1], + ... 'max_depth' : [1, 1], + ... 'chamber_floor_elevation' : [1, 1], + ... 'manhole_area' : [1,1] + ... }) + >>> outfalls = gpd.GeoDataFrame({'id' : ['outfall3'], + ... 'chamber_floor_elevation' : [1], + ... 'x' : [0], + ... 'y' : [0]}) + >>> conduits = gpd.GeoDataFrame({'id' : ['link1','link2'], + ... 'u' : ['node1','node2'], + ... 'v' : ['node2','outfall3'], + ... 'length' : [1,1], + ... 'roughness' : [1,1], + ... 'shape_swmm' : ['CIRCULAR','CIRCULAR'], + ... 'diameter' : [1,1], + ... 'capacity' : [0.1,0.1] + ... }) + >>> subs = gpd.GeoDataFrame({'subcatchment' : ['sub1'], + ... 'rain_gage' : ['1'], + ... 'id' : ['node1'], + ... 'area' : [1], + ... 'rc' : [1], + ... 'width' : [1], + ... 'slope' : [0.001], + ... 'geometry' : [sgeom.Polygon([(0,0), (1,0), + ... (1,1), (0,1)])]}) + >>> rain_fid = os.path.join(os.path.dirname(os.path.abspath(__file__)), + ... '..', + ... 'swmmanywhere', + ... 'defs', + ... 'storm.dat') + >>> event = {'name' : '1', + ... 'unit' : 'mm', + ... 'interval' : 1, + ... 'fid' : rain_fid} + >>> symbol = {'x' : 0, + ... 'y' : 0, + ... 'name' : 'name'} + >>> data_dict = stt.format_to_swmm_dict(nodes, + ... outfalls, + ... conduits, + ... subs, + ... event, + ... symbol) + """ + # Get the directory of the current module + current_dir = os.path.dirname(os.path.abspath(__file__)) + + # TODO use 'load_yaml_from_defs' + # Create the path to iso_converter.yml + iso_path = os.path.join(current_dir, + "defs", + "swmm_conversion.yml") + + + # Load conversion mapping from YAML file + with open(iso_path, "r") as file: + conversion_dict = yaml.safe_load(file) + + ## Create nodes, coordinates and map dimensions + dims = {'x1' : nodes.x.min(), + 'y1' : nodes.y.min(), + 'x2' : nodes.x.max(), + 'y2' : nodes.y.max(), + } + dims = {x : str(y) + ' ' for x,y in dims.items()} + + map_dimensions = pd.Series(dims).reset_index().set_index('index').T + polygons = subs[['subcatchment','geometry']].copy() + + # Format dicts to DataFrames + event = pd.Series(event).reset_index().set_index('index').T + symbol = pd.Series(symbol).reset_index().set_index('index').T + + # Apply the function to each row + polygons = pd.concat(polygons.apply(explode_polygon, axis=1).tolist()) + + ## Specify sections + shps = {'SUBCATCHMENTS' : subs, + 'CONDUITS' : conduits, + 'OUTFALLS' : outfalls, + 'STORAGE' : nodes, + 'XSECTIONS' : conduits, + 'SUBAREAS' : subs, + 'INFILTRATION' : subs, + 'COORDINATES' : pd.concat([nodes, outfalls],axis=0), + 'MAP' : map_dimensions, + 'Polygons' : polygons, + 'PUMPS' : None, + 'ORIFICES' : None, + 'WEIRS' : None, + 'OUTLETS' : None, + 'JUNCTIONS' : None, + 'RAINGAGES' : event, + 'SYMBOLS' : symbol + } + + # Fill backslash columns and store data in data_dict in the correct order + import numpy.typing as npt + + def _fill_backslash_columns(shp: pd.DataFrame | None, + key: str) -> npt.ArrayLike | None: + if shp is None: + return None + + # Extract SWMM order and default values + columns = conversion_dict[key]['iwcolumns'] + shp = shp.fillna(0) + + # Find columns with a default specified + cols_default = [c[1:] for c in columns if c.startswith("/")] + + # Fill columns with defaults + shp[['/' + c for c in cols_default]] = np.array(cols_default).T + return shp[columns].values + + data_dict = {key: _fill_backslash_columns(shp, key) for key, shp in shps.items()} + return data_dict \ No newline at end of file diff --git a/tests/test_post_processing.py b/tests/test_post_processing.py new file mode 100644 index 00000000..16a811e7 --- /dev/null +++ b/tests/test_post_processing.py @@ -0,0 +1,237 @@ +import difflib +import filecmp +import os +import shutil +import tempfile + +import geopandas as gpd +import pandas as pd +import pyswmm +from shapely import geometry as sgeom + +from swmmanywhere import post_processing as stt + + +def test_overwrite_section(): + """Test the overwrite_section function. + + All this tests is that the data is written to the file. + """ + # Copy the example file to a temporary file + fid = os.path.join(os.path.dirname(os.path.abspath(__file__)), + '..', + 'swmmanywhere', + 'defs', + 'basic_drainage_all_bits.inp') + temp_fid = 'temp.inp' + shutil.copy(fid, temp_fid) + try: + data = [["1","1","1","1.166","100","500","0.5","0","empty"], + ["2","1","1","1.1","100","500","0.5","0","empty"], + ["subca_3","1","1","2","100","400","0.5","0","empty"]] + + section = '[SUBCATCHMENTS]' + stt.overwrite_section(data, section, temp_fid) + with open(temp_fid, 'r') as file: + content = file.read() + assert 'subca_3' in content + finally: + os.remove(temp_fid) + +def test_change_flow_routing(): + """Test the change_flow_routing function.""" + # Copy the example file to a temporary file + fid = os.path.join(os.path.dirname(os.path.abspath(__file__)), + '..', + 'swmmanywhere', + 'defs', + 'basic_drainage_all_bits.inp') + temp_fid = 'temp.inp' + shutil.copy(fid, temp_fid) + try: + new_routing = 'STEADY' + stt.change_flow_routing(new_routing, temp_fid) + with open(temp_fid, 'r') as file: + content = file.read() + assert 'STEADY' in content + finally: + os.remove(temp_fid) + +def test_data_input_dict_to_inp(): + """Test the data_input_dict_to_inp function. + + All this tests is that the data is written to a new file. + """ + data_dict = {'SUBCATCHMENTS': + [["1","1","1","1.166","100","500","0.5","0","empty"], + ["2","1","1","1.1","100","500","0.5","0","empty"], + ["subca_3","1","1","2","100","400","0.5","0","empty"]] + } + + fid = os.path.join(os.path.dirname(os.path.abspath(__file__)), + '..', + 'swmmanywhere', + 'defs', + 'basic_drainage_all_bits.inp') + temp_fid = 'temp.inp' + shutil.copy(fid, temp_fid) + try: + stt.data_dict_to_inp(data_dict, + fid, + temp_fid) + with open(temp_fid, 'r') as file: + content = file.read() + assert 'subca_3' in content + finally: + os.remove(temp_fid) + +def test_explode_polygon(): + """Test the explode_polygon function.""" + df = pd.Series({'subcatchment' : '1', + 'geometry' : sgeom.Polygon([(0,0), (1,0), + (1,1), (0,1)])}) + result = stt.explode_polygon(df) + assert result.shape[0] == 5 + assert result.loc[3,'y'] == 1 + +def generate_data_dict(): + """Generate a data dict for testing.""" + nodes = gpd.GeoDataFrame({'id' : ['node1', 'node2'], + 'x' : [0, 1], + 'y' : [0, 1], + 'max_depth' : [1, 1], + 'chamber_floor_elevation' : [1, 1], + 'surface_elevation' : [2, 2], + 'manhole_area' : [0.5,0.5] + }) + outfalls = gpd.GeoDataFrame({'id' : ['node2_outfall'], + 'chamber_floor_elevation' : [0], + 'surface_elevation' : [0], + 'x' : [0], + 'y' : [0]}) + conduits = gpd.GeoDataFrame({'id' : ['node1-node2','node2-node2_outfall'], + 'u' : ['node1','node2'], + 'v' : ['node2','node2_outfall'], + 'length' : [1,1], + 'roughness' : [0.01,0.01], + 'shape_swmm' : ['CIRCULAR','CIRCULAR'], + 'diameter' : [1,15], + 'capacity' : [1E10,1E10] + }) + subs = gpd.GeoDataFrame({'subcatchment' : ['node1-sub'], + 'rain_gage' : ['1'], + 'id' : ['node1'], + 'area' : [1], + 'rc' : [1], + 'width' : [1], + 'slope' : [0.001], + 'geometry' : [sgeom.Polygon([(0,0), (1,0), + (1,1), (0,1)])]}) + rain_fid = 'storm.dat' + event = {'name' : '1', + 'unit' : 'mm', + 'interval' : '01:00', + 'fid' : rain_fid} + symbol = {'x' : 0, + 'y' : 0, + 'name' : '1'} + + return {'nodes' : nodes, + 'outfalls' : outfalls, + 'conduits' : conduits, + 'subs' : subs, + 'event' : event, + 'symbol' : symbol} + +def test_synthetic_write(): + """Test the synthetic_write function.""" + data_dict = generate_data_dict() + model_dir = tempfile.mkdtemp('_1') + try: + # Write the model with synthetic_write + nodes = gpd.GeoDataFrame(data_dict['nodes']) + nodes.geometry = gpd.points_from_xy(nodes.x, nodes.y) + nodes.to_file( + os.path.join( + model_dir, + 'pipe_by_pipe_nodes.geojson')) + nodes = nodes.set_index('id') + edges = gpd.GeoDataFrame(pd.DataFrame(data_dict['conduits']).iloc[[0]]) + edges.geometry = [sgeom.LineString([nodes.loc[u,'geometry'], + nodes.loc[v,'geometry']]) + for u,v in zip(edges.u, edges.v)] + edges.to_file( + os.path.join( + model_dir, + 'pipe_by_pipe_edges.geojson')) + subs = data_dict['subs'].copy() + subs['subcatchment'] = ['node1'] + subs.to_file(os.path.join(model_dir, + 'subcatchments.geojson')) + stt.synthetic_write(model_dir) + + # Write the model with data_dict_to_inp + comparison_file = os.path.join(model_dir, "model_base.inp") + template_fid = os.path.join(os.path.dirname(os.path.abspath(__file__)), + '..', + 'swmmanywhere', + 'defs', + 'basic_drainage_all_bits.inp') + stt.data_dict_to_inp(stt.format_to_swmm_dict(**data_dict), + template_fid, + comparison_file) + + # Compare + new_input_file = os.path.join(model_dir, "model_1.inp") + are_files_identical = filecmp.cmp(new_input_file, + comparison_file, + shallow=False) + if not are_files_identical: + with open(new_input_file, + 'r') as file1, open(comparison_file, + 'r') as file2: + diff = difflib.unified_diff( + file1.readlines(), + file2.readlines(), + fromfile=new_input_file, + tofile=comparison_file, + ) + print(''.join(diff)) + assert are_files_identical, "The files are not identical" + finally: + pass + # shutil.rmtree(model_dir) + +def test_format_to_swmm_dict(): + """Test the format_format_to_swmm_dict function. + + Writes a formatted dict to a model and checks that it runs without + error. + """ + data_dict = generate_data_dict() + data_dict = stt.format_to_swmm_dict(**data_dict) + + + fid = os.path.join(os.path.dirname(os.path.abspath(__file__)), + '..', + 'swmmanywhere', + 'defs', + 'basic_drainage_all_bits.inp') + rain_fid = os.path.join(os.path.dirname(os.path.abspath(__file__)), + '..', + 'swmmanywhere', + 'defs', + 'storm.dat') + temp_dir = tempfile.mkdtemp() + shutil.copy(rain_fid, os.path.join(temp_dir,'storm.dat')) + temp_fid = os.path.join(temp_dir,'temp.inp') + try: + stt.data_dict_to_inp(data_dict, + fid, + temp_fid) + with pyswmm.Simulation(temp_fid) as sim: + for ind, step in enumerate(sim): + pass + finally: + shutil.rmtree(temp_dir) + \ No newline at end of file