
import os

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.getcwd()

import numpy as np
import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame
from shapely.geometry import Point
import folium
import osmnx as ox
import networkx as nx
import math
import momepy
# from funcs_network_FK import roads_type_folium
from shapely import geometry
from shapely.geometry import Point, Polygon
import psycopg2
# import db_connect
import datetime
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from folium_stuff_FK_map_matching import plot_graph_folium_FK
import glob
from funcs_network_FK import cost_assignment
import statistics


# load grafo
# file_graphml = 'Catania__Italy.graphml'
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
# ox.plot_graph(grafo)
place_country = "Catania, Italy"
# cost_assignment(file_graphml, place_country)

####################################################
#### overlay Hexagonal grid  #######################


os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
# load all EDGES obtained from the map-matching algorithm
TIME_EDGES = gpd.read_file("TIME_EDGES.geojson")
# TIME_EDGES.plot()


# load all EDGES obtained from the map-matching algorithm (RECORDS and FREQUENCIES)
RECORDS_CATANIA = gpd.read_file("FREQUENCIES_and_RECORDS_by_EDGES.geojson")

### merge TIME_EDGES with RECORDS_CATANIA referred to the counts (records about number of veicles)
TIME_EDGES = pd.merge(TIME_EDGES, RECORDS_CATANIA, on=['u', 'v', 'DESTINATION', 'ORIGIN', 'access',
                                                'area', 'bridge', 'distance', 'hour',
                                                'id', 'junction', 'key', 'lanes', 'length', 'oneway', 'ref', 'service',
                                                'speed', 'time',  'timedate', 'track_ID', 'tunnel' , 'width',
                                                'length(km)'
                                                ], how='left')

TIME_EDGES=TIME_EDGES.rename(columns = {'geometry_x':'geometry'})
TIME_EDGES=TIME_EDGES.rename(columns = {'color_x':'color'})
TIME_EDGES.drop(['geometry_y'], axis=1, inplace = True)
TIME_EDGES.drop(['color_y'], axis=1, inplace = True)

## fill nans by values after non-nan values (for 'hour' and 'timedate')
TIME_EDGES['hour'] = TIME_EDGES['hour'].ffill()
TIME_EDGES['timedate'] = TIME_EDGES['timedate'].ffill()
## drop nan values from the column "hour"
TIME_EDGES.dropna(subset = ['hour'], inplace= True)

# save file
TIME_EDGES.to_file(filename='TIME_EDGES_complete.geojson', driver='GeoJSON')

# load the file assigning the type of highway to each (u,v) edge
edges_highways = pd.read_csv("edges_highways.csv")
# remove duplicates
edges_highways = edges_highways.drop_duplicates(['u', 'v'])
edges_highways = edges_highways[~edges_highways['highway'].isin(['residential', 'secondary_link', 'tertiary_link',
                                                     'footway', 'primary_link', 'secondary', 'tertiary', 'other'])]

MAIN_TIME_EDGES = pd.merge(TIME_EDGES, edges_highways, on=['u', 'v'], how='inner')
# MAIN_TIME_EDGES.plot()
# TIME_EDGES.plot()
# TIME_EDGES = MAIN_TIME_EDGES

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')


# define distance between two coordinates
def haversine(coord1, coord2):
    # Coordinates in decimal degrees (e.g. 43.60, -79.49)
    lon1, lat1 = coord1
    lon2, lat2 = coord2
    R = 6371000  # radius of Earth in meters
    phi_1 = np.radians(lat1)
    phi_2 = np.radians(lat2)
    delta_phi = np.radians(lat2 - lat1)
    delta_lambda = np.radians(lon2 - lon1)
    a = np.sin(delta_phi / 2.0) ** 2 + np.cos(phi_1) * np.cos(phi_2) * np.sin(delta_lambda / 2.0) ** 2
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    meters = R * c  # output distance in meters
    km = meters / 1000.0  # output distance in kilometers
    meters = round(meters)
    km = round(km, 3)
    # print(f"Distance: {meters} m")
    # print(f"Distance: {km} km")
    return meters


# define the extension of the EDGE file
xmin, ymin, xmax, ymax = TIME_EDGES.total_bounds  # lat-long of 2 corners
# East-West extent of Toronto = 42193 metres
EW = haversine((xmin, ymin), (xmax, ymin))
# North-South extent of Toronto = 30519 metres
NS = haversine((xmin, ymin), (xmin, ymax))
# diamter of each hexagon in the grid = 900 metres
d = 1800
# horizontal width of hexagon = w = d* sin(60)
w = d * np.sin(np.pi / 3)
# Approximate number of hexagons per row = EW/w
n_cols = int(EW / w) + 1
# Approximate number of hexagons per column = NS/d
n_rows = int(NS / d) + 10

# Make a hexagonal grid to cover the entire area
from matplotlib.patches import RegularPolygon

ax = TIME_EDGES.boundary.plot(edgecolor='black', figsize=(20, 60))
w = (xmax - xmin) / n_cols  # width of hexagon
d = w / np.sin(np.pi / 3)  # diameter of hexagon
array_of_hexes = []
for rows in range(0, n_rows):
    hcoord = np.arange(xmin, xmax, w) + (rows % 2) * w / 2
    vcoord = [ymax - rows * d * 0.75] * n_cols
    for x, y in zip(hcoord, vcoord):  # , colors):
        hexes = RegularPolygon((x, y), numVertices=6, radius=d / 2, alpha=0.2, edgecolor='k')
        verts = hexes.get_path().vertices
        trans = hexes.get_patch_transform()
        points = trans.transform(verts)
        array_of_hexes.append(Polygon(points))
        ax.add_patch(hexes)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])

hex_grid = gpd.GeoDataFrame({'geometry': array_of_hexes}, crs={'init': 'epsg:4326'})
# hex_grid = gpd.GeoDataFrame({'geometry':array_of_hexes},crs={'init':'epsg:3035'})
# hex_grid.plot()


#############################################################################################
# create basemap
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
#############################################################################################

# get hexagons containing the edges
egdes_hex = gpd.sjoin(hex_grid, TIME_EDGES, how='inner', op='intersects')
# egdes_hex.plot()
egdes_hex.geometry.to_file(filename='egdes_hex.geojson', driver='GeoJSON')
hex_grid.geometry.to_file(filename='hex_grid.geojson', driver='GeoJSON')

egdes_hex.reset_index(level=0, inplace=True)
hex_grid.reset_index(inplace=True)

style_hex = {'fillColor': '#00000000', 'color': '#00FFFFFF'}
style_edges = {'fillColor': '#00000000', 'color': '#000000'}
style_egdes_hex = {'fillColor': '#00000000', 'color': '#ff0000'}

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
folium.GeoJson('TIME_EDGES.geojson', style_function=lambda x: style_edges).add_to((my_map))
# folium.GeoJson('hex_grid.geojson',  style_function=lambda x: style_hex).add_to((my_map))
# folium.GeoJson('egdes_hex.geojson',  style_function=lambda x: style_egdes_hex).add_to((my_map))

# add 'u' and 'v' as highligths for each edge (in blue)
folium.GeoJson(
    # data to plot
    hex_grid[['index', 'geometry']].to_json(),
    show=True,
    style_function=lambda x: style_hex,
    highlight_function=lambda x: {'weight': 1,
                                  'color': 'yellow',
                                  'fillOpacity': 0.2
                                  },
    # fields to show
    tooltip=folium.features.GeoJsonTooltip(
        fields=['index']
    ),
).add_to(my_map)

folium.GeoJson(
    # data to plot
    egdes_hex[['index', 'geometry']].to_json(),
    show=True,
    style_function=lambda x: style_egdes_hex,
    highlight_function=lambda x: {'weight': 1,
                                  'color': 'yellow',
                                  'fillOpacity': 0.2
                                  },
    # fields to show
    tooltip=folium.features.GeoJsonTooltip(
        fields=['index']
    ),
).add_to(my_map)

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
my_map.save("hex_grid.html")

##############################################################################################
##############################################################################################
###### ORIGIN - DESTINATION TRIPS ############################################################

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')


###########################################
## get edges inside each hexagon ##########

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
# load hex_grid
hex_grid = gpd.read_file("hex_grid.geojson")

# find all edges (links) intersected by the hexagonal grid
## 'links' are the edges (u, v)
## an 'element' is the ensemble of links within a cell or intersecting a cell

elements = {}

hex_grid = hex_grid[hex_grid.index == 533]  # select only one particular cell (when using 1800 m diameter cell)
# hex_grid = hex_grid[hex_grid.index.isin([2359, 2290])]  # select two particular cells (only with 900m hexagon diameter)
for cells in hex_grid.itertuples(index=True):
    edges = []
    index_elements = []
    index_cells = []
    for links in TIME_EDGES.itertuples(index=True):
        if links.geometry.intersects(cells.geometry) is True:
            print("OK=========================================================================IN")
            index_elements.append(links.Index)
            index_cells.append(cells.Index)
            STREET = links.u, links.v
            edges.append(STREET)
            elements[cells.Index] = edges
        else:
            print("NOT====IN")


## filter gdf_edges based on index_edges (edges in the cell)
gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
# save as GeoJson file
gdf_elements.to_file(filename='gdf_elements.geojson', driver='GeoJSON')
gdf_elements.plot()
# single exagon (cell)
hex_gdf_edges = hex_grid[hex_grid.index.isin(index_cells)]
hex_grid.plot()

## make index as column
hex_gdf_edges.reset_index(inplace=True)
hex_grid.reset_index(inplace=True)

#############################################################################################
# create basemap
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
#############################################################################################

folium.GeoJson('gdf_elements.geojson', style_function=lambda x: style_edges).add_to((my_map))

style_edge = {'fillColor': '#00000000', 'color': '#ff0000'}
# add 'u' and 'v' as highligths for each edge (in blue)
folium.GeoJson(
    # data to plot
    hex_gdf_edges[['index', 'geometry']].to_json(),
    show=True,
    style_function=lambda x: style_edge,
    highlight_function=lambda x: {'weight': 1,
                                  'color': 'red',
                                  'fillOpacity': 1
                                  },
    # fields to show
    tooltip=folium.features.GeoJsonTooltip(
        fields=['index']
    ),
).add_to(my_map)

## save map
my_map.save("elements_and_links.html")


'''
# load grafo of all Catania
file_graphml = 'Catania__Italy.graphml'
grafo = ox.load_graphml(file_graphml)
ox.plot_graph(grafo)
# get edges and nodes
gdf_nodes, gdf_edges = ox.graph_to_gdfs(grafo)
# save edges as GeoJson file
# gdf_edges.to_file(filename='Catania_edges_60km.geojson', driver='GeoJSON')
gdf_edges.to_csv('Catania_edges_60km.csv')
'''

