Added the trajectory dataframe that does not work correctly yet

config/RunRimini.sh

@@ -1,4 +1,4 @@
-#./ccm/build.ps1 -UseVCPKG -DisableInteractive -DoNotUpdateTOOL -DoNotDeleteBuildFolder
+./ccm/build.ps1 -UseVCPKG -DisableInteractive -DoNotUpdateTOOL -DoNotDeleteBuildFolder
 ./build_release/city-pro ./bologna_mdt_detailed/2022-05-12/config_bologna.json
 ./build_release/city-pro ./bologna_mdt_detailed/2022-07-01/config_bologna.json
 ./build_release/city-pro ./bologna_mdt_detailed/2022-08-05/config_bologna.json

python/SetRightDirectoriesConfiguration.py

@@ -17,8 +17,8 @@
     config["jump2subnet_analysis"] = True
     config["max_poly_length"] = 6000
     # NOTE: Status -> Try one more class then the usual analysis to see if the Fcm Clustering separates walkers from bikes.
-    config["num_tm"] = config["num_tm"] + 1
-    config["num_tm_subnet"] = config["num_tm_subnet"] + 1
+    config["num_tm"] = config["num_tm"]
+    config["num_tm_subnet"] = config["num_tm_subnet"]
 
     with open(os.path.join(WORKSPACE,"city-pro","bologna_mdt_detailed",StrDate,"config_bologna.json"),"w") as f:
         config = json.dump(config,f,indent=2)

python/work_mdt/script/AnalysisMdt/AnalysisNetworkAllDays.py

@@ -4,7 +4,7 @@
 import numpy as np
 from LatexFunctions import *
 from UsefulStructures import *
-
+import contextily as ctx
 
 
 class NetworkAllDays:
@@ -1092,13 +1092,13 @@ def PlotNumberTrajectoriesGivenClass(self):
         """
         colors = ["blue","red","green","yellow","black","orange","purple","pink","brown","grey"]
         Feature2Label = {"lenght_km":"L","time_hours":"t"}
-        fig,ax = plt.subplots(1,1,figsize = (10,8))
         for Aggregation in ["aggregated"]:
             for Feature in ["lenght_km","time_hours"]:
                 # Compute the alphas that define for a given day fraction of trajectories that are to a class (RhoK)
                 RhoKGivenDay = []
                 RhoxGivenKGivenDay = []
                 RhoxGivenDay = []
+                fig,ax = plt.subplots(1,1,figsize = (10,8))
                 for MobDate in self.ListDailyNetwork:
                     ParametersExpoFitDay = []
                     for StrClass in self.ListStrClassReference:
@@ -1208,28 +1208,89 @@ def PlotComparisonSubnets(self):
             Description:
                 Plot the subnetworks for all days
         """
-        for MobDate in self.ListDailyNetwork:
-            print("Completing the GeoJson for Day: ",MobDate.StrDate)
-            MobDate.CompleteGeoJsonWithClassInfo()
-        NewClasses = np.unique(self.ListDailyNetwork[0].GeoJson["IntClassOrdered_{}".format(self.ListDailyNetwork[0].StrDate)])
-        Colors = ["red","green","yellow","orange","purple","pink","brown","grey"]
+        NewClasses = [0,1,2,3]
+        GeoJson = gpd.read_file("/home/aamad/codice/city-pro/output/bologna_mdt_detailed/BolognaMDTClassInfo.geojson")
+        Colors = ["red","green","purple","orange","yellow","pink","brown","grey"]
         Index2IJ = {0:(0,0),1:(0,1),2:(0,2),3:(1,0),4:(1,1),5:(1,2),6:(2,0),7:(2,1),8:(2,2)}
         for CountClass,NewClass in enumerate(NewClasses):
             fig, ax = plt.subplots(3, 3, figsize=(15, 10))
             for Index,MobDate in enumerate(self.ListDailyNetwork):
                 print("Class {} for Day {}".format(NewClass,MobDate.StrDate))
