How to use the wntr.network.WaterNetworkModel function in wntr

from __future__ import print_function
import wntr

# Create a water network model
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)

# Simulate hydraulics
sim = wntr.sim.EpanetSimulator(wn)
results = sim.run_sim()

# Access node and link panels
print(results.node)
print(results.link)

# Access the pressure and demand at node '123' at 1 hour
print(results.node.loc[['pressure', 'demand'], 3600, '123'])

# Access the pressure for all nodes and times
print(results.node.loc['pressure', :, :])

# Plot time-series

import wntr
from scipy.stats import lognorm
import numpy as np
import matplotlib.pylab as plt

np.random.seed(12343)

# Create a water network model
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)

# Define the earthquake
wn = wntr.morph.node.scale_node_coordinates(wn, 1000)
epicenter = (32000,15000) # x,y location
magnitude = 7 # Richter scale
depth = 10000 # m, shallow depth
earthquake = wntr.scenario.Earthquake(epicenter, magnitude, depth)

# Compute PGA
R = earthquake.distance_to_epicenter(wn, element_type=wntr.network.Pipe)
pga = earthquake.pga_attenuation_model(R)  
wntr.graphics.plot_network(wn, link_attribute=pga,title='Peak Ground Acceleration (g)',
                           node_size=0, link_width=2)
  
# Define fragility curve  
FC = wntr.scenario.FragilityCurve()

import wntr
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

plt.close('all')

# Create a water network model
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)

# Simulate hydraulics
sim = wntr.sim.EpanetSimulator(wn)
results = sim.run_sim()

# Create list of node names
junctions = [name for name, node in wn.junctions()]

# Define pressure lower bound
P_lower = 21.09 # m (30 psi)

# Query pressure
pressure = results.node.loc['pressure', :, junctions]
mask = wntr.metrics.query(pressure, np.greater, P_lower)
pressure_regulation = mask.all(axis=0).sum() # True over all time
print "Fraction of nodes > 30 psi: " + str(pressure_regulation)

from __future__ import print_function
import wntr
import networkx as nx

inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)

# Get the WaterNetworkModel graph
G = wn.get_graph()

# Extract node and link properties
node_name = '123'
print(G.node[node_name])
print(G.adj[node_name])

# Compute betweenness centrality
node_degree = G.degree()
bet_cen = nx.betweenness_centrality(G)
wntr.graphics.plot_network(wn, node_attribute=bet_cen, node_size=30, 
                        title='Betweenness Centrality')

# Convert the digraph to an undirected graph

import wntr
import matplotlib.pyplot as plt

# Create a water network model
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)

# Define WQ scenarios
WQscenario = wntr.scenario.Waterquality('CHEM', '121', 'SETPOINT', 1000, 2*3600, 15*3600)

# Simulate hydraulics and water quality for each scenario
sim = wntr.sim.EpanetSimulator(wn)
results_CHEM = sim.run_sim(WQscenario)

MC = wntr.metrics.mass_contaminant_consumed(results_CHEM.node)
VC = wntr.metrics.volume_contaminant_consumed(results_CHEM.node, 0.001)
EC = wntr.metrics.extent_contaminant(results_CHEM.node, results_CHEM.link, wn, 0.001)

wntr.network.draw_graph(wn, node_attribute=MC.sum(axis=0), node_range = [0,400], node_size=40,
                      title='Total mass consumed')
    
plt.figure()

import wntr

# Generate a water network model
wn = wntr.network.WaterNetworkModel()
wn.add_pattern('pat1', [1])
wn.add_pattern('pat2', [1,2,3,4,5,6,7,8,9,10])
wn.add_junction('node1', base_demand=0.01, demand_pattern_name='pat1',
	elevation=100.0, coordinates=(1,2))
wn.add_junction('node2', base_demand=0.02, demand_pattern_name='pat2',
	elevation=50.0, coordinates=(1,3))
wn.add_pipe('pipe1', 'node1', 'node2', length=304.8, diameter=0.3048, roughness=100,
	minor_loss=0.0, status='OPEN')
wn.add_reservoir('res', base_head=125, head_pattern_name='pat1', coordinates=(0,2))
wn.add_pipe('pipe2', 'node1', 'res', length=100, diameter=0.3048, roughness=100,
	minor_loss=0.0, status='OPEN')
wn.options.duration = 24*3600
wn.options.hydraulic_timestep = 15*60
wn.options.pattern_timestep = 60*60