###################################################################################
## calculate the initial travel timeS and shortest paths between the OD pairs #####
###################################################################################


gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
OD = gdf_elements.drop_duplicates(['ORIGIN', 'DESTINATION'])
len(OD)
OD.reset_index(level=0, inplace=True)
# select one OD pair
i = 1
O = OD.ORIGIN.iloc[i]
D = OD.DESTINATION.iloc[i]
# print('ORIGIN:', O, 'DESTINATION:', D, '================================================')
all_trips_OD = TIME_EDGES[(TIME_EDGES['ORIGIN'] == O) & (TIME_EDGES['DESTINATION'] == D)]
all_trips_OD.plot()  # these are all the ONLY elements intersecting the two cells 2359, 2290
all_trips_OD.geometry.to_file(filename='all_trips_OD.geojson', driver='GeoJSON')


# init_OD_travel_time = sum(all_trips_OD['travel time (min)'])   # not good!!
### Initial OD shortest path
init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length') # using length
path_OD_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path[1:]))
# ox.plot_graph_route(grafo, init_shortest_OD_path, route_color='green', fig_height=12, fig_width=12)
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path.save('OD_shortest_path.html')
## isolate edges in the grafo from 'path_OD_edges'
if len(path_OD_edges) > 1:
    # isolate edges in the grafo from 'all_matched_edges'
    df_nodes = []
    for i in range(len(path_OD_edges)):
        # print(path_OD_edges[i])
        route = path_OD_edges[i]
        df_nodes.append(route)

    df_nodes = pd.DataFrame(df_nodes)
    df_nodes.columns = ['u', 'v']

    # filter TIME_EDGES with df_nodes
    keys = list(df_nodes.columns.values)
    index_TIME_EDGES = TIME_EDGES.set_index(keys).index
    index_df_nodes = df_nodes.set_index(keys).index
    edges_matched_route_OD = TIME_EDGES[index_TIME_EDGES.isin(index_df_nodes)]
    edges_matched_route_OD.plot()
    edges_matched_route_OD.geometry.to_file(filename='edges_matched_route_OD.geojson', driver='GeoJSON')

    ########################################################
    ### Initial OD shortest TIME (based on VIASAT data #####
    ########################################################
    init_OD_travel_time = sum(edges_matched_route_OD['travel time (sec)'])


   # create basemap
    ave_LAT = 37.53988692816245
    ave_LON = 15.044971594798902
    my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
    #############################################################################################

    folium.GeoJson('edges_matched_route_OD.geojson', style_function=lambda x: style_edges).add_to((my_map))
    ## save map
    my_map.save("edges_matched_route_OD.html")

##################################################################################
##################################################################################

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')

# load grafo
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
# ox.plot_graph(grafo)

# make attr "cost" as float!!! (units of "cost" are SECONDS!)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['cost']) >0:
        attr['cost'] = float(attr.get("cost"))


########################################################################################
## replace "cost" with "travel time" got from VIASAT analysis for  given u,v pair ######
########################################################################################
#### LONG TIME RUNNING #################################################################

for U, V, key, attr in grafo.edges(keys=True, data=True):
    # print(attr["cost"])
    attr['VIASAT_cost'] = attr.get("cost")
    zipped = zip(list(TIME_EDGES.u), list(TIME_EDGES.v)) # we should use TIME_EDGES....
    if (U,V) in zipped:
        travel_time = TIME_EDGES[(TIME_EDGES.u == U) &
                                             (TIME_EDGES.v == V)]['travel time (sec)']  # in seconds
        print(float(travel_time))
        attr['VIASAT_cost'] = float(travel_time)
        grafo.add_edge(U, V, key, attr_dict=attr)
ox.save_graphml(grafo, filename= 'CATANIA_VIASAT_cost.graphml')

###########################################################################################################
###########################################################################################################
# create basemap
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=8.5, tiles='cartodbpositron')
###########################################################################################################
###########################################################################################################

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')

file_graphml = 'CATANIA_VIASAT_cost.graphml'
grafo = ox.load_graphml(file_graphml)


# make attr "cost" as float!!! (units of "cost" are SECONDS!)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost']) > 0:
        attr['VIASAT_cost'] = float(attr.get("VIASAT_cost"))
    if len(attr['cost']) >0:
        attr['cost'] = float(attr.get("cost"))

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')

