upload codes

raw_map/rotterdam_citycenter/roads.cpg

+ISO-8859-1
\ No newline at end of file

raw_map/rotterdam_citycenter/roads.prj

+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
\ No newline at end of file

step1_simplify_graph.py

+import geopandas as gpd
+import pandas as pd
+import numpy as np
+from utils import *
+from shapely.geometry import LineString
+
+try:
+    from config_mapmatch import *
+except ModuleNotFoundError:
+    raise Warning('config file for map matchingnot found')
+ 
+
+def transform_osm():
+    '''
+    Input: OpenStreetMap
+    Output: Each segment of the road is a separate edge with source and target node, with its corresponding features
+    The source and target nodes and their corresponding coordinates are stored in a point_save_path
+    '''
+    osm_df = gpd.read_file(raw_osm_file)
+    puntid_dict, idpunt_dict, _, _ = get_map_dict(osm_df, point_save_path)
+    columns_ = [i for i in osm_df.columns if i!="geometry"]
+    data = {i:[] for i in columns_+["source", "target", "km"]}
+
+    geometry = []
+    for index, row in osm_df.iterrows():
+        points = list(row['geometry'].coords)
+        for i in columns_:
+            data[i] += [row[i] for j in range(len(points)-1)]
+        for j in range(len(points)-1):
+            data["source"].append(puntid_dict[str(points[j])])
+            data["target"].append(puntid_dict[str(points[j+1])])
+            data["km"].append(geo_dist_single(points[j], points[j+1])/1000)
+            geometry.append(LineString([points[j], points[j+1]]))
+    data["geometry"] = geometry
+    df = gpd.GeoDataFrame(data, geometry='geometry')
+    df.crs = 'EPSG:4326'
+
+    df["osm_id"] = range(len(df))
+    df.rename(columns={"maxspeed":"kmh", "osm_id":"edge"}, inplace=True)
+    df.drop(columns=["old_osm_id", "name", "ref"], inplace=True)
+    df = df.set_index('edge')
+    return df
+
+def adjacent_similar_edges(raw_df):
+    '''
+    Input: formulated OpenStreetMap dataframe
+    Output: dictionary of adjacent and similar edges d[start_edge]=next_edge
+    '''
+    ### Get (start point -> edge) and (end point -> edge) dictionary
+    edge_pts_dict={i:[j,k] for i,j,k in zip(raw_df.index, raw_df.source, raw_df.target)}
+    all_pts = np.unique(raw_df[['source', 'target']].values)
+    start_pts_edge, end_pts_edge, connection_dict={i:[] for i in all_pts}, {i:[] for i in all_pts}, {i:[] for i in raw_df.index}
+    for edge in edge_pts_dict.keys():
+        start, end=edge_pts_dict[edge]
+        start_pts_edge[start].append(edge)
+        end_pts_edge[end].append(edge)
+
+    ### We only consider two edges are possible to be merged 
+    ### Satify the following conditions: (1) they are not connected by intersection point (2) same road type and same direction information
+    one_way, street_type = raw_df.oneway.to_numpy(), raw_df['type'].to_numpy()
+    for edge in raw_df.index:
+        end_point=edge_pts_dict[edge][-1]
+        if len(start_pts_edge[end_point])==1 and len(end_pts_edge[end_point])==1:
+            if one_way[edge]==one_way[start_pts_edge[end_point][0]] and street_type[edge]==street_type[start_pts_edge[end_point][0]]:
+                connection_dict[edge].append(start_pts_edge[end_point][0])
+    return connection_dict
+
+def merge_edge_features(merged_edges, raw_df):
+    '''
+    Objective: derive the features for the simplified graph
+    '''
+    coarse2full_edge = {i:j for i, j in enumerate(merged_edges)}
+    full2coarse_edge = {}
+    for key in coarse2full_edge.keys():
+        for item in coarse2full_edge[key]:
+            full2coarse_edge[item]=key
+    full2coarse_edge, coarse2full_edge = full2coarse_edge, coarse2full_edge
+
+    # Get the accumulated distance for each edge, the accumulated distance is the distance for all neighbor edges before
+    accumulated_distance, accumulated_index ={i:0 for i in raw_df.index}, {i:0 for i in raw_df.index}
+    kms = raw_df['km'].to_numpy()
+    for edges in merged_edges:
+        for idx in range(1,len(edges)):
+            accumulated_index[edges[idx]]=idx
+            accumulated_distance[edges[idx]]=kms[edges[idx-1]]+accumulated_distance[edges[idx-1]]
+    raw_df['accu_dist']=raw_df.index.map(accumulated_distance)
+    raw_df['accu_idx']=raw_df.index.map(accumulated_index)
+    raw_df['c_edge']=raw_df.index.map(full2coarse_edge)
+    raw_df["edge"]=raw_df.index
+    raw_df.to_file(complete_graph_path, index=False)
+    raw_df = raw_df.sort_index()
+    raw_df.set_index("edge", inplace=True)
+
+    # Build the simplified graph
+    one_way, street_type = raw_df.oneway.to_numpy(), raw_df['type'].to_numpy()
+    target, source = raw_df['target'].to_numpy(), raw_df['source'].to_numpy()
+    simplified_df=pd.DataFrame({'edge':coarse2full_edge.keys()})
+    simplified_df['source']=simplified_df['edge'].map({i:source[coarse2full_edge[i][0]] for i in coarse2full_edge})
+    simplified_df['target']=simplified_df['edge'].map({i:target[coarse2full_edge[i][-1]] for i in coarse2full_edge})
+    simplified_df['km']=simplified_df['edge'].map({i:accumulated_distance[coarse2full_edge[i][-1]]+kms[coarse2full_edge[i][-1]] for i in coarse2full_edge})
+    simplified_df['oneway']=simplified_df['edge'].map({i:one_way[coarse2full_edge[i][0]] for i in coarse2full_edge})
+    simplified_df['type']=simplified_df['edge'].map({i:street_type[coarse2full_edge[i][0]] for i in coarse2full_edge})
+    return raw_df, simplified_df
+
+def define_cost(simplified_df):
+    # Define the travel cost as the distance in meters
+    simplified_df["cost"]=simplified_df["km"]*1000
+    # Define the cost based on rough classification of availability
+    beta_list = [beta0, beta1, beta2, beta3, beta4]
+    class_list = [class0, class1, class2, class3, class4]
+    for class_, beta_ in zip(class_list, beta_list):
+        simplified_df.loc[(simplified_df["type"].isin(class_)), "cost"] = simplified_df.loc[(simplified_df["type"].isin(class_)), "cost"]*(1+beta_)
+    simplified_df["reverse_cost"]=simplified_df["cost"]*(1+revser_coeffi*simplified_df["oneway"])
+    simplified_df = simplified_df.sort_values('edge')
+    return simplified_df
+
+def simplify_graph():
+    def recursive_programming(parent, edge, d):
+        '''
+        Objective: merge all the adjacent edges with same characteristics into one edge
+        '''
+        if d[edge]==-1 or len(d[edge])==0:
+            return parent, d
+        next_edge=d[edge][0]
+        parent=parent+[edge]
+        if next_edge in parent: # get rid of loop
+            return parent, d
+        parent, d=recursive_programming(parent, next_edge, d)
+        if d[next_edge]!=-1:
+            d[edge]+=d[next_edge]
+        d[next_edge]=-1
+        return parent, d
+    
+    raw_df = transform_osm()
+    d = adjacent_similar_edges(raw_df)
+    for edge in raw_df.index:
+        _, d=recursive_programming([], edge, d)
+    merged_edges = [[key]+d[key] for key in d.keys() if d[key]!=-1]
+    raw_df, simplified_df = merge_edge_features(merged_edges, raw_df)
+    simplified_df = define_cost(simplified_df)
+    simplified_df.to_csv(simplified_graph_path, index=False)
+    simplified_df.set_index("edge", inplace=True)
+    return raw_df, simplified_df
+
+    # Separate a linestring into segments
+    
+    
+
+if __name__ == "__main__":
+    '''
+    This file aims to build a simplified graph from the raw OpenStreetMap data
+    Adjacent road segments with similar characteristics are merged into one edge
+    Adjacent road segments means their connecting nodes only have one incoming and one outgoing road segment
+    '''
+    simplify_graph()
+    print("Graph Simplification Completed, Saved at", simplified_graph_path)

