
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)


# load TIME_EDGES file containing times and speed
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()

TIME_EDGES.drop(['level_0'], axis=1, inplace = True)

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

TIME_EDGES.drop(['level_0'], axis=1, inplace = True)
RECORDS_CATANIA.drop(['level_0'], axis=1, inplace = True)

### 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', 'distance', 'hour',
                                                'id', 'junction', 'key',  'length', 'oneway', '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.drop(['index_y'], axis=1, inplace = True)
TIME_EDGES.drop(['index_x'], 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'])

## remove all these fields
# edges_highways = edges_highways[~edges_highways['highway'].isin(['residential', 'secondary_link', 'tertiary_link',
#                                                            'primary_link', 'secondary', 'tertiary'])]
# remove the field 'residential'
# edges_highways = edges_highways[~edges_highways['highway'].isin(['residential', 'pedestrian', 'tertiary', 'tertiary_link',
#                                                                  'secondary', 'secondary_link', 'unclassified',
#                                                                  ['residential', 'unclassified'], ['pedestrian', 'residential'] ])]
## get only.....
edges_highways = edges_highways[edges_highways['highway'].isin(['motorway', 'motorway_link', 'primary', 'trunk',
                                                                'trunk_link', 'road', 'primary_link',
                                                                'secondary', 'secondary_link'])]


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

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


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')


# 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('TIME_EDGES_filtered.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")



###########################################################
#### 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_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)









