import numpy as np
import pandas as pd
from scipy.spatial import distance
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import networkx as nx
import sd
from sd.dist import Distance
from scipy.stats import expon, lognorm
from sd.plot import plot_fragility_curve
np.random.seed(100000)
# read node and link positions
node = pd.read_csv('node_deg_fairfield.csv')
# link = pd.read_csv('edge classify.csv')
link = pd.read_csv('pipe_fairfield.csv')
# Define expicenter coordinates (check with above command)
ex=-13589706.28220519; #this is longitude of the epicenter
ey=4592436.47339457; #this is latitude of the epicenter
M=7 # Richter Magnitude
# trasforming latlong into mercetor coordinates
def tran(data):
utmxy = {}
rx=[]
ry=[]
for index, row in data.iterrows():
x = row['x']
y = row['y']
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:3395')
cz = transform(inProj,outProj,x,y)
# utmxy[int(row['id'])]=(c,z)
r1=cz[0]
r2=cz[1]
rx.append(r1)
ry.append(r2)
rx=np.array(rx)
ry=np.array(ry)
return rx,ry
node_id = list(node.id)
rx = list(rx)
ry = list(ry)
typ = list(node.type)
dn1 = {'id':node_id,'x':rx,'y':ry,'type':typ}
df1=pd.DataFrame(dn1)
# df1.set_index('id');
node=df1
node.head(2)
| id | x | y | type | |
|---|---|---|---|---|
| 0 | 106 | -1.358971e+07 | 4.592436e+06 | 1 |
| 1 | 107 | -1.358947e+07 | 4.592371e+06 | 1 |
# check the distance between two different points
L1=df1.x[3],df1.y[3]
L2=df1.x[5],df1.y[5]
dist = distance.euclidean(L1, L2)
print(dist/1000, "km")
0.11039776364390864 km
link['link_m']=list(Length(node,link))
# node types
a = node[node['type']==1]
b = node[node['type']==2]
c = node[node['type']==3]
# pipe types
p1 = link[link['type']==1]
p2 = link[link['type']==2]
# calculate distances,intensity from epicenter
r1, pga1, pgv1, pos1 = Distance.com_pga_dist(a,ex,ey,M)
r2, pga2, pgv2, pos2 = Distance.com_pga_dist(b,ex,ey,M)
r3, pga3, pgv3, pos3 = Distance.com_pga_dist(c,ex,ey,M)
# plot graps with various components
G1=nx.Graph()
G1.add_nodes_from(pos1)
G2=nx.Graph()
G2.add_nodes_from(pos2)
G3=nx.Graph()
G3.add_nodes_from(pos3)
nx.draw(G1,pos1,node_size=20,node_color='r',with_labels=False)
nx.draw(G2,pos2,node_size=20,node_color='b',with_labels=False)
nx.draw(G3,pos3,node_size=20,node_color='g',with_labels=False)
r, pga, pgv, pos = Distance.com_pga_dist(node,ex,ey,M)
d1 = []
for index, row in p1.iterrows():
stt = str(row['start_node'])+str(row['end_node']).rjust(10)
d1.append(stt)
d2 = []
for index, row in p2.iterrows():
stt = str(row['start_node'])+str(row['end_node']).rjust(10)
d2.append(stt)
# draw Graph of water network
G4 = nx.parse_edgelist(d1,nodetype=int)
G5 = nx.parse_edgelist(d2,nodetype=int)
nx.draw(G4,pos,node_size=10, node_color='k', edge_color='r',width=1.5)
nx.draw(G5,pos,node_size=10, node_color='k', edge_color='b',width=1.5)
nx.draw(G1,pos1,node_size=10,node_color='k',with_labels=False)
nx.draw(G2,pos2,node_size=10,node_color='k',with_labels=False)
nx.draw(G3,pos3,node_size=60,node_color='g',with_labels=False)
plt.savefig('Fairfield_Network.png', dpi = 600,bbox_inches='tight')
/Users/Ram/anaconda3/lib/python3.7/site-packages/networkx/drawing/nx_pylab.py:563: MatplotlibDeprecationWarning: The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead. if not cb.iterable(width): /Users/Ram/anaconda3/lib/python3.7/site-packages/networkx/drawing/nx_pylab.py:611: MatplotlibDeprecationWarning: The is_numlike function was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use isinstance(..., numbers.Number) instead. if cb.is_numlike(alpha):
d = []
for index, row in link.iterrows():
stt = str(row['start_node'])+str(row['end_node']).rjust(10)
d.append(stt)
r, pga, pgv, pos = Distance.com_pga_dist(node,ex,ey,M)
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
node_list = list(G.nodes)
pf = {'node_list':node_list}
n_list=pd.DataFrame(pf)
# n_list.tail(1)
node_pga = list(pga)
node_pgv = list(pgv)
ind= list(pos)
pf = {'node_pga':node_pga,'node_pgv':node_pgv,'pid':ind}
n_int=pd.DataFrame(pf)
# n_int.set_index('pid',inplace=True)
# n_int.tail(1)
n_out=n_list.merge(n_int, left_on='node_list', right_on='pid')
n_out.tail(1)
| node_list | node_pga | node_pgv | pid | |
|---|---|---|---|---|
| 111 | 2 | 0.442408 | 1.804097 | 2 |
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
# jet=plt.cm.jet
nx.draw(G,pos,node_size=20,edge_color = 'k',node_color=n_out.node_pga,width=0.5,with_labels=False)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(n_out.node_pga.min(), n_out.node_pga.max()))
plt.colorbar(sm,label='PGA(g)',shrink=0.7)
plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Node_PGA.png', dpi = 600,bbox_inches='tight')
Min, Mean, Max={n_out.node_pga.min(),n_out.node_pga.mean(), n_out.node_pga.max()}
Min, Mean, Max
(0.42531495339438213, 0.447621579746236, 0.468126922903879)
#convert pgv(m/s) to pgv(in/s)
node_pgv=n_out.node_pgv*39.3701
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
# jet=plt.cm.jet
nx.draw(G,pos,node_size=20,edge_color = 'k',node_color=node_pgv,width=0.5,with_labels=False)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(node_pgv.min(), node_pgv.max()))
plt.colorbar(sm,label='PGV(in/s)',shrink=0.7)
# plt.scatter(ex, ey, s=1000, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Node_PGV.png', dpi = 600,bbox_inches='tight')
dp = []
for index, row in link.iterrows():
stt = str(row['start_node'])+str(row['end_node']).rjust(10)+str(row['PipID']).rjust(10)
dp.append(stt)
# calculate intensity for links
pga_link1, pgv_link1 = Distance.pga_for_link(p1,node,ex,ey,M)
pga_link2, pgv_link2 = Distance.pga_for_link(p2,node,ex,ey,M)
p1_id = list(p1.PipID)
pga_link1 = list(pga_link1)
pgv_link1 = list(pgv_link1)
pf = {'Id':p1_id,'pga_link1':pga_link1,'pgv_link1':pgv_link1}
pipe_1=pd.DataFrame(pf)
pipe_1.set_index('Id',inplace=True);
p2_id = list(p2.PipID)
pga_link2 = list(pga_link2)
pgv_link2 = list(pgv_link2)
pf = {'Id':p2_id,'pga_link2':pga_link2,'pgv_link2':pgv_link2}
pipe_2=pd.DataFrame(pf)
pipe_2.set_index('Id',inplace=True);
p_pga=pipe_1.pga_link1.append(pipe_2.pga_link2)
p_pgv=pipe_1.pgv_link1.append(pipe_2.pgv_link2)
pipe_id=list(p_pga.index)
link_pga=list(p_pga)
link_pgv=list(p_pgv)
pf = {'Id':pipe_id,'link_pga':link_pga,'link_pgv':link_pgv}
pipe=pd.DataFrame(pf)
G = nx.parse_edgelist(dp, nodetype = int, data=(('id',int),))
edges= nx.get_edge_attributes(G,'id')
edge_list=list(edges.values())
pipe.head(2)
| Id | link_pga | link_pgv | |
|---|---|---|---|
| 0 | 119 | 0.466295 | 2.063929 |
| 1 | 121 | 0.463834 | 2.035097 |
pk = {'el':edge_list}
edges=pd.DataFrame(pk)
edges.head(2)
| el | |
|---|---|
| 0 | 119 |
| 1 | 222 |
p_out=edges.merge(pipe, left_on='el', right_on='Id')
p_out.head(2)
| el | Id | link_pga | link_pgv | |
|---|---|---|---|---|
| 0 | 119 | 119 | 0.466295 | 2.063929 |
| 1 | 222 | 222 | 0.467159 | 2.073992 |
pipe_pga=p_out.link_pga
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
jet=plt.cm.jet
nx.draw(G,pos,node_size=5,edge_color =pipe_pga, node_color='k',width=1.5,with_labels=False)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(pipe_pga.min(), pipe_pga.max()))
plt.colorbar(sm,label='PGA(g)',shrink=0.7)
# plt.scatter(ex, ey, s=1000, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Link_PGA.png', dpi = 600,bbox_inches='tight')
/Users/Ram/anaconda3/lib/python3.7/site-packages/networkx/drawing/nx_pylab.py:569: MatplotlibDeprecationWarning: The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead. and cb.iterable(edge_color) \ /Users/Ram/anaconda3/lib/python3.7/site-packages/networkx/drawing/nx_pylab.py:579: MatplotlibDeprecationWarning: The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead. for c in edge_color]):
pipe_pgv=p_out.link_pgv*39.3701
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
jet=plt.cm.jet
nx.draw(G,pos,node_size=5,edge_color =pipe_pgv, node_color='k',width=1.5,with_labels=False)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(pipe_pgv.min(), pipe_pgv.max()))
plt.colorbar(sm,label='PGV(in/s)',shrink=0.7)
plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Link_PGV.png', dpi = 600,bbox_inches='tight')
Min, Mean, Max={pipe_pgv.min(),pipe_pgv.mean(), pipe_pgv.max()}
Min, Mean, Max
(64.94287164997769, 73.20049117037568, 81.65328464383273)
r_rate=0.00187*pipe_pgv
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
jet=plt.cm.jet
nx.draw(G,pos,node_size=5,edge_color =r_rate, node_color='k',width=1.5,with_labels=False)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(r_rate.min(), r_rate.max()))
plt.colorbar(sm,label='Repair Rate/1000 ft',shrink=0.7)
# plt.scatter(ex, ey, s=1000, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_repair_rate.png', dpi = 600,bbox_inches='tight')
value=r_rate
Min, Mean, Max={value.min(),value.mean(), value.max()}
Min, Mean, Max
(0.12144316998545827, 0.1526916422839672, 0.13688491848860246)
Label Plot
d = []
for index, row in link.iterrows():
stt = str(row['start_node'])+str(row['end_node']).rjust(10)+str(row['PipID']).rjust(10)
d.append(stt)
d = []
for index, row in link.iterrows():
stt = str(row['start_node'])+str(row['end_node']).rjust(10)+str(row['PipID']).rjust(10)
# stt = str(row['start node'])+str(row['end node']).rjust(10)
d.append(stt)
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
edges= nx.get_edge_attributes(G,'id')
uu=edges.values()
el = list(uu)
# dl= list(fi)
# kl= list(zi)
pk = {'el':el}
T4=pd.DataFrame(pk)
T4.tail(1)
| el | |
|---|---|
| 125 | 245 |
T4.shape
(126, 1)
mk = list(pga)
ind= list(pos)
pk = {'ds':mk,'pid':ind}
T1=pd.DataFrame(pk)
# nP.set_index('lp1',inplace=True)
dam=T4.merge(T1, left_on='el', right_on='pid')
dam.tail(1)
| el | ds | pid | |
|---|---|---|---|
| 72 | 223 | 0.439259 | 223 |
dam.shape
(73, 3)
# plt.rcParams['figure.figsize'] = [14, 8]
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
# G = nx.parse_edgelist(d,nodetype=int)
jet=plt.cm.jet
# nx.draw(G,pos,node_size=20,edge_color = T4.el,node_color='k',width=1.5,with_labels=False)
# edge_labels = nx.get_edge_attributes(G, 'id')
# nx.draw_networkx_edge_labels(G,pos,edge_labels=edge_labels,label_pos=0.5,alpha=1.0,font_size=8, font_color='k')
# sm = plt.cm.ScalarMappable(norm=plt.Normalize(0, 4))
# plt.colorbar(sm,label='label',shrink=0.5)
# plt.scatter(ex, ey, s=1000, c='r', marker='*', zorder=2)
# plt.savefig('Fairfield_edge_label.png', dpi = 600,bbox_inches='tight')
pipe_characteristics = pd.read_csv('pipe_fairfield.csv',dtype={'index':str})
pipe_characteristics.set_index('PipID',inplace=True)
pipe_characteristics.head(2)
| start_node | end_node | type | Length | dia | M_type | Material | soil_type | age | |
|---|---|---|---|---|---|---|---|---|---|
| PipID | |||||||||
| 119 | 106 | 107 | 1 | 570 | 24 | CI | CIP | L | 1920 |
| 121 | 107 | 108 | 1 | 240 | 16 | CI | CIP | L | 1920 |
C=sd.cf.correction_factor(pipe_characteristics)
link_id=list(C.index)
link_C=list(C)
pf = {'Id':pipe_id,'C':link_C}
c_list=pd.DataFrame(pf)
c_list.tail(2)
| Id | C | |
|---|---|---|
| 124 | 2264 | 0.75 |
| 125 | 150 | 0.75 |
G = nx.parse_edgelist(dp, nodetype = int, data=(('id',int),))
edges= nx.get_edge_attributes(G,'id')
edge_list=list(edges.values())
pk = {'el':edge_list}
c_adj=pd.DataFrame(pk)
c_adj.tail(2)
| el | |
|---|---|
| 124 | 2264 |
| 125 | 245 |
c_fact=c_adj.merge(c_list, left_on='el', right_on='Id')
c_fact.tail(2)
| el | Id | C | |
|---|---|---|---|
| 124 | 2264 | 2264 | 0.75 |
| 125 | 245 | 245 | 0.75 |
c_color=plt.cm.jet
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =c_fact.C, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(c_fact.C.min(), c_fact.C.max()),cmap = c_color)
plt.colorbar(sm,label='Correction Factor',shrink=0.6)
# plt.scatter(ex, ey, s=1000, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Correc_Factor.png', dpi = 600,bbox_inches='tight')
nn=c_adj.merge(link, left_on='el', right_on='PipID')
rkk=nn.merge(p_out, left_on='el', right_on='Id')
out=rkk.merge(c_fact, left_on='el_x', right_on='el')
out.head(2)
| el_x | start_node | end_node | type | PipID | Length | dia | M_type | Material | soil_type | age | link_m | el_y | Id_x | link_pga | link_pgv | el | Id_y | C | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 119 | 106 | 107 | 1 | 119 | 570 | 24 | CI | CIP | L | 1920 | 245.036884 | 119 | 119 | 0.466295 | 2.063929 | 119 | 119 | 1.00 |
| 1 | 222 | 180 | 106 | 2 | 222 | 420 | 24 | DI | DCIP | H | 1960 | 152.826078 | 222 | 222 | 0.467159 | 2.073992 | 222 | 222 | 0.75 |
out.Length.sum()
58733
out.link_m.sum()/1000
20.417016241387323
out.link_m.min()/100
0.00039547885006946886
out.C.min()
0.75
out.link_pgv.min()
1.6495480491534869
out.link_m.min()*out.C.min()*out.link_pgv.min()/1000
4.892710242101675e-05
Mod_RR=0.00187*out.link_pgv*39.3701*out.C
edge_value=Mod_RR
c_color=plt.cm.jet
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()))
plt.colorbar(sm,label='Reair Rate/1000 ft',shrink=0.7)
# plt.scatter(ex, ey, s=1000, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Mod_RR.png', dpi = 600,bbox_inches='tight')
Min, Mean, Max={Mod_RR.min(),Mod_RR.mean(), Mod_RR.max()}
Min, Mean, Max
(0.09361900204720983, 0.2988714939973357, 0.42532195253446836)
C_PGV_L=(out.C)*(out.link_pgv*39.3701)*(out.link_m*3.28084/1000)
edge_value=C_PGV_L
c_color=plt.cm.rainbow
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='Θ',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_C_PGV_L.png', dpi = 600,bbox_inches='tight')
Min, Mean, Max={C_PGV_L.min(),C_PGV_L.mean(), C_PGV_L.max()}
Min, Mean, Max
(0.007501192875703854, 84.02404975329594, 391.1217009986187)
edge_value=2020-out.age
c_color=plt.cm.rainbow
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='Age (Year)',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Age.png', dpi = 600,bbox_inches='tight')
P_leak=1-np.exp(-0.00187*C_PGV_L)
edge_value=P_leak
c_color=plt.cm.jet
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='P(leak)',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_P_leak.png', dpi = 600,bbox_inches='tight')
P_break=P_leak*0.25
edge_value=P_break
c_color=plt.cm.nipy_spectral
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='P(break)',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_P_break.png', dpi = 600,bbox_inches='tight')
np.random.seed(10000)
node_FC=sd.FragilityCurve()
node_FC.add_state('Minor', 1, {'Default': lognorm(0.60, scale=0.30)})
node_FC.add_state('Moderate', 2, {'Default': lognorm(0.60, scale=0.70)})
node_FC.add_state('Extensive', 3, {'Default': lognorm(0.60, scale=1.25)})
node_FC.add_state('Complete', 4, {'Default': lognorm(0.65, scale=1.60)})
plot_fragility_curve(node_FC, xlabel='C∙PGV∙L', ylabel='Prob. of Failure',fill=False,xmax=1.5)
plt.xlabel('PGA(g)',size=20)
plt.ylabel('Prob. of Failure',size=20)
plt.legend(prop={'size':16},loc='upper left')
plt.plot(linewidth=2.0)
plt.rc('font',family='Times New Roman')
plt.rc('xtick', labelsize=16)
plt.rc('ytick', labelsize=16)
plt.grid(False)
plt.savefig('Fragility_Tanks.png', dpi = 600,bbox_inches='tight')
r, pga, pgv, pos = Distance.com_pga_dist(node,ex,ey,M)
node_id = list(node.id)
pga_val = list(pga)
dn1 = {'Id':node_id,'pga1':pga_val}
df1=pd.DataFrame(dn1)
df1.set_index('Id',inplace=True);
df1.head(2)
| pga1 | |
|---|---|
| Id | |
| 106 | 0.468127 |
| 107 | 0.463512 |
np.random.seed(12345)
node_Pr = node_FC.cdf_probability(df1['pga1'])
node_damage_state = node_FC.sample_damage_state(node_Pr)
node_damage_state_map = node_FC.get_priority_map()
node_damage_val = node_damage_state.map(node_damage_state_map)
node_damage_val.value_counts()
1 58 0 32 2 19 4 2 3 1 dtype: int64
node_Pr.tail(2)
| Minor | Moderate | Extensive | Complete | |
|---|---|---|---|---|
| Id | ||||
| 1 | 0.758044 | 0.23819 | 0.0466246 | 0.0268548 |
| 2 | 0.741319 | 0.222211 | 0.0417159 | 0.0239794 |
node_damage_val;
edge_value=node_damage_val
c_color=plt.cm.cool
nx.draw(G,pos,node_size=25,edge_color ='k', node_color=edge_value, width=0.5,with_labels=False,cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='Pf')
plt.savefig('Fairfield_Source_Damage.png', dpi = 600,bbox_inches='tight')
pipe_FC = sd.FragilityCurve()
pipe_FC.add_state('Break', 2, {'Default': expon(scale=1/(0.00187*.25))})
pipe_FC.add_state('Leak', 1, {'Default': expon(scale=1/(0.00187))})
plot_fragility_curve(pipe_FC, xlabel='Θ', ylabel='Prob. of Failure',fill=False,xmax=1000)
plt.xlabel('Θ',size=20)
plt.ylabel('Prob. of Failure',size=20)
plt.legend(prop={'size':16},loc='upper left')
plt.rc('font',family='Times New Roman')
plt.rc('xtick', labelsize=16)
plt.rc('ytick', labelsize=16)
plt.grid(False)
plt.savefig('Fragility.png', dpi = 600,bbox_inches='tight')
pga_link,pgv_link = Distance.pga_for_link(link,node,ex,ey,M)
p11_id = list(link.PipID)
p11_pga=list(pgv_link)
p11_len=list(out.link_m)
CV=list(C_PGV_L)
pk = {'Id':p11_id,'pgv':p11_pga,'len':p11_len,'CV':C_PGV_L}
pip_f=pd.DataFrame(pk)
pip_f.set_index('Id',inplace=True)
pip_f.head()
| pgv | len | CV | |
|---|---|---|---|
| Id | |||
| 119 | 2.063929 | 245.036884 | 65.324751 |
| 121 | 2.035097 | 152.826078 | 30.705590 |
| 122 | 2.030816 | 87.925988 | 23.112874 |
| 123 | 2.023753 | 56.541191 | 14.831570 |
| 124 | 2.014891 | 76.756730 | 20.064369 |
np.random.seed(12345)
pipe_Pr = pipe_FC.cdf_probability(pip_f['CV'])
pipe_damage_state = pipe_FC.sample_damage_state(pipe_Pr)
pipe_damage_state_map = pipe_FC.get_priority_map()
pipe_damage_val = pipe_damage_state.map(pipe_damage_state_map)
pipe_damage_val.value_counts()
0 112 1 10 2 4 dtype: int64
# plt.rcParams['figure.figsize'] = [5.5, 5]
edge_value=pipe_damage_val
c_color=plt.cm.rainbow
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
red_patch = mpatches.Patch(color='red', label='Break')
cyan_patch = mpatches.Patch(color='mediumspringgreen', label='Leak')
blue_patch = mpatches.Patch(color='blueviolet', label='None')
plt.legend(handles=[red_patch,cyan_patch, blue_patch],loc=10,fontsize=13,bbox_to_anchor=(0.8, 0.0, 0.1, 0.4))
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=1)
plt.savefig('Fairfield_Damage_States.png', dpi = 600,bbox_inches='tight')
bc_edge=list(nx.edge_betweenness_centrality(G).values())
ebc=pd.Series(bc_edge)
edge_value=ebc
c_color=plt.cm.OrRd
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='Betweenness Centrality',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_EBC.png', dpi = 600,bbox_inches='tight')
bc_node=list(nx.betweenness_centrality(G).values())
nbc=pd.Series(bc_node)
edge_value=nbc
c_color=plt.cm.Reds
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=20,edge_color ='k', node_color=edge_value,width=.5,with_labels=False,cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='Betweenness Centrality',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_NBC.png', dpi = 600,bbox_inches='tight')
ebc=pd.Series(bc_edge)
norm_ebc=(ebc-ebc.min())/(ebc.max()-ebc.min())
edge_value=norm_ebc
c_color=plt.cm.OrRd
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='Normalized EBC',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_Norm_EBC.png', dpi = 600,bbox_inches='tight')
pc_id = list(link.PipID)
pc_age=list(0.0003*(2020-out.age)**2-0.0003*(2020-out.age)+1)
pk = {'Id':pc_id,'CI':pc_age}
pc_f=pd.DataFrame(pk)
pc_f.set_index('Id',inplace=True)
pc_f.head()
| CI | |
|---|---|
| Id | |
| 119 | 3.970 |
| 121 | 2.062 |
| 122 | 3.970 |
| 123 | 3.970 |
| 124 | 2.896 |
edge_value=pc_f.CI
c_color=plt.cm.jet
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='Condition Index',shrink=0.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_CI.png', dpi = 600,bbox_inches='tight')
Min, Mean, Max={pc_f.CI.min(),pc_f.CI.mean(), pc_f.CI.max()}
Min, Mean, Max
(1.7349999999999999, 3.5346904761904856, 4.5969999999999995)
pc_id = list(link.PipID)
pc_id = list(link.PipID)
CI=list(pc_f.CI)
nbc=list(norm_ebc)
pk = {'Id':pc_id,'CI':CI,'nbc':norm_ebc}
cost=pd.DataFrame(pk)
cost.set_index('Id',inplace=True)
Index=cost.CI+cost.nbc
Index=cost.CI+cost.nbc
pc_id = list(link.PipID)
ind=list(Index)
dl=list(pipe_damage_val)
pk = {'Id':pc_id,'ind':ind,'dl':dl}
need=pd.DataFrame(pk)
need.set_index('Id',inplace=True)
need.head()
| ind | dl | |
|---|---|---|
| Id | ||
| 119 | 4.570193 | 0 |
| 121 | 2.687470 | 0 |
| 122 | 4.010718 | 0 |
| 123 | 4.534392 | 0 |
| 124 | 2.936718 | 0 |
mmm=need[(need.ind>3.5)&(need.dl>0)]
link['link_m']=list(Length(node,link))
need.dl.value_counts()
0 112 1 10 2 4 Name: dl, dtype: int64
# data1 and data2 should be node and link, respectively
def C_Check(data):
maint=[]
for index,row in data.iterrows():
if (row['ind']>=3.5) & (row['dl']>0):
mm=2
elif (row['ind']<3.5) & (row['dl']>0):
mm=1
else:
mm=0
maint.append(mm)
maint=np.array(maint)
return maint
need['MA']=list(C_Check(need))
need.MA.value_counts()
0 112 2 13 1 1 Name: MA, dtype: int64
need.head(2)
| ind | dl | MA | |
|---|---|---|---|
| Id | |||
| 119 | 4.570193 | 0 | 0 |
| 121 | 2.687470 | 0 | 0 |
edge_value=need.MA
c_color=plt.cm.jet
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=5,edge_color =edge_value, node_color='k',width=1.5,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
red_patch = mpatches.Patch(color='firebrick', label='Replace')
cyan_patch = mpatches.Patch(color='mediumspringgreen', label='Repair')
blue_patch = mpatches.Patch(color='navy', label='None')
plt.legend(handles=[red_patch,cyan_patch, blue_patch],loc=10,fontsize=13,bbox_to_anchor=(0.8, 0.0, 0.1, 0.4))
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
plt.savefig('Fairfield_MA.png', dpi = 600,bbox_inches='tight')
need[need.MA>=1].to_csv('pipe_repair.csv')
nw=need.merge(out, left_on='Id', right_on='el_x')
nw['NoB']=(0.00187*nw.link_pgv*39.3701*nw.C)*(nw.link_m*3.28084/1000)
nw.NoB;
nw['nNoB']=np.ceil(nw.NoB).astype(int)
def pipe_annual_cost(data,pipe_cost=None):
"""
"""
nrpr = []
nrpl = []
network_cost = 0
if pipe_cost is None:
diameter = [4, 6, 8, 10, 12, 14, 16, 18, 20, 24, 28, 30, 32, 34, 36] # inch
rpl = [600, 630, 675, 750, 825, 1200, 1950, 2400, 2700, 3450, 4350, 4650, 5250, 5700, 6300] # replace cost/m
rpr = [400, 420, 450, 500, 550, 800, 1300, 1600, 1800, 2300, 2900, 3100, 3500, 3800, 4200] # repair cost
# diameter = np.array(diameter)*0.0254 # m
repair_cost = pd.Series(rpr,diameter)
replace_cost = pd.Series(rpl,diameter)
# Pipe construction cost
for index, row in data.iterrows():
dia = row['dia']
length=row['link_m']
idxrpr = np.argmin([np.abs(repair_cost.index - dia)])
idxrpl = np.argmin([np.abs(replace_cost.index - dia)])
#print(link_name, pipe_cost.iloc[idx], link.length)
repair_C = network_cost + repair_cost.iloc[idxrpr]
replace_C = network_cost + replace_cost.iloc[idxrpl]*length
nrpr.append(repair_C)
nrpl.append(replace_C)
nrpr = np.array(nrpr)
nrpl = np.array(nrpl)
return nrpr,nrpl
xxx,yyy=pipe_annual_cost(nw)
nw['repair_C']=xxx
nw['replace_C']=yyy
def C_Est(data):
cost=[]
for index,row in data.iterrows():
if (row['MA']==1):
mm=row['repair_C']*row['nNoB']
elif (row['MA']==2):
mm=row['replace_C']
else:
mm=0
cost.append(mm)
cost=np.array(cost)
return cost
nw['T_Cost']=list(C_Est(nw))
Cost_for_repair=print(nw[nw.MA==1].T_Cost.sum())
2300.0
Cost_for_replace=print(nw[nw.MA==2].T_Cost.sum())
6370785.43041506
total_maintenace_cost=print(nw.T_Cost.sum())
6373085.43041506
nw[nw.MA>0].link_m.sum()/1000
2.3209501209846217
#Per_Cost_graph
# diameter = [4, 6, 8, 10, 12, 14, 16, 18, 20, 24, 28, 30, 32, 34] # inch
# rpr = [95, 120, 140, 165, 195, 220, 255, 385, 455, 525, 563, 600, 638, 638] # repair cost
# rpl = [1575, 1800, 2270, 2700, 2800, 3000, 4225, 4325, 4450, 4450, 4450, 4450, 4450, 4450] # replace cost
# repair_cost = pd.Series(rpr,diameter)
# replace_cost = pd.Series(rpl,diameter)
# xp=replace_cost.index
# yp1=replace_cost.values
# yp2=repair_cost.values
# plt.plot(xp,yp1)
# plt.plot(xp,yp2)
# plt.xlabel('Diameter(in)',size=16)
# plt.ylabel('Cost $',size=16)
# plt.rc('font',family='Times New Roman')
# plt.savefig('Cost_Graph.png', dpi = 600,bbox_inches='tight')
Costt=nw.T_Cost
edge_value=Costt
c_color=plt.cm.rainbow
G = nx.parse_edgelist(d, nodetype = int, data=(('id',int),))
nx.draw(G,pos,node_size=0.5,edge_color =edge_value, node_color='k',width=2,with_labels=False,edge_cmap = c_color)
sm = plt.cm.ScalarMappable(norm=plt.Normalize(edge_value.min(), edge_value.max()),cmap = c_color)
plt.colorbar(sm,label='$USD',shrink=.7)
# plt.scatter(ex, ey, s=500, c='r', marker='*', zorder=2)
# plt.savefig('Fairfield_Cost.png', dpi = 600,bbox_inches='tight')
<matplotlib.colorbar.Colorbar at 0x1c220442b0>
# nx.shortest_path(G, 1721);
# pd.Series(nx.single_source_shortest_path(G,1721));
n_link=nw[nw.MA==0]
pp = []
for index, row in n_link.iterrows():
stt = str(row['start_node'])+str(row['end_node']).rjust(10)
pp.append(stt)
# draw Graph of water network
G = nx.parse_edgelist(pp,nodetype=int)
nx.draw(G,pos,node_size=5, node_color='k', edge_color='b',width=1.5)
nx.draw(G3,pos3,node_size=60,node_color='g',with_labels=False)
plt.savefig('Isolated_Fairfield.png', dpi = 600,bbox_inches='tight')
nx.number_of_isolates(G)
0
# read node and link positions
node = pd.read_csv('node_deg_fairfield.csv')
# link = pd.read_csv('edge classify.csv')
link = pd.read_csv('pipe_fairfield.csv')
import os
import folium
from folium.plugins.measure_control import MeasureControl
link.head(1)
| start_node | end_node | type | PipID | Length | dia | M_type | Material | soil_type | age | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 106 | 107 | 1 | 119 | 570 | 24 | CI | CIP | L | 1920 |
# data1 and data2 should be node and link, respectively
def Coord(data1,data2):
pos={}
ll=[]
for index, row in data2.iterrows():
sp=data1[data1.id==row['start_node']]
start_x, start_y = (list(sp.x),list(sp.y))
ep=data1[data1.id==row['end_node']]
end_x,end_y=(list(ep.x),list(ep.y))
pos[int(row['PipID'])]=start_y+start_x,end_y+end_x
ll=data2['Length']
return pos, ll
posL,ll=Coord(node,link)
coordinatesL=posL.values()
coL=list(coordinatesL)
bc_node=nx.betweenness_centrality(G)
bc_node;
bc_edge=nx.edge_betweenness_centrality(G)
bc_edge;
link_attribute=bc_edge
link_width=2
link_range=[0,1]
link_cmap_bins='cut'
link_cmap=['cornflowerblue', 'forestgreen', 'gold', 'firebrick']
if node_attribute is not None:
if isinstance(node_attribute, list):
node_cmap=['red']
node_attribute = node_v
node_attribute = pd.Series(node_attribute)
if node_range[0] is not None:
node_attribute[node_attribute < node_range[0]] = node_range[0]
if node_range[1] is not None:
node_attribute[node_attribute > node_range[1]] = node_range[1]
if node_cmap_bins == 'cut':
node_colors, node_bins = pd.cut(node_attribute, len(node_cmap),
labels=node_cmap, retbins =True)
elif node_cmap_bins == 'qcut':
node_colors, node_bins = pd.qcut(node_attribute, len(node_cmap),
labels=node_cmap, retbins =True)
if link_attribute is not None:
if isinstance(link_attribute, list):
link_cmap=['red']
link_attribute = link_v
link_attribute = pd.Series(link_attribute)
if link_range[0] is not None:
link_attribute[link_attribute < link_range[0]] = link_range[0]
if link_range[1] is not None:
link_attribute[link_attribute > link_range[1]] = link_range[1]
if link_cmap_bins == 'cut':
link_colors, link_bins = pd.cut(link_attribute, len(link_cmap),
labels=link_cmap, retbins =True)
elif link_cmap_bins == 'qcut':
link_colors, link_bins = pd.qcut(link_attribute, len(link_cmap),
labels=link_cmap, retbins =True)
center = pd.DataFrame(pos).mean(axis=1)
m = folium.Map(location=[center.iloc[1], center.iloc[0]], zoom_start=14)
folium.TileLayer('cartodbpositron').add_to(m)
folium.TileLayer('Stamen Toner').add_to(m)
folium.TileLayer('openstreetmap').add_to(m)
folium.TileLayer('stamenterrain').add_to(m)
folium.TileLayer('stamenwatercolor').add_to(m)
folium.TileLayer('cartodbdark_matter').add_to(m)
if node_size > 0:
for index, row in node.iterrows():
loc = (row['y'], row['x'])
radius = node_size
color = 'black'
if node_attribute is not None:
if index in node_attribute:
color = node_colors[index]
else:
radius = 0.1
folium.CircleMarker(loc, color=color, fill=True,fill_color=color, radius=1, fill_opacity=0.7, opacity=0.7).add_to(m)
if link_width > 0:
for index, row in link.iterrows():
sp=node[node.id==row['start_node']]
start_x,start_y=(list(sp.x),list(sp.y))
ep=node[node.id==row['end_node']]
end_x,end_y=(list(ep.x),list(ep.y))
posK=start_y+start_x,end_y+end_x
weight = link_width
color='black'
if link_attribute is not None:
if index in link_attribute:
color = link_colors[index]
else:
weight = 0.5
bins = list(link_v.quantile([0, 0.25, 0.5, 0.75, 1]))
folium.PolyLine(posK, color=color,fill=True,fill_color=color,weight=3, opacity=0.7).add_to(m)
# folium.Marker(location=[center.iloc[1], center.iloc[0]],popup='RKM',icon=folium.Icon(color='red')).add_to(m)
folium.LayerControl().add_to(m)
m