-                print("Columns: ",MobDate.GeoJson.columns)
-                GdfClass = MobDate.GeoJson.groupby("IntClassOrdered_{}".format(MobDate.StrDate)).get_group(NewClass)
+                try:
+                    GdfClass = MobDate.GeoJson.groupby("StrClassOrdered_{}".format(MobDate.StrDate)).get_group(MobDate.IntClass2StrClass[NewClass])
+                except:
+                    pass
                 i = Index2IJ[Index][0]
                 j = Index2IJ[Index][1]
-                MobDate.GeoJson.plot(ax=ax[i][j], color="black",alpha = 0.5)
-                GdfClass.plot(ax=ax[i][j], color=Colors[CountClass])
+                GeoJson.plot(ax=ax[i][j], color="black",alpha = 0.01, linewidth=0.1)
+                try:
+                    GdfClass.plot(ax=ax[i][j], color=Colors[CountClass], linewidth=2)
+                except:
+                    pass
+                try:
+                    ctx.add_basemap(ax[i][j], crs="EPSG:32632")
+                except:
+                    print("No Basemap")
                 ax[i][j].set_title(MobDate.StrDate)
+                ax[i][j].set_axis_off()
             # Add legend
-            plt.title('SubNetworks for Class {}'.format(NewClass))
+#            plt.title('SubNetworks for Class {}'.format(MobDate.IntClass2StrClass[NewClass]))
             plt.xlabel('Longitude')
             plt.ylabel('Latitude')
+            fig.suptitle('SubNetworks for Class {}'.format(MobDate.IntClass2StrClass[NewClass]), fontsize=16)
+            plt.tight_layout(rect=[0, 0.03, 1, 0.95])
             plt.savefig(os.path.join(self.PlotDir,"SubNetworks_{}.png".format(NewClass)),dpi = 200) 
             plt.close()       
 
-
+    def PlotUnionSubnets(self):
+        """
+            This snippet works just if I consider 4 classes.
+
+        """
+        UniqueClasses = [0,1,2,3]
+        Colors = ["red","green","purple","orange","yellow","pink","brown","grey"]
+        StrUnion = "OrderedUnion_"
+        GeoJson = gpd.read_file("/home/aamad/codice/city-pro/output/bologna_mdt_detailed/BolognaMDTClassInfo.geojson")
+        fig, ax = plt.subplots(2, 2, figsize=(15, 10))
+        Class2Ax = {0:(0,0),1:(0,1),2:(1,0),3:(1,1)}
+        IntClass2StrClass = {0:"1 slowest",1:"2 slowest",2:"middle velocity class",3:"1 quickest"}
+        for Class in UniqueClasses:
+            i = Class2Ax[Class][0]
+            j = Class2Ax[Class][1]
+            filtered_gdf = GeoJson.loc[GeoJson[StrUnion + IntClass2StrClass[Class]] == True].dropna(subset=['geometry'])
+            filtered_gdf.plot(ax=ax[i][j], color=Colors[Class], linewidth=2)
+            try:
+                ctx.add_basemap(ax[i][j], crs="EPSG:32632")
+            except:
+                print("No Basemap")
+            ax[i][j].set_title("Union All Days Class {}".format(Class))
+            ax[i][j].set_axis_off()
+        plt.savefig(os.path.join(self.PlotDir,"UnionSubNetworks.png"),dpi = 200) 
+        plt.close()
+
+
+    def PlotIntersectionSubnets(self):
+        """
+            This snippet works just if I consider 4 classes.
+
+        """
+        UniqueClasses = [0,1,2,3]
+        Colors = ["red","green","purple","orange","yellow","pink","brown","grey"]
+        StrUnion = "OrderedIntersection_"
+        IntClass2StrClass = {0:"1 slowest",1:"2 slowest",2:"middle velocity class",3:"1 quickest"}
+        GeoJson = gpd.read_file("/home/aamad/codice/city-pro/output/bologna_mdt_detailed/BolognaMDTClassInfo.geojson")
+        fig, ax = plt.subplots(2, 2, figsize=(15, 10))
+        Class2Ax = {0:(0,0),1:(0,1),2:(1,0),3:(1,1)}
+        for Class in UniqueClasses:
+            i = Class2Ax[Class][0]
+            j = Class2Ax[Class][1]
+            filtered_gdf = GeoJson.loc[GeoJson[StrUnion + IntClass2StrClass[Class]] == True].dropna(subset=['geometry'])
+            filtered_gdf.plot(ax=ax[i][j], color=Colors[Class], linewidth=2)
+            try:
+                ctx.add_basemap(ax[i][j], crs="EPSG:32632")
+            except:
+                print("No Basemap")
+            ax[i][j].set_title("Union All Days Class {}".format(Class))
+            ax[i][j].set_axis_off()
+        plt.savefig(os.path.join(self.PlotDir,"UnionSubNetworks.png"),dpi = 200) 
+        plt.close()

