Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def population_impacted_metrics(wn):
# Compute population per node
pop = wntr.metrics.population(wn)
total_population = pop.sum()
print("Total population: " + str(total_population))
wntr.graphics.plot_network(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()
pressure = results.node['pressure'].loc[:,wn.junction_name_list]
pop_impacted = wntr.metrics.population_impacted(pop, pressure, np.less, 40)
plt.figure()
pop_impacted.sum(axis=1).plot(title='Total population with pressure < 40 m')
nodes_impacted = wntr.metrics.query(pressure, np.less, 40)
wntr.graphics.plot_network(wn, node_attribute=nodes_impacted.any(axis=0), node_size=40,
title='Nodes impacted')
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)
wntr.network.draw_graph(wn, node_attribute=nodes_impacted.any(axis=0), node_size=40,
title='Nodes impacted')
# Copute network cost
network_cost = wntr.metrics.cost(wn)
print "Network cost: $" + str(round(network_cost,2))
# COmpute green house gas emissions
network_ghg = wntr.metrics.ghg_emissions(wn)
print "Network GHG emissions: " + str(round(network_ghg,2))
# Isolate node results at consumer nodes (nzd = non-zero demand)
nzd_junctions = wn.query_node_attribute('base_demand', np.greater, 0,
node_type=wntr.network.Junction).keys()
node_results = results.node.loc[:, :, nzd_junctions]
# Compute population per node
pop = wntr.metrics.population(wn)
# Compute FDV, FDD, and FDQ, change average_times and average_nodes to average results
quality_upper_bound = 0.0035 # kg/m3 (3.5 mg/L)
demand_factor = 0.9
average_times = False
average_nodes = False
fdv = wntr.metrics.fdv(node_results, average_times, average_nodes)
fdd = wntr.metrics.fdd(node_results, demand_factor, average_times, average_nodes)
#fdq = wntr.metrics.fdq(node_results, quality_upper_bound, average_times, average_nodes)
# Plot results
if average_times == False and average_nodes == False:
fdv.plot(ylim=(-0.05, 1.05), legend=False, title='FDV for each node and time')
fdd.plot(ylim=(-0.05, 1.05),legend=False, title='FDD for each node and time')
# Fraction of nodes with reduced demand
fdd_fraction_impacted = 1-fdd.sum(axis=1)/fdd.shape[1]
plt.figure()
fdd_fraction_impacted.plot(ylim=(-0.05, 1.05), title='Fraction of nodes not receiving adaquate demand')
# Population impacted by reduced demand
fdd_pop_impacted = wntr.metrics.population_impacted(pop, fdd, np.less, 1)
plt.figure()
fdd_pop_impacted.plot(legend=False, title='Population not receiving adaquate demand')
# Timeseries of fraction of population not meeting demand
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()
EC.sum(axis=1).plot(title='Extent of contamination')
str(time_of_failure) + ', End Time: ' + \
str(time_of_failure+duration_of_failure))
results[i] = sim.run_sim()
# Reload the water network model
f=open('wn.pickle','rb')
wn = pickle.load(f)
f.close()
# Plot water service availability and tank water level for each realization
for i in results.keys():
# Water service availability at each junction and time
expected_demand = wntr.metrics.expected_demand(wn)
demand = results[i].node['demand'].loc[:,wn.junction_name_list]
wsa_nt = wntr.metrics.water_service_availability(expected_demand, demand)
# Average water service availability at each time
wsa_t = wntr.metrics.water_service_availability(expected_demand.sum(axis=1),
demand.sum(axis=1))
# Tank water level
tank_level = results[i].node['pressure'].loc[:,wn.tank_name_list]
# Plot results
plt.figure()
plt.subplot(2,1,1)
wsa_nt.plot(ax=plt.gca(), legend=False)
wsa_t.plot(ax=plt.gca(), label='Average', color='k', linewidth=3.0, legend=False)
plt.ylim( (-0.05, 1.05) )
plt.ylabel('Water service availability')
sim = wntr.sim.WNTRSimulator(wn, mode='PDD')
print('Pipe Breaks: ' + str(pipes_to_fail) + ', Start Time: ' + \
str(time_of_failure) + ', End Time: ' + \
str(time_of_failure+duration_of_failure))
results[i] = sim.run_sim()
# Reload the water network model
f=open('wn.pickle','rb')
wn = pickle.load(f)
f.close()
# Plot water service availability and tank water level for each realization
for i in results.keys():
# Water service availability at each junction and time
expected_demand = wntr.metrics.expected_demand(wn)
demand = results[i].node['demand'].loc[:,wn.junction_name_list]
wsa_nt = wntr.metrics.water_service_availability(expected_demand, demand)
# Average water service availability at each time
wsa_t = wntr.metrics.water_service_availability(expected_demand.sum(axis=1),
demand.sum(axis=1))
# Tank water level
tank_level = results[i].node['pressure'].loc[:,wn.tank_name_list]
# Plot results
plt.figure()
plt.subplot(2,1,1)
wsa_nt.plot(ax=plt.gca(), legend=False)
wsa_t.plot(ax=plt.gca(), label='Average', color='k', linewidth=3.0, legend=False)
for t in np.arange(0, 24*3600+1,3600):
attr = results.link['flowrate'].loc[t,:]
G_flowrate_t.weight_graph(link_attribute=attr)
entropy = wntr.metrics.entropy(G_flowrate_t)
shat.append(entropy[1])
plt.figure()
plt.plot(shat)
plt.ylabel('System Entropy')
plt.xlabel('Time, hr')
print("Entropy")
print(" Mean: " + str(np.mean(shat)))
print(" Max: " + str(np.nanmax(shat)))
print(" Min: " + str(np.nanmin(shat)))
# Compute water service availability, for each junction
expected_demand = wntr.metrics.expected_demand(wn)
demand = results.node['demand'].loc[:,wn.junction_name_list]
wsa = wntr.metrics.water_service_availability(expected_demand.sum(axis=0),
demand.sum(axis=0))
wntr.graphics.plot_network(wn, node_attribute=wsa, node_size=40, node_range=[0,1],
title='Water service availability, averaged over times')
# Simulate hydraulics and water quality
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)
wn.options.duration = 48*3600
for name, node in wn.junctions():
node.nominal_pressure = 60
sim = wntr.sim.ScipySimulator(wn, pressure_driven=True)
results = sim.run_sim()
# Isolate node results at consumer nodes (nzd = non-zero demand)
nzd_junctions = wn.query_node_attribute('base_demand', np.greater, 0,
node_type=wntr.network.Junction).keys()
node_results = results.node.loc[:, :, nzd_junctions]
# Compute population per node
pop = wntr.metrics.population(wn)
# Compute FDV, FDD, and FDQ, change average_times and average_nodes to average results
quality_upper_bound = 0.0035 # kg/m3 (3.5 mg/L)
demand_factor = 0.9
average_times = False
average_nodes = False
fdv = wntr.metrics.fdv(node_results, average_times, average_nodes)
fdd = wntr.metrics.fdd(node_results, demand_factor, average_times, average_nodes)
#fdq = wntr.metrics.fdq(node_results, quality_upper_bound, average_times, average_nodes)
# Plot results
if average_times == False and average_nodes == False:
fdv.plot(ylim=(-0.05, 1.05), legend=False, title='FDV for each node and time')
fdd.plot(ylim=(-0.05, 1.05),legend=False, title='FDD for each node and time')
# Fraction of nodes with reduced demand
def cost_ghg_metrics(wn):
# Compute network cost
network_cost = wntr.metrics.annual_network_cost(wn)
print("Network cost: $" + str(round(network_cost,2)))
# Compute green house gas emissions
network_ghg = wntr.metrics.annual_ghg_emissions(wn)
print("Network GHG emissions: " + str(round(network_ghg,2)))
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)
print "Average node pressure: " +str(pressure.mean().mean()) + " m"
wntr.network.draw_graph(wn, node_attribute=pressure.min(axis=0), node_size=40,
title= 'Min pressure')
# Compute todini index
todini = wntr.metrics.todini(results.node,results.link,wn, P_lower)
plt.figure()
plt.plot(todini)
plt.ylabel('Todini Index')
plt.xlabel('Time, hr')
print "Todini Index"
print " Mean: " + str(np.mean(todini))
print " Max: " + str(np.max(todini))
print " Min: " + str(np.min(todini))
# Create a weighted graph for flowrate at time 36 hours
t = 36*3600
attr = results.link.loc['flowrate', t, :]
G_flowrate_36hrs = wn.get_graph_deep_copy()
G_flowrate_36hrs.weight_graph(link_attribute=attr)
# Compute betweenness-centrality time 36 hours