
import os

os.chdir('D:\\ENEA_CAS_WORK\\Catania_RAFAEL')
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


### FUNCTION to 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
A

# load TIME_EDGES file containing times and speed
os.chdir('D:\\ENEA_CAS_WORK\\Catania_RAFAEL\\viasat_data')

list_MONTHS = ['FEBRUARY', 'AUGUST']
list_DAYS = ['MONDAY', 'TUESDAY', 'WEDNESDAY',
             'THURSDAY', 'FRIDAY', 'SATURDAY', 'SUNDAY',
             'MORNING_PEAK', 'EVENING_PEAK']

# MORNING_PEAK_FEBRUARY_2019_sottorete_speed_counts_all_vehicles.geojson
# TIME_EDGES = gpd.read_file(path + DAY.strip() + "_" + MONTH.strip() + "_2019_sottorete_speed_counts_all_vehicles.geojson")

### load all EDGES obtained from the map-matching algorithm by DAY and by MONTH
for idxm, month in enumerate(list_MONTHS):
    print(month)
    MONTH = str(month)
    for idxd, day in enumerate(list_DAYS):
        DAY = str(day)
        print(DAY)
        path = 'D:/ENEA_CAS_WORK/Catania_RAFAEL/viasat_data/vulnerability_by_days/'
        TIME_EDGES = gpd.read_file(path + DAY.strip() + "_" + MONTH.strip() + "_2019_sottorete_speed_counts_all_vehicles.geojson")
        # TIME_EDGES.plot()

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

        ### make attr "cost" as float!!! (units of "cost" are SECONDS!)
        for u, v, key, attr in grafo.edges(keys=True, data=True):
            print(attr)
            print(attr.get("cost"))
            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))
            if (U,V) in zipped:
                travel_time = TIME_EDGES[(TIME_EDGES.u == U) &
                                                     (TIME_EDGES.v == V)]['travel_time']  # 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_2019_' + DAY.strip() +
                                         '_' + MONTH.strip() + '.graphml')

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

        ### change directory
        path = 'D:/ENEA_CAS_WORK/Catania_RAFAEL/viasat_data/vulnerability_by_days/'

        ### make an HEXAGONAL GRID ---------- ############################
        ### 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 = 400  #900
        # 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])

        ### HEXAGONAL GRID ----------------- ################################################
        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 around Catania
        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()
        ## SAVE HEXAGONAL GRID containing the EDGES
        egdes_hex.geometry.to_file(filename= path + DAY.strip() + "_" + MONTH.strip() + '_egdes_hex_400m.geojson', driver='GeoJSON')
        ## SAVE HEXAGONAL GRID only
        hex_grid.geometry.to_file(filename= path + DAY.strip() + "_" + MONTH.strip() + '_hex_grid_400m.geojson', driver='GeoJSON')

        #### make the MAP
        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'}

        ## add sottorete to my_map
        folium.GeoJson(path + DAY.strip() + "_" + MONTH.strip() + '_2019_sottorete_speed_counts_all_vehicles.geojson', style_function=lambda x: style_edges).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)


        ### ad Hexagonal GRID
        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)

        ### save the .html map
        my_map.save(path + DAY.strip() + "_" + MONTH.strip() + '_hex_grid_400m.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(path + DAY.strip() + "_" + MONTH.strip() + '_elements_400m.npy', elements)