python/work_mdt/script/AnalysisMdt/AnalysisPaper.py

@@ -202,6 +202,9 @@ def MainComparison(ListNetworkDays,PlotDirAggregated,config,verbose):
     NetAllDays.CompareTimePercorrenceAllDays()
     # Compare sub-nets
     NetAllDays.PlotComparisonSubnets()
+    # Compare sub-nets
+    NetAllDays.PlotIntersectionSubnets()
+    NetAllDays.PlotIntersectionSubnets()
 
 
 if __name__ == "__main__":

python/work_mdt/script/AnalysisMdt/GeoPlot.py

@@ -60,7 +60,7 @@ def BuildListStepsGivenDay(GpdClasses,StrDay,StrMetric):
         Col = StrMetric + StartInterval + "_" + StrDay
         GpdClasses[Col] = np.zeros(len(GpdClasses))
         for Road in GpdClasses["poly_lid"].to_numpy().astype(int):
-            if not np.isnan(Road):
+            if not np.isnan(Road) and Road>0:
                 GpdClasses.at[Road,Col] = Hour2Road2Metric[StartInterval][Road]
             else:
                 print("Time Interval: ",StartInterval)

python/work_mdt/script/AnalysisMdt/analysisPlot.py

@@ -951,3 +951,85 @@ def ScatterFitParams(A, b, classes, types,Class2Type2Colors,Class2Type2Shapes,Ti
     ax.legend(handles = legend,loc='center left', bbox_to_anchor=(1, 0.5))  # Place legend outside the plot
     plt.savefig(os.path.join(PlotDir,SaveName + ".png"))
     plt.close()
+
+
+# UNION ALL DAYS SUBNETWORK
+def PlotIntersection(GpdClasses,UniqueClasses,StrClasses2Color,StrIntersection = "Intersection_"):
+    First = True
+    for Class in UniqueClasses:
+        filtered_gdf = GpdClasses.loc[GpdClasses[StrIntersection + Class] == True].dropna(subset=['geometry'])
+        if len(filtered_gdf) == 0:
+            continue
+        else:        
+            if not First:
+                print("class: ",Class," Color: ",StrClasses2Color[Class]," Column: ",StrIntersection + Class)
+                print("Number of roads to Color: ",len([True for i in GpdClasses[StrIntersection + Class] if i if i == True]))
+
+                filtered_gdf.explore(column = StrIntersection + Class,
+                                    color = StrClasses2Color[Class],
+                                    categories = [True,False],
+                                    legend = True,
+                                    legend_kwds = {'loc': 'upper right'},
+                                    figsize = (8,5),
+                                    tooltip = "poly_lid",
+                                    tooltip_kwds = dict(labels= False),
+                                    name = "Intersection All Days " + Class,
+                                    m = m)
+            else:
+                print("First: Class ",Class," Color: ",StrClasses2Color[Class]," Column: ",StrIntersection + Class)
+                print("Number of roads to Color: ",len([True for i in GpdClasses[StrIntersection + Class] if i if i == True]))
+                m = filtered_gdf.explore(column = "Intersection_" + Class,
+                                    color = StrClasses2Color[Class],
+                                    categories = [True,False],
+                                    legend = True,
+                                    legend_kwds = {'loc': 'upper right'},
+                                    figsize = (8,5),
+                                    tooltip = "poly_lid",
+                                    tooltip_kwds = dict(labels= False),
+                                    name = "Intersection All Days " + Class
+                                    )
+                folium.TileLayer("CartoDB positron", show=False).add_to(m)
+                First = False
+    folium.LayerControl().add_to(m)        
+    return m
+
+
+def PlotUnion(GpdClasses,UniqueClasses,StrClasses2Color,StrUnion = "Union_"):
+    First = True
+    for Class in UniqueClasses:
+        print(StrUnion + Class)
+        filtered_gdf = GpdClasses.loc[GpdClasses[StrUnion + Class] == True].dropna(subset=['geometry'])
+        if len(filtered_gdf) == 0:
+            continue
+        else:
+            if not First:
+                print("class: ",Class," Color: ",StrClasses2Color[Class])
+                print("Number of roads to Color: ",len([True for i in GpdClasses[StrUnion + Class] if i if i == True]))
+                filtered_gdf.explore(column = StrUnion + Class,
+                                    color = StrClasses2Color[Class],
+                                    categories = [True,False],
+                                    legend = True,
+                                    legend_kwds = {'loc': 'upper right'},
+                                    figsize = (8,5),
+                                    tooltip = "poly_lid",
+                                    tooltip_kwds = dict(labels= False),
+                                    name = "Union All Days " + Class,
+                                    m = m)
+            else:
+                print("First: Class ",Class," Color: ",StrClasses2Color[Class])
+                print("Number of roads to Color: ",len([True for i in GpdClasses[StrUnion + Class] if i if i == True]))
+                m = filtered_gdf.explore(column = StrUnion + Class,
+                                    color = StrClasses2Color[Class],
+                                    categories = [True,False],
+                                    legend = True,
+                                    legend_kwds = {'loc': 'upper right'},
+                                    figsize = (8,5),
+                                    tooltip = "poly_lid",
+                                    tooltip_kwds = dict(labels= False),
+                                    name = "Union All Days " + Class
+                                    )
+                folium.TileLayer("CartoDB positron", show=False).add_to(m)
+                First = False
+    folium.LayerControl().add_to(m)        
+    return m
+

source/data_analysis.cpp

@@ -958,6 +958,22 @@ void dump_longest_traj(std::vector<traj_base> &traj)
     outtraj1.close();
 }
 