step2_main_road_evaluation.py

+import os
+import time
+
+import glob
+import datetime
+import networkx as nx
+
+import geopandas as gpd
+import pandas as pd
+import numpy as np
+import scipy.sparse as sp
+from scipy.sparse import csr_matrix
+
+from utils import *
+from scipy import stats, sparse
+from geopy.distance import geodesic
+from shapely.geometry import LineString
+
+try:
+    from config_mapmatch import *
+except ModuleNotFoundError:
+    raise Warning('config file for map matchingnot found')
+
+
+
+
+def complete_twoway_edges(edgesDf):
+    '''
+    Objective: for oneway==0 (bidirectional roads), add both direction explicitly in the dataframe
+    '''
+    edgesDf["edge"] = edgesDf.index
+    twoway_ = edgesDf[edgesDf["oneway"]==0]
+    arr1, arr2 = twoway_.source.to_numpy(), twoway_.target.to_numpy()
+    twoway = twoway_.copy()
+    twoway['source'], twoway['target'] = arr2, arr1
+    df = pd.concat([twoway, edgesDf]).reset_index(drop=True)
+    # Add 10% to the distance of footway and residential roads
+    df.loc[df["type"].isin([i for i in bikable_type if i!="cycleway"]), "km"] = df.loc[df["type"].isin(["footway", "residential"]), "km"]*1.1
+    return df
+
+def build_graph(input_df):    
+    '''
+    strict_G: directed complete road network graph
+    bikeable_G: directed road network graph with only bikeable roads
+    full_G: undirected complete road network graph 
+    '''                    
+    strict_G, bikable_G, full_G = nx.DiGraph(), nx.DiGraph(), nx.DiGraph()
+    for idx, row in input_df.iterrows():
+        source = row['source']
+        target = row['target']
+        cost = row['cost']
+        km = row["km"]
+        type_ = row["type"]
+        strict_G.add_edge(source, target, cost=cost, distance=km, roadtype=type_, label=idx)
+        full_G.add_edge(source, target, cost=cost, distance=km, roadtype=type_, label=idx)               
+        if row["oneway"]==1:
+            full_G.add_edge(target, source, cost=cost, distance=km, roadtype=type_, label=idx)
+        if type_ in bikable_type:
+            bikable_G.add_edge(source, target, cost=cost, distance=km, roadtype=type_, label=idx)
+    return strict_G, bikable_G, full_G
+
+def enlarge_downstream_scope(current_node, graph, depth, max_km, raw_km, loop_tolerance_coeff,visited, 
+                    inverse_km, ver_km, 
+                    browsed_km, new_added_km, new_edge_idx,
+                    browsed_km_dict, browsed_trajectory_dict):
+    '''
+    This function records the distance and path for each downstream node of the current node
+    Special case (a) is considered
+    browsed_trajectory_dict: dict[e_4] = [e_1, e_1, e_2, e_3], shows the path from the initial starting node to the current node
+    browsed_km_dict: dict[e_4] = 0.1, shows the distance from the initial starting node to the current node
+    '''
+    if browsed_km>=max_km or depth>=max_depth:
+        return browsed_km_dict, browsed_trajectory_dict
+    
+    # Update the browsed_trajectory_dict
+    browsed_trajectory_dict[current_node]=browsed_trajectory_dict[visited[-1]]+[new_edge_idx]
+
+    if len(visited)>=2: 
+        # Special Case (a), cacluate the vertical and inverse distance between the original start node and the current node
+        original_start_node = str2lst(idpunt_dict[visited[0]])
+        original_end_node = str2lst(idpunt_dict[visited[1]])
+        new_node = str2lst(idpunt_dict[current_node])
+        last_added_node = str2lst(idpunt_dict[visited[-1]])
+        vec1, vec2 = original_start_node-original_end_node, new_node-last_added_node
+        angle = np.arccos(np.sum(vec1*vec2)/np.sqrt(np.sum(vec1*vec1)*np.sum(vec2*vec2)))
+        ver_km += (np.sin(angle)*new_added_km)
+        inverse_km += np.maximum((np.cos(angle)*new_added_km), 0)
+        if (inverse_km<=1e-10 and ver_km>=0.1):
+            pass
+        elif (ver_km>=min(loop_tolerance_coeff*(browsed_km+raw_km), 0.04) or inverse_km>=min(loop_tolerance_coeff*(browsed_km+raw_km), 0.04)):
+            return browsed_km_dict, browsed_trajectory_dict
+        
+    # Find the incoming and outgoing bikeable edges of the current node, as well as all the outgoing unbikeable edges
+    outgoing_edges = [i for i in graph.out_edges(current_node, data=True) if i[1] not in visited and i[2].get('roadtype') in unbikeable_type]
+    out_bikeable = [i for i in graph.out_edges(current_node, data=True) if i[2].get('roadtype') in bikable_type]
+    in_bikeable = [i for i in graph.in_edges(current_node, data=True) if i[2].get('roadtype') in bikable_type]
+    
+    # Update the browsed_km_dict with browsed_km if the current node if connected to the bikeable roads
+    if len(out_bikeable) or len(in_bikeable)>0:
+        browsed_km_dict[current_node]=browsed_km   
+    # If the current node only have no outgoing roads, put the end to it
+    if len(out_bikeable)+len(outgoing_edges)==0:
+        browsed_km_dict[current_node]=-1
+
+    # Recursively search the downstream nodes
+    for next_edge in outgoing_edges:
+        browsed_km_dict, browsed_trajectory_dict = enlarge_downstream_scope(next_edge[1],graph=graph, depth=depth+1, max_km=max_km, raw_km=raw_km, loop_tolerance_coeff=loop_tolerance_coeff,visited=visited+[current_node], 
+                                                                      ver_km=ver_km, inverse_km=inverse_km,
+                                                                      browsed_km=browsed_km+next_edge[2]["distance"], new_added_km=next_edge[2]["distance"],new_edge_idx=next_edge[2]["label"],
+                                                                      browsed_km_dict=browsed_km_dict, browsed_trajectory_dict=browsed_trajectory_dict, 
+                                                )
+    return browsed_km_dict, browsed_trajectory_dict
+
+def enlarge_upstream_scope(current_node, graph, depth, max_km, raw_km, loop_tolerance_coeff,visited, 
+                    inverse_km, ver_km, 
+                    browsed_km, new_added_km, new_edge_idx,
+                    browsed_km_dict, browsed_trajectory_dict):
+
+    if browsed_km>=max_km or depth>=max_depth:
+        return browsed_km_dict, browsed_trajectory_dict
+        
+    browsed_trajectory_dict[current_node]=browsed_trajectory_dict[visited[-1]]+[new_edge_idx]
+    
+    if len(visited)>=2:# and raw_km<0.1:
+        original_start_node = str2lst(idpunt_dict[visited[0]])
+        original_end_node = str2lst(idpunt_dict[visited[1]])
+        new_node = str2lst(idpunt_dict[current_node])
+        last_added_node = str2lst(idpunt_dict[visited[-1]])
+        vec1, vec2 = original_start_node-original_end_node, new_node-last_added_node
+        angle = np.arccos(np.sum(vec1*vec2)/np.sqrt(np.sum(vec1*vec1)*np.sum(vec2*vec2)))            
+        ver_km += (np.sin(angle)*new_added_km)
+        inverse_km += np.maximum((np.cos(angle)*new_added_km), 0)
+        
+        if (inverse_km<=1e-10 and ver_km>=0.1):
+            pass
+        elif (ver_km>=min(loop_tolerance_coeff*(browsed_km+raw_km), 0.04) or inverse_km>=min(loop_tolerance_coeff*(browsed_km+raw_km), 0.04)):# or (ver_km+inverse_km)/(km+raw_km)>=0.3):#if (ver_km/(km+raw_km)>=0.2 or inverse_km/(km+raw_km)>=0.2  or (ver_km+inverse_km)/(km+raw_km)>=0.3):# or ver_km>=0.02 or inverse_km>=0.02):
+            return browsed_km_dict, browsed_trajectory_dict
+    
+    incoming_edges = [i for i in graph.in_edges(current_node, data=True) if i[0] not in visited and i[2].get('roadtype') in unbikeable_type]
+    # print(incoming_edges)
+    out_bikeable = [i for i in graph.out_edges(current_node, data=True) if i[2].get('roadtype') in bikable_type]
+    in_bikeable = [i for i in graph.in_edges(current_node, data=True) if i[2].get('roadtype') in bikable_type]
+    if len(out_bikeable) or len(in_bikeable)>0: 
+        browsed_km_dict[current_node]=browsed_km   
+    if len(incoming_edges)+len(in_bikeable)==0:
+        browsed_km_dict[current_node]=-1
+
+    for next_edge in incoming_edges:
+        browsed_km_dict, browsed_trajectory_dict = enlarge_upstream_scope(next_edge[0], graph=graph, depth=depth+1, max_km=max_km, raw_km=raw_km, loop_tolerance_coeff=loop_tolerance_coeff,visited=visited+[current_node], 
+                                                                      ver_km=ver_km, inverse_km=inverse_km,
+                                                                      browsed_km=browsed_km+next_edge[2]["distance"], new_added_km=next_edge[2]["distance"],  new_edge_idx=next_edge[2]["label"],
+                                                                    browsed_km_dict=browsed_km_dict, browsed_trajectory_dict=browsed_trajectory_dict)
+    return browsed_km_dict, browsed_trajectory_dict
+
+def clean_candidate(bikable_G, full_G, expanded_nodes):
+    '''
+    This function considers special case (b), ensuring the connectivity of the bikeable roads
+    '''
+    # We only need to consider nodes that connects to the bikeable roads
+    candidates = [i for i in expanded_nodes if bikable_G.has_node(i)]
+    bad_expanded_nodes = []
+    for i in candidates:
+        for j in candidates:
+            if i==j:
+                continue
+            # Calculate the shortest path distance for the full road network
+            shortest_path_fullG=nx.shortest_path(full_G, source=i, target=j, weight='distance', method='dijkstra')
+            shortest_dist_fullG = sum(full_G[u][v]["distance"] for u, v in zip(shortest_path_fullG[:-1], shortest_path_fullG[1:]))
+            try:
+                shortest_path_bikeableG=nx.shortest_path(bikable_G, source=i, target=j, weight='distance', method='dijkstra')
+            except:
+                # If the bikeable road network is not connected, we need to remove the corresponding road
+                if j in full_G[i].keys():
+                    bad_expanded_nodes.append(full_G[i][j]["label"])
+                continue
+            shortest_dist_bikeableG = sum(bikable_G[u][v]["distance"] for u, v in zip(shortest_path_bikeableG[:-1], shortest_path_bikeableG[1:]))
+            # If it causes significant detour to the bikeable road network
+            if shortest_dist_bikeableG/shortest_dist_fullG>=network_interrupt_threshold: 
+                bad_expanded_nodes+=[full_G[source_node][target_node]['label'] for source_node, target_node in zip(shortest_path_fullG[:-1], shortest_path_fullG[1:])]
+                break   
+    return bad_expanded_nodes   
+   
+def access_candidate_roads(bikable_G, upstream_browsed_km, downstream_browsed_km, raw_km, upstream_browsed_path, downstream_browsed_path,upstream_bad_expanded_nodes, downstream_bad_expanded_nodes):
+    '''
+    raw_km: the length of our considered unconfirmed car road
+    This function accesses if the found alternative roads are ligible
+    '''
+    all_possibilities = []
+    
+    for upstream_node in upstream_browsed_km.keys():
+        km1 = upstream_browsed_km[upstream_node]
+        for downstream_node in downstream_browsed_km.keys():
+            km2 = downstream_browsed_km[downstream_node]
+
+            if upstream_node==downstream_node:
+                continue
+            try:
+                shortest_path_list=nx.shortest_path(bikable_G, source=upstream_node, target=downstream_node, weight='distance', method='dijkstra')
+            except:
+                # If the bikeable road network is not connected
+                # We could still keep the expanded upstream and downstream nodes
+                # If they are at the frontier of the network besides those listed in frontier_exception_nodes
+                if km1==-1 or km2==-1 and downstream_node not in frontier_exception_nodes:
+                    all_possibilities += upstream_browsed_path[upstream_node]
+                    all_possibilities += downstream_browsed_path[downstream_node]
+                continue
+
+            bike_km = sum([bikable_G[shortest_path_list[i]][shortest_path_list[i + 1]]['distance'] for i in range(0,len(shortest_path_list)-1)])
+            unbikeable_km = raw_km+km1+km2
+
+            if bike_km/unbikeable_km<=case1_theta or (abs(unbikeable_km-bike_km)<=case3_delta and bike_km/unbikeable_km<=case2_theta):
+                for i in upstream_browsed_path[upstream_node]:
+                    if i in upstream_bad_expanded_nodes:
+                        break
+                    all_possibilities.append(i)
+                for i in downstream_browsed_path[downstream_node]:
+                    if i in downstream_bad_expanded_nodes:
+                        break
+                    all_possibilities.append(i)
+    return all_possibilities
+    
+def access_by_row(idx, row, strict_G, bikable_G, full_G):
+    '''
+    This function aims to find the alternative roads for the unconfirmed roads
+    '''
+    source_node, target_node = row["target"], row["source"]
+    loop_tolerance_coeff = ver_invser_km_dict[row["type"]]
+    downstream_browsed_km, downstream_browsed_path = enlarge_downstream_scope(current_node=source_node, graph=strict_G, depth=0, max_km=row["km"]*max_km_coefficient, raw_km=row["km"], loop_tolerance_coeff=loop_tolerance_coeff,visited=[target_node],
+                                                  inverse_km=0, ver_km=0,  
+                                                  browsed_km=0, new_added_km=0, new_edge_idx=idx,
+                                                  browsed_km_dict={}, browsed_trajectory_dict={target_node:[]} )
+    downstream_bad_expanded_nodes = clean_candidate(bikable_G, full_G, downstream_browsed_km.keys())
+    if len(downstream_browsed_km) == 0:
+        return [],[]
+    upstream_browsed_km, upstream_browsed_path = enlarge_upstream_scope(current_node=target_node, graph=strict_G, depth=0, max_km=row["km"]*max_km_coefficient, raw_km=row["km"], loop_tolerance_coeff=loop_tolerance_coeff,visited=[source_node],
+                                                              inverse_km=0, ver_km=0, 
+                                                              browsed_km=0, new_added_km=0, new_edge_idx=idx,
+                                                            browsed_km_dict={}, browsed_trajectory_dict={source_node:[]} )
+    upstream_bad_expanded_nodes = clean_candidate(bikable_G, full_G, upstream_browsed_km.keys())
+    if len(upstream_browsed_km)==0:  
+        return [],[]
+    new_confirmed = access_candidate_roads(bikable_G=bikable_G, upstream_browsed_km=upstream_browsed_km, downstream_browsed_km=downstream_browsed_km, 
+                                                                  raw_km=row["km"], 
+                                                                  upstream_browsed_path=upstream_browsed_path, downstream_browsed_path=downstream_browsed_path,
+                                                                  upstream_bad_expanded_nodes=upstream_bad_expanded_nodes, downstream_bad_expanded_nodes=downstream_bad_expanded_nodes)
+    return new_confirmed
+    
+def evaluate_unconfirmed_main_car_roads():
+    edgesDf = pd.read_csv(simplified_graph_path)
+    edgesDf.set_index("edge", inplace=True)
+    df = complete_twoway_edges(edgesDf)
+    unconfirmed_roads_df = df[df["type"].isin(unconfirmed_roads_type)]
+
+    strict_G, bikable_G, full_G = build_graph(df)
+
+    confirmed_car_roads_idx_lst = []
+    for edge in unconfirmed_roads_df["edge"].unique():
+        bidirection_expanded_roads = []
+        # For each edge, it could have two rows corresponding to two directions
+        # Only when these two directions both have alternative bikeable roads, we can assign it with low availability
+        for idx, row in unconfirmed_roads_df[unconfirmed_roads_df["edge"]==edge].iterrows():
+            candidate = access_by_row(idx, row, strict_G, bikable_G, full_G)
+            if idx in candidate:
+                bidirection_expanded_roads.append(candidate)
+        if len(unconfirmed_roads_df[unconfirmed_roads_df["edge"]==edge])==len(bidirection_expanded_roads):
+            confirmed_car_roads_idx_lst+=[i for j in bidirection_expanded_roads for i in j]
+
+    good_edges = df.drop(confirmed_car_roads_idx_lst)["edge"].unique()
+    confirmaed_car_roads_lst = list(set(df["edge"].unique())-set(good_edges))
+    edgesDf.loc[confirmaed_car_roads_lst, "penality"] = 1
+    edgesDf.loc[edgesDf["type"].isin(confirmed_car_roads_type), "penality"] = 1
+    all_confirmed_car_roads = ",".join([str(i) for i in edgesDf[edgesDf["penality"]==1].index.to_list()])
+
+    os.system("python visual_map_matching.py --function_name=penality_roads_plot --input_file="+all_confirmed_car_roads+" --output_path="+confirmed_roads_html)
+    print("The confirmed car roads are visualized in "+confirmed_roads_html+".\nPlease input the roads that you want to keep, separated with \",\", q to exit.")
+    input_str = input()
+    while input_str != "q":
+        wrongly_selected_ones = input_str.split(",")
+        all_confirmed_car_roads = ",".join([i for i in all_confirmed_car_roads.split(",") if i not in wrongly_selected_ones])
+        os.system("python visual_map_matching.py --function_name=penality_roads_plot --input_file="+all_confirmed_car_roads+" --output_path="+confirmed_roads_html)
+        print("The visualization is updated, do you want to continue? q to exit.")
+        input_str = input()    
+    
+    print("Do you want to manually add confirmed roads?, separated with \",\", q to exit.")
+    input_str = input()
+    while input_str != "q":
+        all_confirmed_car_roads = all_confirmed_car_roads+","+input_str
+        os.system("python visual_map_matching.py --function_name=penality_roads_plot --input_file="+all_confirmed_car_roads+" --output_path="+confirmed_roads_html)
+        print("The visualization is updated, do you want to continue? q to exit.")
+        input_str = input()           
+
+    with open(confirmed_roads_txt, 'w') as file:
+        file.write(all_confirmed_car_roads)
+    
+    all_confirmed_car_roads=[int(i) for i in all_confirmed_car_roads.split(",")]
+    edgesDf["edge"]=edgesDf.index
+    edgesDf["penality"]=0
+    edgesDf.loc[all_confirmed_car_roads, "penality"]=1
+
+    edgesDf["cost"]=edgesDf["km"]*1000
+    beta_list = [beta0, beta1, beta2, beta3]
+    class_list = [class0, class1, class2, class3]
+    for class_, beta_ in zip(class_list, beta_list):
+        edgesDf.loc[(edgesDf["type"].isin(class_))&(edgesDf["penality"]==0), "cost"] = edgesDf.loc[(edgesDf["type"].isin(class_))&(edgesDf["penality"]==0), "cost"]*(1+beta_)
+    edgesDf.loc[edgesDf["penality"]==1, "cost"] = edgesDf.loc[edgesDf["penality"]==1, "cost"]*(1+beta4)
+
+    edgesDf["reverse_cost"]=edgesDf["cost"]*(1+revser_coeffi*edgesDf["oneway"])
+    edgesDf.to_csv(confirmed_coarse_path, index=False)
+    edgesDf.set_index("edge", inplace=True)
+    return edgesDf
+
+if __name__ == "__main__":
+    '''
+    This file aims to evaluate the availability for roads in config_mapmatch/unconfirmed_roads_type
+    Which is based on the existence of nearby alternative roads of config_mapmatch/bikable_type
+    The road network file is stored in config_mapmatch/confirmed_coarse_path
+    The visualtion is available in config_mapmatch/confirmed_roads_html
+    '''
+    _, idpunt_dict, _, _ = get_map_dict(None, point_save_path)
+    evaluate_unconfirmed_main_car_roads()