init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length')  # using cost (time)
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path_edges = list(zip(init_shortest_OD_path,init_shortest_OD_path[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='cost')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path.save('OD_shortest_path_length.html')

## find shortest path based on the "cost" (time)
init_shortest_OD_path_cost = nx.shortest_path(grafo, O, D, weight='cost')  # using cost (time)
path_cost = ox.plot_route_folium(grafo, init_shortest_OD_path_cost, route_color='green')
path_edges = list(zip(init_shortest_OD_path,init_shortest_OD_path_cost[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='cost')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path_cost.save('OD_shortest_path_cost.html')

# add shortest path (by "cost" to the my_map
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']

# filter 'TIME_EDGES' with 'df_nodes'
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_cost = TIME_EDGES[index_gdf_edges.isin(index_df_nodes)]

## add plot in Folium map
# save first as geojson file
style_edge = {'fillColor': '#00000000', 'color': '#000000'}
edges_shortest_route_cost.geometry.to_file(filename='edges_shortest_route_cost.geojson',
                                 driver='GeoJSON')
folium.GeoJson('edges_shortest_route_cost.geojson', style_function=lambda x: style_edges).add_to((my_map))
my_map.save("shortest_routes_examples.html")


# find shortet path by "VIASAT cost" that is the "time" got from the VIASAT data
init_shortest_OD_path_VIASAT_cost = nx.shortest_path(grafo, O, D, weight='VIASAT_cost')  # using cost (time)
path_VIASAT_cost = ox.plot_route_folium(grafo, init_shortest_OD_path_VIASAT_cost, route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_cost,init_shortest_OD_path_VIASAT_cost[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='VIASAT_cost')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path_VIASAT_cost.save('OD_VIASAT_shortest_path_cost.html')
# total time from O--D at normal conditions
time_OD_NULL = lr/3600

# add shortest path (by "VIASAT cost" to the my_map
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']

# filter 'TIME_EDGES' with 'df_nodes'
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_cost = TIME_EDGES[index_gdf_edges.isin(index_df_nodes)]

## add plot in Folium map
# save first as geojson file
style_VIASAT = {'fillColor': '#00000000', 'color': '#ff0000'}
edges_shortest_route_VIASAT_cost.geometry.to_file(filename='edges_shortest_route_VIASAT_cost.geojson',
                                 driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_cost.geojson', style_function=lambda x: style_VIASAT).add_to((my_map))
my_map.save("shortest_routes_examples.html")


# get an edge that is on the same route....
U = 2811353264
V = 2811352855


os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
# close a LINK (u,v pair) by adding sufficiently large penalty M (time in seconds ~ 5 hours = 18000 secs)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost_penalty'] = attr.get("VIASAT_cost")
    if u == U and v == V:
        penalty = 18000  # seconds of closure of the link (u,v)
        print("gotta!=============================================================================================")
        # break
        attr['VIASAT_cost_penalty'] = float(attr['VIASAT_cost']) + penalty
        grafo.add_edge(u, v, key, attr_dict=attr)
ox.save_graphml(grafo, filename= 'CATANIA_VIASAT_cost_penalty.graphml')


## load grafo again...
file_graphml = 'CATANIA_VIASAT_cost_penalty.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost_penalty']) > 0:
        attr['VIASAT_cost_penalty'] = float(attr.get("VIASAT_cost_penalty"))


# check...
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if u == U and v == V:
        print("gotta!=============================================================================================")
        print(attr['VIASAT_cost_penalty'])
        break

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
# get shortest path again...
init_shortest_OD_path_VIASAT_penalty = nx.shortest_path(grafo, O, D, weight='VIASAT_cost_penalty')  # using cost (time)
path = ox.plot_route_folium(grafo, init_shortest_OD_path_VIASAT_penalty, route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_penalty,init_shortest_OD_path_VIASAT_penalty[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='VIASAT_cost_penalty')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path.save('OD_VIASAT_shortest_path_penalty.html')
# total time from O--D with penalty
time_OD_penalty = lr/3600


gdf_nodes, gdf_edges = ox.graph_to_gdfs(grafo)
# add shortest path (by "VIASAT cost" to the my_map
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']

# filter 'TIME_EDGES' with 'df_nodes'
keys = list(df_nodes.columns.values)
index_gdf_edges = gdf_edges.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_penalty = gdf_edges[index_gdf_edges.isin(index_df_nodes)]

## add plot in Folium map
# save first as geojson file
style_VIASAT_penalty = {'fillColor': '#00000000', 'color': '#0000ff'}
edges_shortest_route_VIASAT_penalty.geometry.to_file(filename='edges_shortest_route_VIASAT_penalty.geojson',
                                 driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_penalty.geojson', style_function=lambda x: style_VIASAT_penalty).add_to((my_map))
my_map.save("shortest_routes_examples.html")


#################################################################################
#################################################################################
#################################################################################
#################################################################################

#############################
### closure IMPACT ##########
#############################

## flux [N vehicles/hour] per OD pair (travel demand per unit time)
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')

df_flux = edges_matched_route_OD
df_flux['timedate'] = pd.to_datetime(df_flux['timedate'])

# print(df_flux.timedate.astype('<M8[h]'))
# BBB = (df_flux.groupby([df_flux.timedate.astype('<M8[h]')]).mean()).records
# GGG = (df_flux.groupby([df_flux.hour.astype('<M8[h]')]).mean())
# BBB = pd.DataFrame(BBB)
# TRAVEL_DEMAND_OD = BBB.records.mean()

## use 15 minutes average to calculate FLUX [vehicles/hour]
df_flux = df_flux.resample(rule='15Min', on='timedate').mean()
df_flux['records'] = df_flux['records'].ffill()  # 15 minutes flux
df_flux['records'] = df_flux['records']/0.25
BBB = pd.DataFrame(df_flux['records'])
TRAVEL_DEMAND_OD = BBB.records.mean()  # veichles/hour

## difference time between closure and normal conditions
DT = time_OD_penalty - time_OD_NULL  # (hours)

## closure impact
if DT < penalty/3600:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD*DT*( (penalty/3600) - DT/2)  # hours
else:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD* (((penalty/3600)**2)/2)  # hours



###########################################################
#### save all the elements into a dictionary  #############
###########################################################
elements = {}

# get elements (all edges) within each hexagonal cell)
for idx in range(len(hex_grid)):
    print(idx)
    # hex_grid = hex_grid[hex_grid.index.isin([533])]  # select one particular cells
    hex_grid_idx = hex_grid[hex_grid.index.isin([idx])]  # select one cells by time
    for cells in hex_grid_idx.itertuples(index=True):
        edges = []
        index_elements = []
        index_cells = []
        for links in TIME_EDGES.itertuples(index=True):
            if links.geometry.intersects(cells.geometry) is True:
                print("OK=========================================================================IN")
                index_elements.append(links.Index)
                index_cells.append(cells.Index)
                STREET = links.u, links.v
                edges.append(STREET)
                # elements[cells.Index] = edges
    gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
    if len(gdf_elements) > 0:
        elements[cells.Index] = gdf_elements
        # gdf_elements.plot()


import numpy as np
# Save
np.save('elements.npy', elements)


e0dc5ce4325fe03afdcddce56d18384b9d31952b
<mask token>

<mask token>
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.getcwd()
<mask token>
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
<mask token>
TIME_EDGES.drop(['geometry_y'], axis=1, inplace=True)
TIME_EDGES.drop(['color_y'], axis=1, inplace=True)
<mask token>
TIME_EDGES.dropna(subset=['hour'], inplace=True)
TIME_EDGES.to_file(filename='TIME_EDGES_complete.geojson', driver='GeoJSON')
<mask token>
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')


def haversine(coord1, coord2):
    lon1, lat1 = coord1
    lon2, lat2 = coord2
    R = 6371000
    phi_1 = np.radians(lat1)
    phi_2 = np.radians(lat2)
    delta_phi = np.radians(lat2 - lat1)
    delta_lambda = np.radians(lon2 - lon1)
    a = np.sin(delta_phi / 2.0) ** 2 + np.cos(phi_1) * np.cos(phi_2) * np.sin(
        delta_lambda / 2.0) ** 2
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    meters = R * c
    km = meters / 1000.0
    meters = round(meters)
    km = round(km, 3)
    return meters


<mask token>
for rows in range(0, n_rows):
    hcoord = np.arange(xmin, xmax, w) + rows % 2 * w / 2
    vcoord = [ymax - rows * d * 0.75] * n_cols
    for x, y in zip(hcoord, vcoord):
        hexes = RegularPolygon((x, y), numVertices=6, radius=d / 2, alpha=
            0.2, edgecolor='k')
        verts = hexes.get_path().vertices
        trans = hexes.get_patch_transform()
        points = trans.transform(verts)
        array_of_hexes.append(Polygon(points))
        ax.add_patch(hexes)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
<mask token>
egdes_hex.geometry.to_file(filename='egdes_hex.geojson', driver='GeoJSON')
hex_grid.geometry.to_file(filename='hex_grid.geojson', driver='GeoJSON')
egdes_hex.reset_index(level=0, inplace=True)
hex_grid.reset_index(inplace=True)
<mask token>
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
folium.GeoJson('TIME_EDGES.geojson', style_function=lambda x: style_edges
    ).add_to(my_map)
folium.GeoJson(hex_grid[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_hex, highlight_function=lambda x: {
    'weight': 1, 'color': 'yellow', 'fillOpacity': 0.2}, tooltip=folium.
    features.GeoJsonTooltip(fields=['index'])).add_to(my_map)
folium.GeoJson(egdes_hex[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_egdes_hex, highlight_function=lambda x:
    {'weight': 1, 'color': 'yellow', 'fillOpacity': 0.2}, tooltip=folium.
    features.GeoJsonTooltip(fields=['index'])).add_to(my_map)
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
my_map.save('hex_grid.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
<mask token>
for cells in hex_grid.itertuples(index=True):
    edges = []
    index_elements = []
    index_cells = []
    for links in TIME_EDGES.itertuples(index=True):
        if links.geometry.intersects(cells.geometry) is True:
            print(
                'OK=========================================================================IN'
                )
            index_elements.append(links.Index)
            index_cells.append(cells.Index)
            STREET = links.u, links.v
            edges.append(STREET)
            elements[cells.Index] = edges
        else:
            print('NOT====IN')
<mask token>
gdf_elements.to_file(filename='gdf_elements.geojson', driver='GeoJSON')
gdf_elements.plot()
<mask token>
hex_grid.plot()
hex_gdf_edges.reset_index(inplace=True)
hex_grid.reset_index(inplace=True)
<mask token>
folium.GeoJson('gdf_elements.geojson', style_function=lambda x: style_edges
    ).add_to(my_map)
<mask token>
folium.GeoJson(hex_gdf_edges[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_edge, highlight_function=lambda x: {
    'weight': 1, 'color': 'red', 'fillOpacity': 1}, tooltip=folium.features
    .GeoJsonTooltip(fields=['index'])).add_to(my_map)
my_map.save('elements_and_links.html')
<mask token>
len(OD)
OD.reset_index(level=0, inplace=True)
<mask token>
all_trips_OD.plot()
all_trips_OD.geometry.to_file(filename='all_trips_OD.geojson', driver='GeoJSON'
    )
<mask token>
path.save('OD_shortest_path.html')
if len(path_OD_edges) > 1:
    df_nodes = []
    for i in range(len(path_OD_edges)):
        route = path_OD_edges[i]
        df_nodes.append(route)
    df_nodes = pd.DataFrame(df_nodes)
    df_nodes.columns = ['u', 'v']
    keys = list(df_nodes.columns.values)
    index_TIME_EDGES = TIME_EDGES.set_index(keys).index
    index_df_nodes = df_nodes.set_index(keys).index
    edges_matched_route_OD = TIME_EDGES[index_TIME_EDGES.isin(index_df_nodes)]
    edges_matched_route_OD.plot()
    edges_matched_route_OD.geometry.to_file(filename=
        'edges_matched_route_OD.geojson', driver='GeoJSON')
    init_OD_travel_time = sum(edges_matched_route_OD['travel time (sec)'])
    ave_LAT = 37.53988692816245
    ave_LON = 15.044971594798902
    my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles=
        'cartodbpositron')
    folium.GeoJson('edges_matched_route_OD.geojson', style_function=lambda
        x: style_edges).add_to(my_map)
    my_map.save('edges_matched_route_OD.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
<mask token>
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['cost']) > 0:
        attr['cost'] = float(attr.get('cost'))
for U, V, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost'] = attr.get('cost')
    zipped = zip(list(TIME_EDGES.u), list(TIME_EDGES.v))
    if (U, V) in zipped:
        travel_time = TIME_EDGES[(TIME_EDGES.u == U) & (TIME_EDGES.v == V)][
            'travel time (sec)']
        print(float(travel_time))
        attr['VIASAT_cost'] = float(travel_time)
        grafo.add_edge(U, V, key, attr_dict=attr)
ox.save_graphml(grafo, filename='CATANIA_VIASAT_cost.graphml')
<mask token>
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
<mask token>
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost']) > 0:
        attr['VIASAT_cost'] = float(attr.get('VIASAT_cost'))
    if len(attr['cost']) > 0:
        attr['cost'] = float(attr.get('cost'))
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
<mask token>
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path.save('OD_shortest_path_length.html')
<mask token>
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path_cost.save('OD_shortest_path_cost.html')
<mask token>
edges_shortest_route_cost.geometry.to_file(filename=
    'edges_shortest_route_cost.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_cost.geojson', style_function=lambda x:
    style_edges).add_to(my_map)
my_map.save('shortest_routes_examples.html')
<mask token>
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path_VIASAT_cost.save('OD_VIASAT_shortest_path_cost.html')
<mask token>
edges_shortest_route_VIASAT_cost.geometry.to_file(filename=
    'edges_shortest_route_VIASAT_cost.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_cost.geojson', style_function=
    lambda x: style_VIASAT).add_to(my_map)
my_map.save('shortest_routes_examples.html')
<mask token>
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
for u, v, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost_penalty'] = attr.get('VIASAT_cost')
    if u == U and v == V:
        penalty = 18000
        print(
            'gotta!============================================================================================='
            )
        attr['VIASAT_cost_penalty'] = float(attr['VIASAT_cost']) + penalty
        grafo.add_edge(u, v, key, attr_dict=attr)
ox.save_graphml(grafo, filename='CATANIA_VIASAT_cost_penalty.graphml')
<mask token>
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost_penalty']) > 0:
        attr['VIASAT_cost_penalty'] = float(attr.get('VIASAT_cost_penalty'))
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if u == U and v == V:
        print(
            'gotta!============================================================================================='
            )
        print(attr['VIASAT_cost_penalty'])
        break
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
<mask token>
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path.save('OD_VIASAT_shortest_path_penalty.html')
<mask token>
edges_shortest_route_VIASAT_penalty.geometry.to_file(filename=
    'edges_shortest_route_VIASAT_penalty.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_penalty.geojson',
    style_function=lambda x: style_VIASAT_penalty).add_to(my_map)
my_map.save('shortest_routes_examples.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
<mask token>
if DT < penalty / 3600:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD * DT * (penalty / 3600 - DT / 2)
else:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD * ((penalty / 3600) ** 2 / 2)
<mask token>
for idx in range(len(hex_grid)):
    print(idx)
    hex_grid_idx = hex_grid[hex_grid.index.isin([idx])]
    for cells in hex_grid_idx.itertuples(index=True):
        edges = []
        index_elements = []
        index_cells = []
        for links in TIME_EDGES.itertuples(index=True):
            if links.geometry.intersects(cells.geometry) is True:
                print(
                    'OK=========================================================================IN'
                    )
                index_elements.append(links.Index)
                index_cells.append(cells.Index)
                STREET = links.u, links.v
                edges.append(STREET)
    gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
    if len(gdf_elements) > 0:
        elements[cells.Index] = gdf_elements
<mask token>
np.save('elements.npy', elements)

<mask token>
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.getcwd()
<mask token>
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
place_country = 'Catania, Italy'
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
TIME_EDGES = gpd.read_file('TIME_EDGES.geojson')
RECORDS_CATANIA = gpd.read_file('FREQUENCIES_and_RECORDS_by_EDGES.geojson')
TIME_EDGES = pd.merge(TIME_EDGES, RECORDS_CATANIA, on=['u', 'v',
    'DESTINATION', 'ORIGIN', 'access', 'area', 'bridge', 'distance', 'hour',
    'id', 'junction', 'key', 'lanes', 'length', 'oneway', 'ref', 'service',
    'speed', 'time', 'timedate', 'track_ID', 'tunnel', 'width',
    'length(km)'], how='left')
TIME_EDGES = TIME_EDGES.rename(columns={'geometry_x': 'geometry'})
TIME_EDGES = TIME_EDGES.rename(columns={'color_x': 'color'})
TIME_EDGES.drop(['geometry_y'], axis=1, inplace=True)
TIME_EDGES.drop(['color_y'], axis=1, inplace=True)
TIME_EDGES['hour'] = TIME_EDGES['hour'].ffill()
TIME_EDGES['timedate'] = TIME_EDGES['timedate'].ffill()
TIME_EDGES.dropna(subset=['hour'], inplace=True)
TIME_EDGES.to_file(filename='TIME_EDGES_complete.geojson', driver='GeoJSON')
edges_highways = pd.read_csv('edges_highways.csv')
edges_highways = edges_highways.drop_duplicates(['u', 'v'])
edges_highways = edges_highways[~edges_highways['highway'].isin([
    'residential', 'secondary_link', 'tertiary_link', 'footway',
    'primary_link', 'secondary', 'tertiary', 'other'])]
MAIN_TIME_EDGES = pd.merge(TIME_EDGES, edges_highways, on=['u', 'v'], how=
    'inner')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')


def haversine(coord1, coord2):
    lon1, lat1 = coord1
    lon2, lat2 = coord2
    R = 6371000
    phi_1 = np.radians(lat1)
    phi_2 = np.radians(lat2)
    delta_phi = np.radians(lat2 - lat1)
    delta_lambda = np.radians(lon2 - lon1)
    a = np.sin(delta_phi / 2.0) ** 2 + np.cos(phi_1) * np.cos(phi_2) * np.sin(
        delta_lambda / 2.0) ** 2
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    meters = R * c
    km = meters / 1000.0
    meters = round(meters)
    km = round(km, 3)
    return meters


xmin, ymin, xmax, ymax = TIME_EDGES.total_bounds
EW = haversine((xmin, ymin), (xmax, ymin))
NS = haversine((xmin, ymin), (xmin, ymax))
d = 1800
w = d * np.sin(np.pi / 3)
n_cols = int(EW / w) + 1
n_rows = int(NS / d) + 10
<mask token>
ax = TIME_EDGES.boundary.plot(edgecolor='black', figsize=(20, 60))
w = (xmax - xmin) / n_cols
d = w / np.sin(np.pi / 3)
array_of_hexes = []
for rows in range(0, n_rows):
    hcoord = np.arange(xmin, xmax, w) + rows % 2 * w / 2
    vcoord = [ymax - rows * d * 0.75] * n_cols
    for x, y in zip(hcoord, vcoord):
        hexes = RegularPolygon((x, y), numVertices=6, radius=d / 2, alpha=
            0.2, edgecolor='k')
        verts = hexes.get_path().vertices
        trans = hexes.get_patch_transform()
        points = trans.transform(verts)
        array_of_hexes.append(Polygon(points))
        ax.add_patch(hexes)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
hex_grid = gpd.GeoDataFrame({'geometry': array_of_hexes}, crs={'init':
    'epsg:4326'})
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
egdes_hex = gpd.sjoin(hex_grid, TIME_EDGES, how='inner', op='intersects')
egdes_hex.geometry.to_file(filename='egdes_hex.geojson', driver='GeoJSON')
hex_grid.geometry.to_file(filename='hex_grid.geojson', driver='GeoJSON')
egdes_hex.reset_index(level=0, inplace=True)
hex_grid.reset_index(inplace=True)
style_hex = {'fillColor': '#00000000', 'color': '#00FFFFFF'}
style_edges = {'fillColor': '#00000000', 'color': '#000000'}
style_egdes_hex = {'fillColor': '#00000000', 'color': '#ff0000'}
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
folium.GeoJson('TIME_EDGES.geojson', style_function=lambda x: style_edges
    ).add_to(my_map)
folium.GeoJson(hex_grid[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_hex, highlight_function=lambda x: {
    'weight': 1, 'color': 'yellow', 'fillOpacity': 0.2}, tooltip=folium.
    features.GeoJsonTooltip(fields=['index'])).add_to(my_map)
folium.GeoJson(egdes_hex[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_egdes_hex, highlight_function=lambda x:
    {'weight': 1, 'color': 'yellow', 'fillOpacity': 0.2}, tooltip=folium.
    features.GeoJsonTooltip(fields=['index'])).add_to(my_map)
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
my_map.save('hex_grid.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
hex_grid = gpd.read_file('hex_grid.geojson')
elements = {}
hex_grid = hex_grid[hex_grid.index == 533]
for cells in hex_grid.itertuples(index=True):
    edges = []
    index_elements = []
    index_cells = []
    for links in TIME_EDGES.itertuples(index=True):
        if links.geometry.intersects(cells.geometry) is True:
            print(
                'OK=========================================================================IN'
                )
            index_elements.append(links.Index)
            index_cells.append(cells.Index)
            STREET = links.u, links.v
            edges.append(STREET)
            elements[cells.Index] = edges
        else:
            print('NOT====IN')
gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
gdf_elements.to_file(filename='gdf_elements.geojson', driver='GeoJSON')
gdf_elements.plot()
hex_gdf_edges = hex_grid[hex_grid.index.isin(index_cells)]
hex_grid.plot()
hex_gdf_edges.reset_index(inplace=True)
hex_grid.reset_index(inplace=True)
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
folium.GeoJson('gdf_elements.geojson', style_function=lambda x: style_edges
    ).add_to(my_map)
style_edge = {'fillColor': '#00000000', 'color': '#ff0000'}
folium.GeoJson(hex_gdf_edges[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_edge, highlight_function=lambda x: {
    'weight': 1, 'color': 'red', 'fillOpacity': 1}, tooltip=folium.features
    .GeoJsonTooltip(fields=['index'])).add_to(my_map)
my_map.save('elements_and_links.html')
<mask token>
gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
OD = gdf_elements.drop_duplicates(['ORIGIN', 'DESTINATION'])
len(OD)
OD.reset_index(level=0, inplace=True)
i = 1
O = OD.ORIGIN.iloc[i]
D = OD.DESTINATION.iloc[i]
all_trips_OD = TIME_EDGES[(TIME_EDGES['ORIGIN'] == O) & (TIME_EDGES[
    'DESTINATION'] == D)]
all_trips_OD.plot()
all_trips_OD.geometry.to_file(filename='all_trips_OD.geojson', driver='GeoJSON'
    )
init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length')
path_OD_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path[1:]))
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path.save('OD_shortest_path.html')
if len(path_OD_edges) > 1:
    df_nodes = []
    for i in range(len(path_OD_edges)):
        route = path_OD_edges[i]
        df_nodes.append(route)
    df_nodes = pd.DataFrame(df_nodes)
    df_nodes.columns = ['u', 'v']
    keys = list(df_nodes.columns.values)
    index_TIME_EDGES = TIME_EDGES.set_index(keys).index
    index_df_nodes = df_nodes.set_index(keys).index
    edges_matched_route_OD = TIME_EDGES[index_TIME_EDGES.isin(index_df_nodes)]
    edges_matched_route_OD.plot()
    edges_matched_route_OD.geometry.to_file(filename=
        'edges_matched_route_OD.geojson', driver='GeoJSON')
    init_OD_travel_time = sum(edges_matched_route_OD['travel time (sec)'])
    ave_LAT = 37.53988692816245
    ave_LON = 15.044971594798902
    my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles=
        'cartodbpositron')
    folium.GeoJson('edges_matched_route_OD.geojson', style_function=lambda
        x: style_edges).add_to(my_map)
    my_map.save('edges_matched_route_OD.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['cost']) > 0:
        attr['cost'] = float(attr.get('cost'))
for U, V, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost'] = attr.get('cost')
    zipped = zip(list(TIME_EDGES.u), list(TIME_EDGES.v))
    if (U, V) in zipped:
        travel_time = TIME_EDGES[(TIME_EDGES.u == U) & (TIME_EDGES.v == V)][
            'travel time (sec)']
        print(float(travel_time))
        attr['VIASAT_cost'] = float(travel_time)
        grafo.add_edge(U, V, key, attr_dict=attr)
ox.save_graphml(grafo, filename='CATANIA_VIASAT_cost.graphml')
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=8.5, tiles='cartodbpositron'
    )
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
file_graphml = 'CATANIA_VIASAT_cost.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost']) > 0:
        attr['VIASAT_cost'] = float(attr.get('VIASAT_cost'))
    if len(attr['cost']) > 0:
        attr['cost'] = float(attr.get('cost'))
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length')
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='cost')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path.save('OD_shortest_path_length.html')
init_shortest_OD_path_cost = nx.shortest_path(grafo, O, D, weight='cost')
path_cost = ox.plot_route_folium(grafo, init_shortest_OD_path_cost,
    route_color='green')
path_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path_cost[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='cost')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path_cost.save('OD_shortest_path_cost.html')
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_cost = TIME_EDGES[index_gdf_edges.isin(index_df_nodes)]
style_edge = {'fillColor': '#00000000', 'color': '#000000'}
edges_shortest_route_cost.geometry.to_file(filename=
    'edges_shortest_route_cost.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_cost.geojson', style_function=lambda x:
    style_edges).add_to(my_map)
my_map.save('shortest_routes_examples.html')
init_shortest_OD_path_VIASAT_cost = nx.shortest_path(grafo, O, D, weight=
    'VIASAT_cost')
path_VIASAT_cost = ox.plot_route_folium(grafo,
    init_shortest_OD_path_VIASAT_cost, route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_cost,
    init_shortest_OD_path_VIASAT_cost[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='VIASAT_cost')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path_VIASAT_cost.save('OD_VIASAT_shortest_path_cost.html')
time_OD_NULL = lr / 3600
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_cost = TIME_EDGES[index_gdf_edges.isin(
    index_df_nodes)]
style_VIASAT = {'fillColor': '#00000000', 'color': '#ff0000'}
edges_shortest_route_VIASAT_cost.geometry.to_file(filename=
    'edges_shortest_route_VIASAT_cost.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_cost.geojson', style_function=
    lambda x: style_VIASAT).add_to(my_map)
my_map.save('shortest_routes_examples.html')
U = 2811353264
V = 2811352855
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
for u, v, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost_penalty'] = attr.get('VIASAT_cost')
    if u == U and v == V:
        penalty = 18000
        print(
            'gotta!============================================================================================='
            )
        attr['VIASAT_cost_penalty'] = float(attr['VIASAT_cost']) + penalty
        grafo.add_edge(u, v, key, attr_dict=attr)
ox.save_graphml(grafo, filename='CATANIA_VIASAT_cost_penalty.graphml')
file_graphml = 'CATANIA_VIASAT_cost_penalty.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost_penalty']) > 0:
        attr['VIASAT_cost_penalty'] = float(attr.get('VIASAT_cost_penalty'))
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if u == U and v == V:
        print(
            'gotta!============================================================================================='
            )
        print(attr['VIASAT_cost_penalty'])
        break
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
init_shortest_OD_path_VIASAT_penalty = nx.shortest_path(grafo, O, D, weight
    ='VIASAT_cost_penalty')
path = ox.plot_route_folium(grafo, init_shortest_OD_path_VIASAT_penalty,
    route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_penalty,
    init_shortest_OD_path_VIASAT_penalty[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='VIASAT_cost_penalty')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path.save('OD_VIASAT_shortest_path_penalty.html')
time_OD_penalty = lr / 3600
gdf_nodes, gdf_edges = ox.graph_to_gdfs(grafo)
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']
keys = list(df_nodes.columns.values)
index_gdf_edges = gdf_edges.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_penalty = gdf_edges[index_gdf_edges.isin(
    index_df_nodes)]
style_VIASAT_penalty = {'fillColor': '#00000000', 'color': '#0000ff'}
edges_shortest_route_VIASAT_penalty.geometry.to_file(filename=
    'edges_shortest_route_VIASAT_penalty.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_penalty.geojson',
    style_function=lambda x: style_VIASAT_penalty).add_to(my_map)
my_map.save('shortest_routes_examples.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
df_flux = edges_matched_route_OD
df_flux['timedate'] = pd.to_datetime(df_flux['timedate'])
df_flux = df_flux.resample(rule='15Min', on='timedate').mean()
df_flux['records'] = df_flux['records'].ffill()
df_flux['records'] = df_flux['records'] / 0.25
BBB = pd.DataFrame(df_flux['records'])
TRAVEL_DEMAND_OD = BBB.records.mean()
DT = time_OD_penalty - time_OD_NULL
if DT < penalty / 3600:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD * DT * (penalty / 3600 - DT / 2)
else:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD * ((penalty / 3600) ** 2 / 2)
elements = {}
for idx in range(len(hex_grid)):
    print(idx)
    hex_grid_idx = hex_grid[hex_grid.index.isin([idx])]
    for cells in hex_grid_idx.itertuples(index=True):
        edges = []
        index_elements = []
        index_cells = []
        for links in TIME_EDGES.itertuples(index=True):
            if links.geometry.intersects(cells.geometry) is True:
                print(
                    'OK=========================================================================IN'
                    )
                index_elements.append(links.Index)
                index_cells.append(cells.Index)
                STREET = links.u, links.v
                edges.append(STREET)
    gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
    if len(gdf_elements) > 0:
        elements[cells.Index] = gdf_elements
<mask token>
np.save('elements.npy', elements)

import os
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.getcwd()
import numpy as np
import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame
from shapely.geometry import Point
import folium
import osmnx as ox
import networkx as nx
import math
import momepy
from shapely import geometry
from shapely.geometry import Point, Polygon
import psycopg2
import datetime
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from folium_stuff_FK_map_matching import plot_graph_folium_FK
import glob
from funcs_network_FK import cost_assignment
import statistics
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
place_country = 'Catania, Italy'
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
TIME_EDGES = gpd.read_file('TIME_EDGES.geojson')
RECORDS_CATANIA = gpd.read_file('FREQUENCIES_and_RECORDS_by_EDGES.geojson')
TIME_EDGES = pd.merge(TIME_EDGES, RECORDS_CATANIA, on=['u', 'v',
    'DESTINATION', 'ORIGIN', 'access', 'area', 'bridge', 'distance', 'hour',
    'id', 'junction', 'key', 'lanes', 'length', 'oneway', 'ref', 'service',
    'speed', 'time', 'timedate', 'track_ID', 'tunnel', 'width',
    'length(km)'], how='left')
TIME_EDGES = TIME_EDGES.rename(columns={'geometry_x': 'geometry'})
TIME_EDGES = TIME_EDGES.rename(columns={'color_x': 'color'})
TIME_EDGES.drop(['geometry_y'], axis=1, inplace=True)
TIME_EDGES.drop(['color_y'], axis=1, inplace=True)
TIME_EDGES['hour'] = TIME_EDGES['hour'].ffill()
TIME_EDGES['timedate'] = TIME_EDGES['timedate'].ffill()
TIME_EDGES.dropna(subset=['hour'], inplace=True)
TIME_EDGES.to_file(filename='TIME_EDGES_complete.geojson', driver='GeoJSON')
edges_highways = pd.read_csv('edges_highways.csv')
edges_highways = edges_highways.drop_duplicates(['u', 'v'])
edges_highways = edges_highways[~edges_highways['highway'].isin([
    'residential', 'secondary_link', 'tertiary_link', 'footway',
    'primary_link', 'secondary', 'tertiary', 'other'])]
MAIN_TIME_EDGES = pd.merge(TIME_EDGES, edges_highways, on=['u', 'v'], how=
    'inner')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')


def haversine(coord1, coord2):
    lon1, lat1 = coord1
    lon2, lat2 = coord2
    R = 6371000
    phi_1 = np.radians(lat1)
    phi_2 = np.radians(lat2)
    delta_phi = np.radians(lat2 - lat1)
    delta_lambda = np.radians(lon2 - lon1)
    a = np.sin(delta_phi / 2.0) ** 2 + np.cos(phi_1) * np.cos(phi_2) * np.sin(
        delta_lambda / 2.0) ** 2
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    meters = R * c
    km = meters / 1000.0
    meters = round(meters)
    km = round(km, 3)
    return meters


xmin, ymin, xmax, ymax = TIME_EDGES.total_bounds
EW = haversine((xmin, ymin), (xmax, ymin))
NS = haversine((xmin, ymin), (xmin, ymax))
d = 1800
w = d * np.sin(np.pi / 3)
n_cols = int(EW / w) + 1
n_rows = int(NS / d) + 10
from matplotlib.patches import RegularPolygon
ax = TIME_EDGES.boundary.plot(edgecolor='black', figsize=(20, 60))
w = (xmax - xmin) / n_cols
d = w / np.sin(np.pi / 3)
array_of_hexes = []
for rows in range(0, n_rows):
    hcoord = np.arange(xmin, xmax, w) + rows % 2 * w / 2
    vcoord = [ymax - rows * d * 0.75] * n_cols
    for x, y in zip(hcoord, vcoord):
        hexes = RegularPolygon((x, y), numVertices=6, radius=d / 2, alpha=
            0.2, edgecolor='k')
        verts = hexes.get_path().vertices
        trans = hexes.get_patch_transform()
        points = trans.transform(verts)
        array_of_hexes.append(Polygon(points))
        ax.add_patch(hexes)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
hex_grid = gpd.GeoDataFrame({'geometry': array_of_hexes}, crs={'init':
    'epsg:4326'})
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
egdes_hex = gpd.sjoin(hex_grid, TIME_EDGES, how='inner', op='intersects')
egdes_hex.geometry.to_file(filename='egdes_hex.geojson', driver='GeoJSON')
hex_grid.geometry.to_file(filename='hex_grid.geojson', driver='GeoJSON')
egdes_hex.reset_index(level=0, inplace=True)
hex_grid.reset_index(inplace=True)
style_hex = {'fillColor': '#00000000', 'color': '#00FFFFFF'}
style_edges = {'fillColor': '#00000000', 'color': '#000000'}
style_egdes_hex = {'fillColor': '#00000000', 'color': '#ff0000'}
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
folium.GeoJson('TIME_EDGES.geojson', style_function=lambda x: style_edges
    ).add_to(my_map)
folium.GeoJson(hex_grid[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_hex, highlight_function=lambda x: {
    'weight': 1, 'color': 'yellow', 'fillOpacity': 0.2}, tooltip=folium.
    features.GeoJsonTooltip(fields=['index'])).add_to(my_map)
folium.GeoJson(egdes_hex[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_egdes_hex, highlight_function=lambda x:
    {'weight': 1, 'color': 'yellow', 'fillOpacity': 0.2}, tooltip=folium.
    features.GeoJsonTooltip(fields=['index'])).add_to(my_map)
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
my_map.save('hex_grid.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
hex_grid = gpd.read_file('hex_grid.geojson')
elements = {}
hex_grid = hex_grid[hex_grid.index == 533]
for cells in hex_grid.itertuples(index=True):
    edges = []
    index_elements = []
    index_cells = []
    for links in TIME_EDGES.itertuples(index=True):
        if links.geometry.intersects(cells.geometry) is True:
            print(
                'OK=========================================================================IN'
                )
            index_elements.append(links.Index)
            index_cells.append(cells.Index)
            STREET = links.u, links.v
            edges.append(STREET)
            elements[cells.Index] = edges
        else:
            print('NOT====IN')
gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
gdf_elements.to_file(filename='gdf_elements.geojson', driver='GeoJSON')
gdf_elements.plot()
hex_gdf_edges = hex_grid[hex_grid.index.isin(index_cells)]
hex_grid.plot()
hex_gdf_edges.reset_index(inplace=True)
hex_grid.reset_index(inplace=True)
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
folium.GeoJson('gdf_elements.geojson', style_function=lambda x: style_edges
    ).add_to(my_map)
style_edge = {'fillColor': '#00000000', 'color': '#ff0000'}
folium.GeoJson(hex_gdf_edges[['index', 'geometry']].to_json(), show=True,
    style_function=lambda x: style_edge, highlight_function=lambda x: {
    'weight': 1, 'color': 'red', 'fillOpacity': 1}, tooltip=folium.features
    .GeoJsonTooltip(fields=['index'])).add_to(my_map)
my_map.save('elements_and_links.html')
<mask token>
gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
OD = gdf_elements.drop_duplicates(['ORIGIN', 'DESTINATION'])
len(OD)
OD.reset_index(level=0, inplace=True)
i = 1
O = OD.ORIGIN.iloc[i]
D = OD.DESTINATION.iloc[i]
all_trips_OD = TIME_EDGES[(TIME_EDGES['ORIGIN'] == O) & (TIME_EDGES[
    'DESTINATION'] == D)]
all_trips_OD.plot()
all_trips_OD.geometry.to_file(filename='all_trips_OD.geojson', driver='GeoJSON'
    )
init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length')
path_OD_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path[1:]))
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path.save('OD_shortest_path.html')
if len(path_OD_edges) > 1:
    df_nodes = []
    for i in range(len(path_OD_edges)):
        route = path_OD_edges[i]
        df_nodes.append(route)
    df_nodes = pd.DataFrame(df_nodes)
    df_nodes.columns = ['u', 'v']
    keys = list(df_nodes.columns.values)
    index_TIME_EDGES = TIME_EDGES.set_index(keys).index
    index_df_nodes = df_nodes.set_index(keys).index
    edges_matched_route_OD = TIME_EDGES[index_TIME_EDGES.isin(index_df_nodes)]
    edges_matched_route_OD.plot()
    edges_matched_route_OD.geometry.to_file(filename=
        'edges_matched_route_OD.geojson', driver='GeoJSON')
    init_OD_travel_time = sum(edges_matched_route_OD['travel time (sec)'])
    ave_LAT = 37.53988692816245
    ave_LON = 15.044971594798902
    my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles=
        'cartodbpositron')
    folium.GeoJson('edges_matched_route_OD.geojson', style_function=lambda
        x: style_edges).add_to(my_map)
    my_map.save('edges_matched_route_OD.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['cost']) > 0:
        attr['cost'] = float(attr.get('cost'))
for U, V, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost'] = attr.get('cost')
    zipped = zip(list(TIME_EDGES.u), list(TIME_EDGES.v))
    if (U, V) in zipped:
        travel_time = TIME_EDGES[(TIME_EDGES.u == U) & (TIME_EDGES.v == V)][
            'travel time (sec)']
        print(float(travel_time))
        attr['VIASAT_cost'] = float(travel_time)
        grafo.add_edge(U, V, key, attr_dict=attr)
ox.save_graphml(grafo, filename='CATANIA_VIASAT_cost.graphml')
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=8.5, tiles='cartodbpositron'
    )
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
file_graphml = 'CATANIA_VIASAT_cost.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost']) > 0:
        attr['VIASAT_cost'] = float(attr.get('VIASAT_cost'))
    if len(attr['cost']) > 0:
        attr['cost'] = float(attr.get('cost'))
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length')
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='cost')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path.save('OD_shortest_path_length.html')
init_shortest_OD_path_cost = nx.shortest_path(grafo, O, D, weight='cost')
path_cost = ox.plot_route_folium(grafo, init_shortest_OD_path_cost,
    route_color='green')
path_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path_cost[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='cost')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path_cost.save('OD_shortest_path_cost.html')
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_cost = TIME_EDGES[index_gdf_edges.isin(index_df_nodes)]
style_edge = {'fillColor': '#00000000', 'color': '#000000'}
edges_shortest_route_cost.geometry.to_file(filename=
    'edges_shortest_route_cost.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_cost.geojson', style_function=lambda x:
    style_edges).add_to(my_map)
my_map.save('shortest_routes_examples.html')
init_shortest_OD_path_VIASAT_cost = nx.shortest_path(grafo, O, D, weight=
    'VIASAT_cost')
path_VIASAT_cost = ox.plot_route_folium(grafo,
    init_shortest_OD_path_VIASAT_cost, route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_cost,
    init_shortest_OD_path_VIASAT_cost[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='VIASAT_cost')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path_VIASAT_cost.save('OD_VIASAT_shortest_path_cost.html')
time_OD_NULL = lr / 3600
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_cost = TIME_EDGES[index_gdf_edges.isin(
    index_df_nodes)]
style_VIASAT = {'fillColor': '#00000000', 'color': '#ff0000'}
edges_shortest_route_VIASAT_cost.geometry.to_file(filename=
    'edges_shortest_route_VIASAT_cost.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_cost.geojson', style_function=
    lambda x: style_VIASAT).add_to(my_map)
my_map.save('shortest_routes_examples.html')
U = 2811353264
V = 2811352855
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
for u, v, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost_penalty'] = attr.get('VIASAT_cost')
    if u == U and v == V:
        penalty = 18000
        print(
            'gotta!============================================================================================='
            )
        attr['VIASAT_cost_penalty'] = float(attr['VIASAT_cost']) + penalty
        grafo.add_edge(u, v, key, attr_dict=attr)
ox.save_graphml(grafo, filename='CATANIA_VIASAT_cost_penalty.graphml')
file_graphml = 'CATANIA_VIASAT_cost_penalty.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost_penalty']) > 0:
        attr['VIASAT_cost_penalty'] = float(attr.get('VIASAT_cost_penalty'))
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if u == U and v == V:
        print(
            'gotta!============================================================================================='
            )
        print(attr['VIASAT_cost_penalty'])
        break
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
init_shortest_OD_path_VIASAT_penalty = nx.shortest_path(grafo, O, D, weight
    ='VIASAT_cost_penalty')
path = ox.plot_route_folium(grafo, init_shortest_OD_path_VIASAT_penalty,
    route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_penalty,
    init_shortest_OD_path_VIASAT_penalty[1:]))
lr = nx.shortest_path_length(grafo, O, D, weight='VIASAT_cost_penalty')
lunghezza = []
for l in path_edges:
    lunghezza.append(grafo[l[0]][l[1]][0]['length'])
print('km:{0:.3f} h:{1:.3f} vm:{2:.0f}'.format(sum(lunghezza) / 1000, lr / 
    3600, sum(lunghezza) / 1000 / lr * 3600))
path.save('OD_VIASAT_shortest_path_penalty.html')
time_OD_penalty = lr / 3600
gdf_nodes, gdf_edges = ox.graph_to_gdfs(grafo)
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']
keys = list(df_nodes.columns.values)
index_gdf_edges = gdf_edges.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_penalty = gdf_edges[index_gdf_edges.isin(
    index_df_nodes)]
style_VIASAT_penalty = {'fillColor': '#00000000', 'color': '#0000ff'}
edges_shortest_route_VIASAT_penalty.geometry.to_file(filename=
    'edges_shortest_route_VIASAT_penalty.geojson', driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_penalty.geojson',
    style_function=lambda x: style_VIASAT_penalty).add_to(my_map)
my_map.save('shortest_routes_examples.html')
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
df_flux = edges_matched_route_OD
df_flux['timedate'] = pd.to_datetime(df_flux['timedate'])
df_flux = df_flux.resample(rule='15Min', on='timedate').mean()
df_flux['records'] = df_flux['records'].ffill()
df_flux['records'] = df_flux['records'] / 0.25
BBB = pd.DataFrame(df_flux['records'])
TRAVEL_DEMAND_OD = BBB.records.mean()
DT = time_OD_penalty - time_OD_NULL
if DT < penalty / 3600:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD * DT * (penalty / 3600 - DT / 2)
else:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD * ((penalty / 3600) ** 2 / 2)
elements = {}
for idx in range(len(hex_grid)):
    print(idx)
    hex_grid_idx = hex_grid[hex_grid.index.isin([idx])]
    for cells in hex_grid_idx.itertuples(index=True):
        edges = []
        index_elements = []
        index_cells = []
        for links in TIME_EDGES.itertuples(index=True):
            if links.geometry.intersects(cells.geometry) is True:
                print(
                    'OK=========================================================================IN'
                    )
                index_elements.append(links.Index)
                index_cells.append(cells.Index)
                STREET = links.u, links.v
                edges.append(STREET)
    gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
    if len(gdf_elements) > 0:
        elements[cells.Index] = gdf_elements
import numpy as np
np.save('elements.npy', elements)


import os

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
os.getcwd()

import numpy as np
import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame
from shapely.geometry import Point
import folium
import osmnx as ox
import networkx as nx
import math
import momepy
# from funcs_network_FK import roads_type_folium
from shapely import geometry
from shapely.geometry import Point, Polygon
import psycopg2
# import db_connect
import datetime
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from folium_stuff_FK_map_matching import plot_graph_folium_FK
import glob
from funcs_network_FK import cost_assignment
import statistics


# load grafo
# file_graphml = 'Catania__Italy.graphml'
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
# ox.plot_graph(grafo)
place_country = "Catania, Italy"
# cost_assignment(file_graphml, place_country)

####################################################
#### overlay Hexagonal grid  #######################


os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
# load all EDGES obtained from the map-matching algorithm
TIME_EDGES = gpd.read_file("TIME_EDGES.geojson")
# TIME_EDGES.plot()


# load all EDGES obtained from the map-matching algorithm (RECORDS and FREQUENCIES)
RECORDS_CATANIA = gpd.read_file("FREQUENCIES_and_RECORDS_by_EDGES.geojson")

### merge TIME_EDGES with RECORDS_CATANIA referred to the counts (records about number of veicles)
TIME_EDGES = pd.merge(TIME_EDGES, RECORDS_CATANIA, on=['u', 'v', 'DESTINATION', 'ORIGIN', 'access',
                                                'area', 'bridge', 'distance', 'hour',
                                                'id', 'junction', 'key', 'lanes', 'length', 'oneway', 'ref', 'service',
                                                'speed', 'time',  'timedate', 'track_ID', 'tunnel' , 'width',
                                                'length(km)'
                                                ], how='left')

TIME_EDGES=TIME_EDGES.rename(columns = {'geometry_x':'geometry'})
TIME_EDGES=TIME_EDGES.rename(columns = {'color_x':'color'})
TIME_EDGES.drop(['geometry_y'], axis=1, inplace = True)
TIME_EDGES.drop(['color_y'], axis=1, inplace = True)

## fill nans by values after non-nan values (for 'hour' and 'timedate')
TIME_EDGES['hour'] = TIME_EDGES['hour'].ffill()
TIME_EDGES['timedate'] = TIME_EDGES['timedate'].ffill()
## drop nan values from the column "hour"
TIME_EDGES.dropna(subset = ['hour'], inplace= True)

# save file
TIME_EDGES.to_file(filename='TIME_EDGES_complete.geojson', driver='GeoJSON')

# load the file assigning the type of highway to each (u,v) edge
edges_highways = pd.read_csv("edges_highways.csv")
# remove duplicates
edges_highways = edges_highways.drop_duplicates(['u', 'v'])
edges_highways = edges_highways[~edges_highways['highway'].isin(['residential', 'secondary_link', 'tertiary_link',
                                                     'footway', 'primary_link', 'secondary', 'tertiary', 'other'])]

MAIN_TIME_EDGES = pd.merge(TIME_EDGES, edges_highways, on=['u', 'v'], how='inner')
# MAIN_TIME_EDGES.plot()
# TIME_EDGES.plot()
# TIME_EDGES = MAIN_TIME_EDGES

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')


# define distance between two coordinates
def haversine(coord1, coord2):
    # Coordinates in decimal degrees (e.g. 43.60, -79.49)
    lon1, lat1 = coord1
    lon2, lat2 = coord2
    R = 6371000  # radius of Earth in meters
    phi_1 = np.radians(lat1)
    phi_2 = np.radians(lat2)
    delta_phi = np.radians(lat2 - lat1)
    delta_lambda = np.radians(lon2 - lon1)
    a = np.sin(delta_phi / 2.0) ** 2 + np.cos(phi_1) * np.cos(phi_2) * np.sin(delta_lambda / 2.0) ** 2
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    meters = R * c  # output distance in meters
    km = meters / 1000.0  # output distance in kilometers
    meters = round(meters)
    km = round(km, 3)
    # print(f"Distance: {meters} m")
    # print(f"Distance: {km} km")
    return meters


# define the extension of the EDGE file
xmin, ymin, xmax, ymax = TIME_EDGES.total_bounds  # lat-long of 2 corners
# East-West extent of Toronto = 42193 metres
EW = haversine((xmin, ymin), (xmax, ymin))
# North-South extent of Toronto = 30519 metres
NS = haversine((xmin, ymin), (xmin, ymax))
# diamter of each hexagon in the grid = 900 metres
d = 1800
# horizontal width of hexagon = w = d* sin(60)
w = d * np.sin(np.pi / 3)
# Approximate number of hexagons per row = EW/w
n_cols = int(EW / w) + 1
# Approximate number of hexagons per column = NS/d
n_rows = int(NS / d) + 10

# Make a hexagonal grid to cover the entire area
from matplotlib.patches import RegularPolygon

ax = TIME_EDGES.boundary.plot(edgecolor='black', figsize=(20, 60))
w = (xmax - xmin) / n_cols  # width of hexagon
d = w / np.sin(np.pi / 3)  # diameter of hexagon
array_of_hexes = []
for rows in range(0, n_rows):
    hcoord = np.arange(xmin, xmax, w) + (rows % 2) * w / 2
    vcoord = [ymax - rows * d * 0.75] * n_cols
    for x, y in zip(hcoord, vcoord):  # , colors):
        hexes = RegularPolygon((x, y), numVertices=6, radius=d / 2, alpha=0.2, edgecolor='k')
        verts = hexes.get_path().vertices
        trans = hexes.get_patch_transform()
        points = trans.transform(verts)
        array_of_hexes.append(Polygon(points))
        ax.add_patch(hexes)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])

hex_grid = gpd.GeoDataFrame({'geometry': array_of_hexes}, crs={'init': 'epsg:4326'})
# hex_grid = gpd.GeoDataFrame({'geometry':array_of_hexes},crs={'init':'epsg:3035'})
# hex_grid.plot()


#############################################################################################
# create basemap
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
#############################################################################################

# get hexagons containing the edges
egdes_hex = gpd.sjoin(hex_grid, TIME_EDGES, how='inner', op='intersects')
# egdes_hex.plot()
egdes_hex.geometry.to_file(filename='egdes_hex.geojson', driver='GeoJSON')
hex_grid.geometry.to_file(filename='hex_grid.geojson', driver='GeoJSON')

egdes_hex.reset_index(level=0, inplace=True)
hex_grid.reset_index(inplace=True)

style_hex = {'fillColor': '#00000000', 'color': '#00FFFFFF'}
style_edges = {'fillColor': '#00000000', 'color': '#000000'}
style_egdes_hex = {'fillColor': '#00000000', 'color': '#ff0000'}

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
folium.GeoJson('TIME_EDGES.geojson', style_function=lambda x: style_edges).add_to((my_map))
# folium.GeoJson('hex_grid.geojson',  style_function=lambda x: style_hex).add_to((my_map))
# folium.GeoJson('egdes_hex.geojson',  style_function=lambda x: style_egdes_hex).add_to((my_map))

# add 'u' and 'v' as highligths for each edge (in blue)
folium.GeoJson(
    # data to plot
    hex_grid[['index', 'geometry']].to_json(),
    show=True,
    style_function=lambda x: style_hex,
    highlight_function=lambda x: {'weight': 1,
                                  'color': 'yellow',
                                  'fillOpacity': 0.2
                                  },
    # fields to show
    tooltip=folium.features.GeoJsonTooltip(
        fields=['index']
    ),
).add_to(my_map)

folium.GeoJson(
    # data to plot
    egdes_hex[['index', 'geometry']].to_json(),
    show=True,
    style_function=lambda x: style_egdes_hex,
    highlight_function=lambda x: {'weight': 1,
                                  'color': 'yellow',
                                  'fillOpacity': 0.2
                                  },
    # fields to show
    tooltip=folium.features.GeoJsonTooltip(
        fields=['index']
    ),
).add_to(my_map)

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
my_map.save("hex_grid.html")

##############################################################################################
##############################################################################################
###### ORIGIN - DESTINATION TRIPS ############################################################

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')


###########################################
## get edges inside each hexagon ##########

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
# load hex_grid
hex_grid = gpd.read_file("hex_grid.geojson")

# find all edges (links) intersected by the hexagonal grid
## 'links' are the edges (u, v)
## an 'element' is the ensemble of links within a cell or intersecting a cell

elements = {}

hex_grid = hex_grid[hex_grid.index == 533]  # select only one particular cell (when using 1800 m diameter cell)
# hex_grid = hex_grid[hex_grid.index.isin([2359, 2290])]  # select two particular cells (only with 900m hexagon diameter)
for cells in hex_grid.itertuples(index=True):
    edges = []
    index_elements = []
    index_cells = []
    for links in TIME_EDGES.itertuples(index=True):
        if links.geometry.intersects(cells.geometry) is True:
            print("OK=========================================================================IN")
            index_elements.append(links.Index)
            index_cells.append(cells.Index)
            STREET = links.u, links.v
            edges.append(STREET)
            elements[cells.Index] = edges
        else:
            print("NOT====IN")


## filter gdf_edges based on index_edges (edges in the cell)
gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
# save as GeoJson file
gdf_elements.to_file(filename='gdf_elements.geojson', driver='GeoJSON')
gdf_elements.plot()
# single exagon (cell)
hex_gdf_edges = hex_grid[hex_grid.index.isin(index_cells)]
hex_grid.plot()

## make index as column
hex_gdf_edges.reset_index(inplace=True)
hex_grid.reset_index(inplace=True)

#############################################################################################
# create basemap
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
#############################################################################################

folium.GeoJson('gdf_elements.geojson', style_function=lambda x: style_edges).add_to((my_map))

style_edge = {'fillColor': '#00000000', 'color': '#ff0000'}
# add 'u' and 'v' as highligths for each edge (in blue)
folium.GeoJson(
    # data to plot
    hex_gdf_edges[['index', 'geometry']].to_json(),
    show=True,
    style_function=lambda x: style_edge,
    highlight_function=lambda x: {'weight': 1,
                                  'color': 'red',
                                  'fillOpacity': 1
                                  },
    # fields to show
    tooltip=folium.features.GeoJsonTooltip(
        fields=['index']
    ),
).add_to(my_map)

## save map
my_map.save("elements_and_links.html")


'''
# load grafo of all Catania
file_graphml = 'Catania__Italy.graphml'
grafo = ox.load_graphml(file_graphml)
ox.plot_graph(grafo)
# get edges and nodes
gdf_nodes, gdf_edges = ox.graph_to_gdfs(grafo)
# save edges as GeoJson file
# gdf_edges.to_file(filename='Catania_edges_60km.geojson', driver='GeoJSON')
gdf_edges.to_csv('Catania_edges_60km.csv')
'''

###################################################################################
## calculate the initial travel timeS and shortest paths between the OD pairs #####
###################################################################################


gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
OD = gdf_elements.drop_duplicates(['ORIGIN', 'DESTINATION'])
len(OD)
OD.reset_index(level=0, inplace=True)
# select one OD pair
i = 1
O = OD.ORIGIN.iloc[i]
D = OD.DESTINATION.iloc[i]
# print('ORIGIN:', O, 'DESTINATION:', D, '================================================')
all_trips_OD = TIME_EDGES[(TIME_EDGES['ORIGIN'] == O) & (TIME_EDGES['DESTINATION'] == D)]
all_trips_OD.plot()  # these are all the ONLY elements intersecting the two cells 2359, 2290
all_trips_OD.geometry.to_file(filename='all_trips_OD.geojson', driver='GeoJSON')


# init_OD_travel_time = sum(all_trips_OD['travel time (min)'])   # not good!!
### Initial OD shortest path
init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length') # using length
path_OD_edges = list(zip(init_shortest_OD_path, init_shortest_OD_path[1:]))
# ox.plot_graph_route(grafo, init_shortest_OD_path, route_color='green', fig_height=12, fig_width=12)
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path.save('OD_shortest_path.html')
## isolate edges in the grafo from 'path_OD_edges'
if len(path_OD_edges) > 1:
    # isolate edges in the grafo from 'all_matched_edges'
    df_nodes = []
    for i in range(len(path_OD_edges)):
        # print(path_OD_edges[i])
        route = path_OD_edges[i]
        df_nodes.append(route)

    df_nodes = pd.DataFrame(df_nodes)
    df_nodes.columns = ['u', 'v']

    # filter TIME_EDGES with df_nodes
    keys = list(df_nodes.columns.values)
    index_TIME_EDGES = TIME_EDGES.set_index(keys).index
    index_df_nodes = df_nodes.set_index(keys).index
    edges_matched_route_OD = TIME_EDGES[index_TIME_EDGES.isin(index_df_nodes)]
    edges_matched_route_OD.plot()
    edges_matched_route_OD.geometry.to_file(filename='edges_matched_route_OD.geojson', driver='GeoJSON')

    ########################################################
    ### Initial OD shortest TIME (based on VIASAT data #####
    ########################################################
    init_OD_travel_time = sum(edges_matched_route_OD['travel time (sec)'])


   # create basemap
    ave_LAT = 37.53988692816245
    ave_LON = 15.044971594798902
    my_map = folium.Map([ave_LAT, ave_LON], zoom_start=11, tiles='cartodbpositron')
    #############################################################################################

    folium.GeoJson('edges_matched_route_OD.geojson', style_function=lambda x: style_edges).add_to((my_map))
    ## save map
    my_map.save("edges_matched_route_OD.html")

##################################################################################
##################################################################################

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')

# load grafo
file_graphml = 'Catania__Italy_cost.graphml'
grafo = ox.load_graphml(file_graphml)
# ox.plot_graph(grafo)

# make attr "cost" as float!!! (units of "cost" are SECONDS!)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['cost']) >0:
        attr['cost'] = float(attr.get("cost"))


########################################################################################
## replace "cost" with "travel time" got from VIASAT analysis for  given u,v pair ######
########################################################################################
#### LONG TIME RUNNING #################################################################

for U, V, key, attr in grafo.edges(keys=True, data=True):
    # print(attr["cost"])
    attr['VIASAT_cost'] = attr.get("cost")
    zipped = zip(list(TIME_EDGES.u), list(TIME_EDGES.v)) # we should use TIME_EDGES....
    if (U,V) in zipped:
        travel_time = TIME_EDGES[(TIME_EDGES.u == U) &
                                             (TIME_EDGES.v == V)]['travel time (sec)']  # in seconds
        print(float(travel_time))
        attr['VIASAT_cost'] = float(travel_time)
        grafo.add_edge(U, V, key, attr_dict=attr)
ox.save_graphml(grafo, filename= 'CATANIA_VIASAT_cost.graphml')

###########################################################################################################
###########################################################################################################
# create basemap
ave_LAT = 37.53988692816245
ave_LON = 15.044971594798902
my_map = folium.Map([ave_LAT, ave_LON], zoom_start=8.5, tiles='cartodbpositron')
###########################################################################################################
###########################################################################################################

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')

file_graphml = 'CATANIA_VIASAT_cost.graphml'
grafo = ox.load_graphml(file_graphml)


# make attr "cost" as float!!! (units of "cost" are SECONDS!)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost']) > 0:
        attr['VIASAT_cost'] = float(attr.get("VIASAT_cost"))
    if len(attr['cost']) >0:
        attr['cost'] = float(attr.get("cost"))

# change directory
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')

init_shortest_OD_path = nx.shortest_path(grafo, O, D, weight='length')  # using cost (time)
path = ox.plot_route_folium(grafo, init_shortest_OD_path, route_color='green')
path_edges = list(zip(init_shortest_OD_path,init_shortest_OD_path[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='cost')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path.save('OD_shortest_path_length.html')

## find shortest path based on the "cost" (time)
init_shortest_OD_path_cost = nx.shortest_path(grafo, O, D, weight='cost')  # using cost (time)
path_cost = ox.plot_route_folium(grafo, init_shortest_OD_path_cost, route_color='green')
path_edges = list(zip(init_shortest_OD_path,init_shortest_OD_path_cost[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='cost')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path_cost.save('OD_shortest_path_cost.html')

# add shortest path (by "cost" to the my_map
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']

# filter 'TIME_EDGES' with 'df_nodes'
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_cost = TIME_EDGES[index_gdf_edges.isin(index_df_nodes)]

## add plot in Folium map
# save first as geojson file
style_edge = {'fillColor': '#00000000', 'color': '#000000'}
edges_shortest_route_cost.geometry.to_file(filename='edges_shortest_route_cost.geojson',
                                 driver='GeoJSON')
folium.GeoJson('edges_shortest_route_cost.geojson', style_function=lambda x: style_edges).add_to((my_map))
my_map.save("shortest_routes_examples.html")


# find shortet path by "VIASAT cost" that is the "time" got from the VIASAT data
init_shortest_OD_path_VIASAT_cost = nx.shortest_path(grafo, O, D, weight='VIASAT_cost')  # using cost (time)
path_VIASAT_cost = ox.plot_route_folium(grafo, init_shortest_OD_path_VIASAT_cost, route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_cost,init_shortest_OD_path_VIASAT_cost[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='VIASAT_cost')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path_VIASAT_cost.save('OD_VIASAT_shortest_path_cost.html')
# total time from O--D at normal conditions
time_OD_NULL = lr/3600

# add shortest path (by "VIASAT cost" to the my_map
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']

# filter 'TIME_EDGES' with 'df_nodes'
keys = list(df_nodes.columns.values)
index_gdf_edges = TIME_EDGES.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_cost = TIME_EDGES[index_gdf_edges.isin(index_df_nodes)]

## add plot in Folium map
# save first as geojson file
style_VIASAT = {'fillColor': '#00000000', 'color': '#ff0000'}
edges_shortest_route_VIASAT_cost.geometry.to_file(filename='edges_shortest_route_VIASAT_cost.geojson',
                                 driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_cost.geojson', style_function=lambda x: style_VIASAT).add_to((my_map))
my_map.save("shortest_routes_examples.html")


# get an edge that is on the same route....
U = 2811353264
V = 2811352855


os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')
# close a LINK (u,v pair) by adding sufficiently large penalty M (time in seconds ~ 5 hours = 18000 secs)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    attr['VIASAT_cost_penalty'] = attr.get("VIASAT_cost")
    if u == U and v == V:
        penalty = 18000  # seconds of closure of the link (u,v)
        print("gotta!=============================================================================================")
        # break
        attr['VIASAT_cost_penalty'] = float(attr['VIASAT_cost']) + penalty
        grafo.add_edge(u, v, key, attr_dict=attr)
ox.save_graphml(grafo, filename= 'CATANIA_VIASAT_cost_penalty.graphml')


## load grafo again...
file_graphml = 'CATANIA_VIASAT_cost_penalty.graphml'
grafo = ox.load_graphml(file_graphml)
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if len(attr['VIASAT_cost_penalty']) > 0:
        attr['VIASAT_cost_penalty'] = float(attr.get("VIASAT_cost_penalty"))


# check...
for u, v, key, attr in grafo.edges(keys=True, data=True):
    if u == U and v == V:
        print("gotta!=============================================================================================")
        print(attr['VIASAT_cost_penalty'])
        break

os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing\\vulnerability')
# get shortest path again...
init_shortest_OD_path_VIASAT_penalty = nx.shortest_path(grafo, O, D, weight='VIASAT_cost_penalty')  # using cost (time)
path = ox.plot_route_folium(grafo, init_shortest_OD_path_VIASAT_penalty, route_color='green')
path_edges = list(zip(init_shortest_OD_path_VIASAT_penalty,init_shortest_OD_path_VIASAT_penalty[1:]))
lr = nx.shortest_path_length(grafo, O,D, weight='VIASAT_cost_penalty')
lunghezza=[]
for l in path_edges:
  lunghezza.append(grafo [l[0]] [l[1]] [0]['length'])  # get only the length for each arch between 2 path edges, [0] it the key = 0
print("km:{0:.3f} h:{1:.3f} vm:{2:.0f}".format(sum(lunghezza)/1000, lr/3600, sum(lunghezza)/1000/lr*3600))  # units == km
path.save('OD_VIASAT_shortest_path_penalty.html')
# total time from O--D with penalty
time_OD_penalty = lr/3600


gdf_nodes, gdf_edges = ox.graph_to_gdfs(grafo)
# add shortest path (by "VIASAT cost" to the my_map
df_nodes = pd.DataFrame(path_edges)
df_nodes.columns = ['u', 'v']

# filter 'TIME_EDGES' with 'df_nodes'
keys = list(df_nodes.columns.values)
index_gdf_edges = gdf_edges.set_index(keys).index
index_df_nodes = df_nodes.set_index(keys).index
edges_shortest_route_VIASAT_penalty = gdf_edges[index_gdf_edges.isin(index_df_nodes)]

## add plot in Folium map
# save first as geojson file
style_VIASAT_penalty = {'fillColor': '#00000000', 'color': '#0000ff'}
edges_shortest_route_VIASAT_penalty.geometry.to_file(filename='edges_shortest_route_VIASAT_penalty.geojson',
                                 driver='GeoJSON')
folium.GeoJson('edges_shortest_route_VIASAT_penalty.geojson', style_function=lambda x: style_VIASAT_penalty).add_to((my_map))
my_map.save("shortest_routes_examples.html")


#################################################################################
#################################################################################
#################################################################################
#################################################################################

#############################
### closure IMPACT ##########
#############################

## flux [N vehicles/hour] per OD pair (travel demand per unit time)
os.chdir('C:\\ENEA_CAS_WORK\\Catania_RAFAEL\\postprocessing')

df_flux = edges_matched_route_OD
df_flux['timedate'] = pd.to_datetime(df_flux['timedate'])

# print(df_flux.timedate.astype('<M8[h]'))
# BBB = (df_flux.groupby([df_flux.timedate.astype('<M8[h]')]).mean()).records
# GGG = (df_flux.groupby([df_flux.hour.astype('<M8[h]')]).mean())
# BBB = pd.DataFrame(BBB)
# TRAVEL_DEMAND_OD = BBB.records.mean()

## use 15 minutes average to calculate FLUX [vehicles/hour]
df_flux = df_flux.resample(rule='15Min', on='timedate').mean()
df_flux['records'] = df_flux['records'].ffill()  # 15 minutes flux
df_flux['records'] = df_flux['records']/0.25
BBB = pd.DataFrame(df_flux['records'])
TRAVEL_DEMAND_OD = BBB.records.mean()  # veichles/hour

## difference time between closure and normal conditions
DT = time_OD_penalty - time_OD_NULL  # (hours)

## closure impact
if DT < penalty/3600:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD*DT*( (penalty/3600) - DT/2)  # hours
else:
    DELAY = 0
    DELAY = DELAY + TRAVEL_DEMAND_OD* (((penalty/3600)**2)/2)  # hours



###########################################################
#### save all the elements into a dictionary  #############
###########################################################
elements = {}

# get elements (all edges) within each hexagonal cell)
for idx in range(len(hex_grid)):
    print(idx)
    # hex_grid = hex_grid[hex_grid.index.isin([533])]  # select one particular cells
    hex_grid_idx = hex_grid[hex_grid.index.isin([idx])]  # select one cells by time
    for cells in hex_grid_idx.itertuples(index=True):
        edges = []
        index_elements = []
        index_cells = []
        for links in TIME_EDGES.itertuples(index=True):
            if links.geometry.intersects(cells.geometry) is True:
                print("OK=========================================================================IN")
                index_elements.append(links.Index)
                index_cells.append(cells.Index)
                STREET = links.u, links.v
                edges.append(STREET)
                # elements[cells.Index] = edges
    gdf_elements = TIME_EDGES[TIME_EDGES.index.isin(index_elements)]
    if len(gdf_elements) > 0:
        elements[cells.Index] = gdf_elements
        # gdf_elements.plot()


import numpy as np
# Save
np.save('elements.npy', elements)