+/*
+    make_bp_traj: compute the best path trajectories
+    @param traj: vector of trajectories
+    @param config_: configuration parameters
+    @param sigma: standard deviation DEPRECATED
+    @param dataloss: data loss -> counts the number of trajectories that are not considered due to threshold of time
+    @param poly: vector of polies (the roads)
+    @param centers_fcm: vector of fcm centers
+    @param node: vector of nodes
+    @param classes_flux: vector of classes fluxes
+    Description: 
+        compute the best path trajectories
+
+*/
+
+
 void make_bp_traj(std::vector<traj_base> &traj,config &config_,double &sigma,data_loss &dataloss,std::vector<poly_base> &poly,std::vector<centers_fcm_base> &centers_fcm,std::vector<node_base> &node,std::vector<std::map<int,int>> &classes_flux)
 {
     std::cout << "Computing Best Path Trajectories..." << std::endl;
@@ -2668,3 +2684,55 @@ double measure_representativity(const string &label,std::vector<poly_base> &poly
     return sub / all;
 }
 //----------------------------------------------------------------------------------------------------
+
+/*
+*   @brief make_traj_dataframe: Computes the trajectory dataframe:
+*   @param traj vector of trajectories
+*   @details
+*   - Open the file 'config_.cartout_basename + "/" + config_.name_pro + "_traj_dataframe.csv"'
+*   - Write the header "user_id;latitude;longitude;timestamp;speed;road_id"
+*   - For each trajectory in 'traj'
+*       - For each stop_point in the trajectory
+*           - Write the 'id_act' of the trajectory
+*           - Write the 'latitude' of the stop_point
+*           - Write the 'longitude' of the stop_point
+*           - Write the 'itime' of the stop_point
+*           - Write the 'inst_speed' of the stop_point
+*           - Write the 'first' element of the 'path' of the stop_point
+*   - Close the file
+
+*/
+void make_traj_dataframe(std::vector<traj_base> &traj,std::vector<poly_base> &poly){
+    ofstream out_traj(config_.cartout_basename + "/" + config_.name_pro + "_traj_dataframe.csv");
+    ofstream size_path(config_.cartout_basename + "/" + config_.name_pro + "_size_path.csv");
+    out_traj << "user_id;latitude;longitude;timestamp;speed;node_front;node_tail;road_id" << std::endl;
+    for (auto &t : traj)
+    {
+        auto it = t.path.begin();
+        size_path << t.id_act << ";" << t.path.size() << std::endl;
+        size_path.close();
+        for (int sp = 0; sp < t.stop_point.size(); ++sp)
+        {
+//            if (it == t.path.end()) {
+//                break;  // Ensure we don't go past the end of the list
+//            }    
+
+            if (it->first<0){
+                    int index_ = -it->first;
+                    int PolyFront = poly[index_].cid_Fjnct; 
+                    int PolyTail = poly[index_].cid_Tjnct;
+                    int LidPoly = poly[index_].id_local;
+                    out_traj << t.id_act << ";" << t.stop_point[sp].centroid.lat << ";" << t.stop_point[sp].centroid.lon << ";" << t.stop_point[sp].points.front().itime << ";" << t.stop_point[sp].inst_speed << ";" << it->first << ";"<< PolyFront << ";"<< PolyTail << ";"<< LidPoly << std::endl;
+                }
+            else{
+                    int index_ = it->first;
+                    int PolyFront = poly[index_].cid_Fjnct; 
+                    int PolyTail = poly[index_].cid_Tjnct;
+                    int LidPoly = poly[index_].id_local;
+                    out_traj << t.id_act << ";" << t.stop_point[sp].centroid.lat << ";" << t.stop_point[sp].centroid.lon << ";" << t.stop_point[sp].points.front().itime << ";" << t.stop_point[sp].inst_speed << ";" << it->first << ";"<< PolyFront << ";"<< PolyTail << ";"<< LidPoly << std::endl;
+                }
+            ++it;
+        }
+    }
+    out_traj.close();
+}

source/data_analysis.h

@@ -103,4 +103,5 @@ void analysis_subnets(std::vector<traj_base> &traj,std::vector<poly_base> &poly,
 void dump_longest_traj(std::vector<traj_base> &traj);
 std::map<std::string,std::vector<int>> hierarchical_deletion_of_intersection(std::map<std::string, std::vector<int>> subnets80);
 void assign_new_class(std::vector<traj_base> &traj,std::vector<poly_base> &poly,std::map<std::string, std::vector<int>> subnets80);
+void make_traj_dataframe(std::vector<traj_base> &traj,std::vector<poly_base> &poly);
 // ALBI

source/main_city-pro.cpp

@@ -154,6 +154,7 @@ int main(int argc, char **argv) {
 
 //  ACTIVITY -> TRAJECTORY
     make_traj(analysis_container.activity,analysis_container.dataloss,analysis_container.traj,analysis_container.data_notoncarto,analysis_container.presence);//for t in traj_temp 1) activity->traj_temp 2) t.record -> t.stop_point 3)t.stop_point = sp_on_carto (vector clster_base)
+    make_traj_dataframe(analysis_container.traj,analysis_container.poly);
 // FUZZY ALGORITHM -> fcm
     if (config_.enable_multimodality) make_multimodality(analysis_container.traj,analysis_container.config_,analysis_container.centers_fcm);
 // PROJECT TRAJECTORIES on Cartography