step3_supercomputer_run.sh

+#!/bin/sh
+#
+#SBATCH --job-name=auto
+#SBATCH --output=auto.log
+#SBATCH --partition=compute
+#SBATCH --time=10:00:00
+#SBATCH --ntasks=59
+#SBATCH --cpus-per-task=1
+#SBATCH --mem-per-cpu=1GB
+#SBATCH --account=research-ceg-tp
+#SBATCH --exclusive
+
+
+# Load modules:
+module load miniconda3
+module load 2022r2
+module load intel/oneapi-all
+export I_MPI_PMI_LIBRARY=/cm/shared/apps/slurm/current/lib64/libpmi2.so
+
+# Set conda env:
+unset CONDA_SHLVL
+source "$(conda info --base)/etc/profile.d/conda.sh"
+
+# Activate conda, run job, deactivate conda
+# Change the name of the conda environment to the one you created
+conda activate network
+
+srun python delftblue_parallel.py DelftBlue_getsPts data/stepII.csv mapMatch_result/ptsDf_/
\ No newline at end of file

step4_supercomputer_run.sh

+#!/bin/sh
+#
+#SBATCH --job-name=auto
+#SBATCH --output=auto.log
+#SBATCH --partition=compute
+#SBATCH --time=10:00:00
+#SBATCH --ntasks=59
+#SBATCH --cpus-per-task=1
+#SBATCH --mem-per-cpu=1GB
+#SBATCH --account=research-ceg-tp
+#SBATCH --exclusive
+
+
+# Load modules:
+module load miniconda3
+module load 2022r2
+module load intel/oneapi-all
+export I_MPI_PMI_LIBRARY=/cm/shared/apps/slurm/current/lib64/libpmi2.so
+
+# Set conda env:
+unset CONDA_SHLVL
+source "$(conda info --base)/etc/profile.d/conda.sh"
+
+# Activate conda, run job, deactivate conda
+# Change the name of the conda environment to the one you created
+conda activate network
+
+srun python delftblue_parallel.py DelftBlue_matchTrace mapMatch_result/rlt_/
\ No newline at end of file