# Simulate hydraulics

print('running last 14 hours')
last_14_hours_of_results = sim.run_sim()
node_results = {}
link_results = {}
for key in first_10_hours_of_results.node.keys():
    node_results[key] = pd.concat([first_10_hours_of_results.node[key],
                                   last_14_hours_of_results.node[key]],axis=0)
for key in first_10_hours_of_results.link.keys():
    link_results[key] = pd.concat([first_10_hours_of_results.link[key],
                                   last_14_hours_of_results.link[key]],axis=0)
# node_results now has the exact same results as res1.node
# link_results now has the exact same results as res1.link

# Stop a simulation and save the water network model to a file
# Open the file and continue running
wn = wntr.network.WaterNetworkModel(inp_file)
wn.options.time.duration = 10*3600
wn.options.time.hydraulic_timestep = 3600
sim = wntr.sim.WNTRSimulator(wn)
first_10_hours_of_results = sim.run_sim()
f=open('pickle_example.pickle','wb')
pickle.dump(wn,f)
f.close()
f=open('pickle_example.pickle','rb')
new_wn = pickle.load(f)
f.close()
new_wn.options.time.duration = 24*3600
sim = wntr.sim.WNTRSimulator(new_wn)
print('running last 14 hours')
last_14_hours_of_results = sim.run_sim()
node_results = {}
link_results = {}

# Simulate using WNTR
wntr_sim = wntr.sim.WNTRSimulator(wn)
wntr_sim_results = wntr_sim.run_sim()

# Compare link flowrate results
plt.figure(figsize=(6,10))
plt.subplot(2,1,1)
plt.plot(epanet_sim_results.link['flowrate'])
plt.title('EPANET, Link Flowrate')
plt.subplot(2,1,2)
plt.plot(epanet_sim_results.link['flowrate'] - wntr_sim_results.link['flowrate'])
plt.title('EPANET - WNTR, Link Flowrate')

# Reset the water network and run again
wn = wntr.network.WaterNetworkModel(inp_file)
sim = wntr.sim.WNTRSimulator(wn)
res1 = sim.run_sim()
wn.reset_initial_values()
res2 = sim.run_sim()
# res2 has the exact same results as res1

# Stop and restart a simulation
wn = wntr.network.WaterNetworkModel(inp_file)
wn.options.time.duration = 10*3600
wn.options.time.hydraulic_timestep = 3600
sim = wntr.sim.WNTRSimulator(wn)
first_10_hours_of_results = sim.run_sim()
wn.options.time.duration = 24*3600
print('running last 14 hours')
last_14_hours_of_results = sim.run_sim()
node_results = {}

# Modify Network Stucture/Operations/Controls and simulate hydraulics
import wntr

# Create a water network model
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)

# Modify the water network model
wn.options.duration = 24*3600
wn.options.hydraulic_timestep = 1800
wn.options.report_timestep = 1800
for junction_name, junction in wn.nodes(wntr.network.Junction):
    junction.minimum_pressure = 0.0
    junction.nominal_pressure = 15.0

# Define pipe leaks
wn.split_pipe_with_junction('123','123A','123B','leak1')
leak1 = wn.get_node('leak1')
leak1.add_leak(area=3.14159/4.0*0.05**2, start_time=5*3600, end_time=20*3600)
wn.split_pipe_with_junction('225','225A','225B','leak2')
leak2 = wn.get_node('leak2')
leak2.add_leak(area=3.14159/4.0*0.1**2, start_time=7*3600, end_time=15*3600)

import wntr
import numpy as np
import matplotlib.pyplot as plt

# Create a water network model
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)
                   
# Compute population per node
pop = wntr.metrics.population(wn)
total_population = pop.sum()
print "Total population: " + str(total_population)
wntr.network.draw_graph(wn, node_attribute=pop, node_range = [0,400], node_size=40,
                      title='Population, Total = ' + str(total_population))
                      
# Find population and nodes impacted by pressure less than 40 m
sim = wntr.sim.EpanetSimulator(wn)
results = sim.run_sim()
junctions = [name for name, node in wn.junctions()]
pop_impacted = wntr.metrics.population_impacted(pop, results.node['pressure',:,junctions], np.less, 40)
plt.figure()
pop_impacted.sum(axis=1).plot(title='Total population with pressure < 40 m') 
nodes_impacted = wntr.metrics.query(results.node['pressure',:,junctions], np.less, 40)