visual_map_tools.py

+import glob
+import datetime
+import argparse
+import random
+import math
+import numpy as np
+import pandas as pd
+import matplotlib.colors as mcolors
+from dateutil.relativedelta import relativedelta
+
+from utils import *
+
+
+def match_one_txtfile(h5_file, txt_folder="mapMatch_result/rlt_/", sample_size=400, group_size=10, selected_tripIDs=None, outputfolder="visualization/map_match/"):
+    nb = h5_file.split("/")[-1].split(".h5")[0]
+    txt_file=txt_folder+str(nb)+".txt"
+    with open(txt_file, 'r') as file:
+        lines = file.readlines()
+        if not os.path.exists(outputfolder+str(nb)+"/"):
+            os.makedirs(outputfolder+str(nb)+"/")
+        if len(glob.glob(outputfolder+str(nb)+"/*.html"))!=math.ceil(sample_size/group_size):
+            print(txt_file+"is visited.")
+            return
+
+        all_df = read_h5(h5_file, time_converse=False)
+        all_df["speed_"]=all_df['speed_'].round(3)
+        all_df["timestamp"] = pd.to_timedelta(all_df.secs, unit='s')+datetime.datetime(2020, 10, 1)
+        all_df.drop(["secs"], axis=1, inplace=True)
+        
+        multipleTraces, multipleRoutes, multipleProjection={}, {}, {}
+        starttrip, rd = int(nb)*ptsDf_group_size, 0
+        line_nb = len(lines)
+        random_integers = random.sample(range(1, line_nb), sample_size)
+
+        # for line in tqdm(file.readlines()[1:], desc="Plot "+txt_file, unit="item"):
+        for i in tqdm(random_integers, desc="Plot "+txt_file, unit="item"):
+            l = lines[i].strip()
+            tripID = int(l.split(",")[0])
+            if selected_tripIDs is not None and tripID not in selected_tripIDs:
+                continue
+            multipleTraces[tripID] = [[], [], [], []] #multipleTraces[tripID][0], obs_lat, obs_color, obs_info
+            multipleRoutes[tripID] = [[], [], None, []] #lons, lats, None, info
+            
+            df = all_df[all_df.tripID==tripID]
+            _, _, projected_lon, projected_lat = get_accurate_start_end_point(df, streetmap, edgesDf)
+            multipleProjection[tripID] = [projected_lon, projected_lat]
+
+            raw_df = all_raw_df[all_raw_df.tripID==tripID]          
+
+            multipleTraces[tripID][1]=multipleTraces[tripID][1]+raw_df.lat.values.tolist()+[None]
+            multipleTraces[tripID][0]=multipleTraces[tripID][0]+raw_df.lon.values.tolist()+[None]
+            for index in range(len(raw_df)):
+                if index in df.index.unique():
+                    row = df.loc[index]
+                    multipleTraces[tripID][3]=multipleTraces[tripID][3]+['tripID: '+str(row["tripID"])+
+                        "<br>Timestamps: "+str(row["timestamp"])+
+                        "<br>Coarse edge: "+str(row["edge"])+
+                        "<br>Uturn: "+str(row["uturn"])+
+                        "<br>Speed: "+str(row["speed_"])+
+                        "<br>Direction: "+str(row["dir"])+
+                        '<br>Fraction: '+str(row["frcalong"])]
+                else:
+                    row = raw_df.iloc[index]
+                    multipleTraces[tripID][3]=multipleTraces[tripID][3]+['tripID: '+str(row["tripID"])+
+                        "<br>Timestamps: "+str(row["timestamp"])]
+            multipleTraces[tripID][3]=multipleTraces[tripID][3]+[' ']
+            multipleTraces[tripID][2]= multipleTraces[tripID][2]+colorFader(raw_df.timestamp, c2='#FDC9C9', c1='#920808')  + ["#000000"]
+
+            edges = [int(i) for i in l.split(",")[1].split(" ")]
+            unique_full_edges = np.concatenate([coarse2full_edge[i] for i in edges])
+            selected_edge=streetmap.loc[unique_full_edges]
+            geo_list = []
+            for index, row in selected_edge.iterrows():
+                for j in row.geometry.coords:
+                    geo_list += [list(j)]
+                multipleRoutes[tripID][3] = multipleRoutes[tripID][3] + ['Road type:'+row["type"]+
+                                '<br>tripID: '+str(tripID)+
+                                '<br>Coarse edge: '+str(row["c_edge"]) for i in range(len(row.geometry.coords))]+[' ']
+                geo_list += [[None, None]]
+            geo_list = np.asarray(geo_list)
+            multipleRoutes[tripID][0], multipleRoutes[tripID][1] = geo_list[:,0], geo_list[:,1]
+            if (rd+1)%group_size==0:
+                plot_tool.plot_trace(outputpath=outputfolder+str(nb)+"/"+str(rd//group_size)+".html", background=background, multipleRoutes=multipleRoutes, multipleTraces=multipleTraces, multipleProjection=multipleProjection)
+                starttrip = tripID+1
+                multipleTraces, multipleRoutes, multipleProjection={}, {}, {}
+            rd += 1
+        if rd%group_size!=0:
+            plot_tool.plot_trace(outputpath=outputfolder+str(nb)+"/"+str(rd//group_size+1)+".html", background=background, multipleRoutes=multipleRoutes, multipleTraces=multipleTraces, multipleProjection=multipleProjection)
+
+
+        #     if (rd+1)%group_size==0:
+        #         plot_tool.plot_trace(outputpath=outputfolder+str(nb)+"/"+str(rd//group_size)+"_"+str(starttrip)+".html", background=background, multipleRoutes=multipleRoutes, multipleTraces=multipleTraces, multipleProjection=multipleProjection)
+        #         starttrip = tripID+1
+        #         multipleTraces, multipleRoutes, multipleProjection={}, {}, {}
+        #     rd += 1
+        # if rd%group_size!=0:
+        #     plot_tool.plot_trace(outputpath=outputfolder+str(nb)+"/"+str(rd//group_size+1)+"_"+str(starttrip)+".html", background=background, multipleRoutes=multipleRoutes, multipleTraces=multipleTraces, multipleProjection=multipleProjection)
+
+def selected_roads_plot(input_folder="mapMatch_result/rlt_/",
+                        outputfolder="visualization/map_match/",
+                        all_tripID=None):
+    # parse the map matching result to get unique coarse roads and full roads
+    edges= []
+    print("Read text files.")
+    edge_day_dict = {i:[] for i in all_tripID}
+    for txt_file in tqdm(glob.glob(input_folder+"*.txt")):
+        if txt_file.split("/")[-1]=="visited.txt":
+            continue
+        h5_file = input_folder+"viterbi/"+txt_file.split("/")[-1].split(".txt")[0]+".h5"        
+        all_df = read_h5(h5_file, time_converse=False)
+        all_df["timestamp"] = pd.to_timedelta(all_df.secs, unit='s')+datetime.datetime(2020, 10, 1)
+        all_df["day"]= (all_df['timestamp'] - pd.Timestamp('2020-10-01')).dt.days
+        ## don't consider the same tripID in different days
+        all_df = all_df.drop_duplicates(subset=["tripID"], keep="first")
+        tripID_time_dict = dict(zip(all_df.tripID, all_df.day))
+        
+        with open(txt_file, 'r') as file:
+            for line in file.readlines()[1:]:
+                l = line.strip()
+                day_ = tripID_time_dict[int(l.split(",")[0])]
+                edges_to_add = [int(i) for i in l.split(",")[1].split(" ")]
+                for edge in edges_to_add:
+                    edge_day_dict[edge].append(day_)
+                edges += edges_to_add
+    edge_nunique_day = {i:len(np.unique(edge_day_dict[i])) for i in edge_day_dict.keys()}
+
+    unique_coarse_edges, counts=np.unique(edges, return_counts=True)
+    counts = [counts[i]/edge_nunique_day[unique_coarse_edges[i]] for i in range(len(unique_coarse_edges))]
+
+    coarse_count = dict(zip(unique_coarse_edges, counts))
+    normalized_counts = (counts - np.min(counts)) / (np.max(counts) - np.min(counts))
+    cmap = plt.get_cmap('cool')
+    rgb_colors = [cmap(norm) for norm in normalized_counts]
+    hex_colors = [mcolors.to_hex(rgb) for rgb in rgb_colors]
+    coarse_color = dict(zip(unique_coarse_edges, hex_colors))
+
+    full_edges = np.concatenate([coarse2full_edge[i] for i in unique_coarse_edges])
+    unique_full_edges = np.unique(full_edges)
+
+    geo_list, info, color = [], [], []
+    selected_edge=streetmap.loc[unique_full_edges]
+
+    print("Plotting heatmap for selected roads.")
+    for index, row in tqdm(selected_edge.iterrows()):
+        lst = [np.asarray(j) for j in row.geometry.coords]
+        lst = add_intermediate_coords(lst)
+        geo_list += lst+[[None, None]]
+        nodes = [row["source"]]+[None for i in range(len(lst)-2)]+[row["target"]]
+        info = info + ['Road type:'+row["type"]+
+                        '<br>Node: '+str(nodes[i])+
+                        '<br>Average Daily Count: '+str(coarse_count[full2coarse_edge[index]])+
+                        '<br>Source: '+str(row['source'])+
+                        '<br>Target: '+str(row['target']) for i in range(len(lst))] + [' ']
+        color += [coarse_color[full2coarse_edge[index]] for i in range(len(lst))] + ["#000000"]
+    geo_list = np.asarray(geo_list)
+    lons, lats = geo_list[:,0], geo_list[:,1]
+    roads = [lons, lats, color, info]
+
+    plot_tool.plot_trace(outputpath=outputfolder+"heatmap_selected_roads.html", marker_size=5, line_marker_size=5, background=background, routes=roads)
+
+    print("Plotting selected roads.")
+    geo_list, info, color = [], [], []
+    for index, row in tqdm(selected_edge.iterrows()):
+        lst = [np.asarray(j) for j in row.geometry.coords]
+        geo_list += lst+[[None, None]]
+        nodes = [row["source"], row["target"]]
+        info = info + ['Road type:'+row["type"]+
+                        '<br>Node: '+str(nodes[i])+
+                        '<br>Count: '+str(coarse_count[full2coarse_edge[index]])+
+                        '<br>Source: '+str(row['source'])+
+                        '<br>Target: '+str(row['target']) for i in range(len(lst))] + [' ']
+    geo_list = np.asarray(geo_list)
+    lons, lats = geo_list[:,0], geo_list[:,1]
+    roads = [lons, lats, None, info]
+    plot_tool.plot_trace(outputpath=outputfolder+"selected_roads.html", marker_size=5, line_marker_size=5, background=background, routes=roads)
+
+    print("Plotting projected location.")
+    lon, lat = [],[]
+    for h5_file in glob.glob(input_folder+"viterbi/*.h5"):
+        df = read_h5(h5_file, time_converse=False)
+        _, _, projected_lon, projected_lat = get_accurate_start_end_point(df, streetmap, edgesDf)
+        lon.append(projected_lon)
+        lat.append(projected_lat)
+    lon, lat=np.concatenate(lon), np.concatenate(lat)
+    plot_tool.heatmap_plot(data=pd.DataFrame({"lon":lon, "lat":lat}), outputpath=outputfolder+"projection_heatmap.html")
+    return 
+
+def get_background(street_map, coarse_street_map,):
+    streetmap = gpd.read_file(street_map)
+    streetmap.set_index("edge", inplace=True)
+    edgesDf = pd.read_csv(coarse_street_map)
+    edgesDf.set_index("edge", inplace=True)
+    coarse2full_edge = {i:[] for i in edgesDf.index}
+    full2coarse_edge = dict(streetmap.c_edge)
+    for full_edge in full2coarse_edge:
+        coarse_edge = full2coarse_edge[full_edge]
+        coarse2full_edge[coarse_edge].append(full_edge)
+
+    geo_list, info = [], []
+    for index, row in streetmap.iterrows():
+        for j in row.geometry.coords:
+            geo_list += [list(j)]
+        nodes = [row["source"], row["target"]]
+        geo_list += [[None, None]]
+        info = info + ['Road type:'+row["type"]+
+                        '<br>Node: '+str(nodes[i])+
+                        '<br>Full edge: '+str(index)+
+                        '<br>Coarse edge: '+str(row["c_edge"])+
+                        '<br>Oneway:'+str(row['oneway'])+
+                        '<br>Source: '+str(row['source'])+
+                        '<br>Target: '+str(row['target'])   for i in range(len(row.geometry.coords))] + [' ']
+    geo_list = np.asarray(geo_list)
+    lons, lats = geo_list[:,0], geo_list[:,1]
+    background = [lons, lats, None, info]
+    return background, coarse2full_edge, full2coarse_edge, streetmap, edgesDf
+
+def penality_roads_plot(input_str, outputpath):
+    unique_coarse_edges = [int(i) for i in input_str.split(",")]
+    full_edges = np.concatenate([coarse2full_edge[i] for i in unique_coarse_edges])
+    unique_full_edges = np.unique(full_edges)
+
+    geo_list, info = [], []
+    selected_edge=streetmap.loc[unique_full_edges]
+
+    for index, row in selected_edge.iterrows():
+        lst = [np.asarray(j) for j in row.geometry.coords]
+        geo_list += lst+[[None, None]]
+        info = info + ['Road type:'+row["type"]+
+                       '<br>Coarse edge: '+str(row["c_edge"])+
+                       '<br>Full edge: '+str(index)+
+                       '<br>Source: '+str(row["source"])+
+                        '<br>Target: '+str(row["target"]) for i in range(len(lst))] + [' ']
+    geo_list = np.asarray(geo_list)
+    lons, lats = geo_list[:,0], geo_list[:,1]
+    roads = [lons, lats, None, info]
+    plot_tool.plot_trace(outputpath=outputpath, background=background, routes=roads)
+
+def count_trajectory(data:pd.DataFrame):
+    trip_dic = {}
+    for df_id, df in data.groupby("tripID"):
+        trip_dic[df_id] = len(df)
+    return np.asarray([min(300,trip_dic[i]) for i in trip_dic.keys()])
+
+def count_time_interval(data:pd.DataFrame):
+    df = data.sort_values(by=['tripID', 'timestamp'])
+    df['time_difference'] = df.groupby('tripID')['timestamp'].diff()
+    arr =  df["time_difference"].to_numpy()/60
+    arr[arr>10]=10
+    return arr
+
+def count_duration(data:pd.DataFrame):
+    trip_dic = {}
+    data = data.sort_values(by=['tripID', 'timestamp'])
+    for df_id, df in data.groupby("tripID"):
+        trip_dic[df_id] = (df["timestamp"].max()-df["timestamp"].min())/60
+    return np.asarray([min(300,trip_dic[i]) for i in trip_dic.keys()])
+
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument("--input_file", type=str)
+    parser.add_argument("--function_name", type=str)
+    parser.add_argument("--ptsDf_group_size", type=int, default=2000)
+    parser.add_argument("--output_path", type=str, default="visualization/map_match/punish_roads.html")
+    parser.add_argument("--outputfolder", type=str, default="visualization/map_match/")
+    parser.add_argument("--raw_file", type=str, default="data/stepII.h5")
+    parser.add_argument("--street_map", type=str, default="mapMatch_result/full_roads.shp")
+    parser.add_argument("--coarse_street_map", type=str, default="mapMatch_result/coarse_roads.csv")
+    
+    args = parser.parse_args()
+    outputfolder=args.outputfolder
+
+    if not os.path.exists(outputfolder):
+        os.mkdir(outputfolder)
+
+    plot_tool = Plot_html()
+    plot_plt = Plot_plt()
+    # get full roads and coarse roads
+    
+    if args.function_name=="data_plot":
+        data = tracetable(args.input_file)
+        traj_duration = count_duration(data=data)
+        plot_plt.hist_density_plot(data=traj_duration, x_label="Duration of routes (min)", y_label="distribution", title=None, bin=100, outputpath=args.outputfolder+"duration.png")
+        traj_count = count_trajectory(data=data)
+        plot_plt.hist_density_plot(data=traj_count, x_label="Nb of data points per route", y_label="distribution", title=None, bin=100, outputpath=args.outputfolder+"nb_datapoint.png")
+        interval_count = count_time_interval(data=data)
+        plot_plt.hist_density_plot(data=interval_count, x_label="Data point interval (min)", y_label="distribution", title=None, bin=100, outputpath=args.outputfolder+"time_interval.png")
+
+    if args.function_name=="shapefile_plot":
+        shapefile = gpd.read_file(args.input_file)
+        geo_list, info = [], []
+        print(shapefile)
+        box = plot_tool.shp_plot_box(shapefile=shapefile, colorby="Layer")
+        # print(box)
+        plot_tool.plot_map_objs(outputpath=args.outputfolder+args.input_file.split("/")[-1][:-4]+".html", line_box=box)
+        exit()
+
+    background, coarse2full_edge, full2coarse_edge, streetmap, edgesDf = get_background(args.street_map, args.coarse_street_map)
+
+    if args.function_name=="penality_roads_plot":
+        penality_roads_plot(input_str=args.input_file, outputpath=args.output_path)
+
+    if args.function_name=="selected_roads_plot":
+        all_raw_df=read_h5(args.raw_file)
+        all_tripID = all_raw_df.tripID.unique()
+        if not os.path.exists(outputfolder):
+            os.mkdir(outputfolder)
+        selected_roads_plot(input_folder=args.input_file, outputfolder=args.outputfolder, all_tripID=all_tripID)
\ No newline at end of file

visualize_tool.py

+import sys
+sys.path.append('../')
+from visual_map_matching import get_background
+from utils import read_h5, get_accurate_start_end_point, colorFader
+import argparse
+import glob
+import os
+import datetime
+from tqdm import tqdm
+import pandas as pd
+import numpy as np
+import plotly.graph_objects as go
+
+def find_trip_in_auto(tripID, inputfolder="../mapMatch_result/rlt_pg/"):
+    for h5_file in glob.glob(inputfolder+"viterbi/*.h5"):
+        all_df = read_h5(h5_file, time_converse=False)
+        if not tripID in all_df['tripID'].values:
+            continue 
+        all_df["speed_"]=all_df['speed_'].round(3)
+        all_df["timestamp"] = pd.to_timedelta(all_df.secs, unit='s')+datetime.datetime(2020, 10, 1)
+        all_df.drop(["secs"], axis=1, inplace=True)   
+        nb = h5_file.split("/")[-1].split(".h5")[0]
+        txt_file=inputfolder+str(nb)+".txt"
+        with open(txt_file, 'r') as file:
+            lines = file.readlines()
+            for line in lines:
+                if str(tripID) == line.split(",")[0]:
+                    break
+        return all_df[all_df['tripID']==tripID], [int(i) for i in line.split(",")[1].split(" ")]
+
+def plot_elements(df, edges, streetmap, edgesDf, tripID, raw_df, coarse2full_edge):
+    Traces, Routes = [[], [], [], []], [[], [], None, []]
+    _, _, projected_lon, projected_lat = get_accurate_start_end_point(df, streetmap, edgesDf)
+    Projection = [projected_lon, projected_lat]
+
+
+    raw_df = raw_df[raw_df['tripID']==tripID]
+    Traces[0], Traces[1] = Traces[0]+raw_df.lon.values.tolist()+[None], Traces[1]+raw_df.lat.values.tolist()+[None]
+
+    for index in range(len(raw_df)):
+        if index in df.index.unique():
+            row = df.loc[index]
+            Traces[3]=Traces[3]+['tripID: '+str(row["tripID"])+
+                "<br>Timestamps: "+str(row["timestamp"])+
+                "<br>Coarse edge: "+str(row["edge"])+
+                "<br>Uturn: "+str(row["uturn"])+
+                "<br>Speed: "+str(row["speed_"])+
+                "<br>Direction: "+str(row["dir"])+
+                '<br>Fraction: '+str(row["frcalong"])]
+        else:
+            row = raw_df.iloc[index]
+            Traces[3]=Traces[3]+['tripID: '+str(row["tripID"])+
+                "<br>Timestamps: "+str(row["timestamp"])]
+    Traces[3]=Traces[3]+[' ']
+    raw_df = raw_df.copy()
+    raw_df.loc[:,"timestamp"] = (raw_df.loc[:,"timestamp"]-datetime.datetime(2020, 10, 1)).dt.total_seconds()
+    Traces[2]= Traces[2]+colorFader(raw_df.timestamp, c2='#FDC9C9', c1='#920808')  + ["#000000"]
+    
+    unique_full_edges = np.concatenate([coarse2full_edge[i] for i in edges])
+    selected_edge=streetmap.loc[unique_full_edges]
+    geo_list = []
+    for index, row in selected_edge.iterrows():
+        for j in row.geometry.coords:
+            geo_list += [list(j)]
+        Routes[3] = Routes[3] + ['Road type:'+row["type"]+
+                        '<br>tripID: '+str(tripID)+
+                        '<br>Coarse edge: '+str(row["c_edge"]) for i in range(len(row.geometry.coords))]+[' ']
+        geo_list += [[None, None]]
+    geo_list = np.asarray(geo_list)
+    Routes[0], Routes[1] = geo_list[:,0], geo_list[:,1]
+    return Traces, Routes, Projection
+
+def plot_html(traces, routes, projection, background, tripID, outputfolder,suffix="", 
+              line_marker_size=20, marker_size=10,):
+    zoom_center = {"lat": 51.9248025, "lon":4.5}
+    fig = go.Figure()
+    # add background
+    lons, lats, color, info = background
+    line_trace0 = go.Scattermapbox(      
+                mode = "markers+lines+text",
+                lon = lons, lat = lats,
+                marker = {'size': line_marker_size, 'color': '#CCFFFF'},line={'width':10, 'color':'#CCFFFF'},
+                name = "Background",
+                text=info)
+    fig.add_trace(line_trace0)
+
+    # add traces
+    lons, lats, color, info = traces
+    line_trace1 = go.Scattermapbox(      
+                mode = "markers+lines+text",
+                lon = lons, lat = lats,
+                marker = {'size': line_marker_size, 'color': color},line={'width':5, 'color':'#FFCCE5'},
+                name = "Trace",
+                text=info)
+    fig.add_trace(line_trace1)
+
+    # add routes
+    lons, lats, color, info = routes
+    if color is None:
+        color, edge_color="#5CA961", "#5CA961"
+    else:  
+        edge_color="#CCFFFF"
+    line_trace2 = go.Scattermapbox(      
+                mode = "markers+lines+text",
+                lon = lons, lat = lats,
+                marker = {'size': line_marker_size, 'color': color},line={'width':10, 'color':edge_color},
+                name = "Selected_edge",
+                text=info)
+    fig.add_trace(line_trace2)
+
+    # add projection
+    lons, lats = projection
+    point_trace = go.Scattermapbox(
+                    mode = "markers",
+                    lon = lons, lat = lats,
+                    marker = {'size': marker_size, 'color': "#FF0000"},
+                    name = str(tripID))
+    fig.add_trace(point_trace) 
+
+    fig.update_layout(mapbox_style="open-street-map")
+    fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0}, mapbox = {
+        'center': zoom_center,
+        'zoom': 13})
+    fig.update_layout(legend={"orientation":"h"})
+    fig.update_layout(height=800, width=2500)
+    fig.for_each_trace(lambda trace: trace.update(visible="legendonly"))
+    fig.write_html(outputfolder+str(tripID)+suffix+".html")  
+    
+if __name__ == "__main__":
+    
+    parser = argparse.ArgumentParser(description='Visualize trips')
+    parser.add_argument('--trip', default="", type=str, help='Trip to visualize')
+    args = parser.parse_args()
+    raw_df = pd.read_csv("selected_stepII.csv")
+    raw_df = read_h5("../data/stepII.h5", True).sort_values(by=["tripID", "timestamp"])
+    raw_df.timestamp = pd.to_datetime(raw_df.timestamp) #, format="%m/%d/%y %H:%M:%S")
+    outputfolder = "trips_visualization/"
+
+    if args.trip == "":
+        all_tripIDs = raw_df.tripID.unique().tolist()
+    else:
+        all_tripIDs = [int(args.trip)]
+
+    os.system("rm -r "+outputfolder)
+    os.mkdir(outputfolder)
+
+    background, coarse2full_edge, full2coarse_edge, streetmap, edgesDf = get_background("../mapMatch_result/full_roads.shp", "../mapMatch_result/coarse_roads.csv")
+    # for tripID in tqdm(all_tripIDs):
+    
+    # for tripID in [int(i) for i in "54259 91200 8409 38839 28415 96376 83565 98376 18947 41151".split(" ")]:
+    for tripID in pd.read_csv("selected_stepII.csv").tripID.unique():
+    # for tripID in [4627, 13205, 16469, 18992, 30041, 31319, 36726, 45671, 74128, 80909, 89819]:
+        mapped_trip_df, mapped_roads = find_trip_in_auto(tripID, "../mapMatch_result/rlt_test/")
+        Traces, Routes, Projection = plot_elements(mapped_trip_df, mapped_roads, streetmap, edgesDf, tripID, raw_df, coarse2full_edge)
+        plot_html(Traces, Routes, Projection, background, tripID, outputfolder, suffix="_auto")
+
+        # mapped_trip_df, mapped_roads = find_trip_in_auto(tripID, "../mapMatch_result/rlt_comp_test/")
+        # Traces, Routes, Projection = plot_elements(mapped_trip_df, mapped_roads, streetmap, edgesDf, tripID, raw_df, coarse2full_edge)
+        # plot_html(Traces, Routes, Projection, background, tripID, outputfolder, suffix="_comp")
+
+        # mapped_trip_df, mapped_roads = find_trip_in_auto(tripID, "../mapMatch_result/rlt_comp_test/")
+        # Traces, Routes, Projection = plot_elements(mapped_trip_df, mapped_roads, streetmap, edgesDf, tripID, raw_df, coarse2full_edge)
+        # plot_html(Traces, Routes, Projection, background, tripID, outputfolder, suffix="_comp")
+
+
+    
+
+