diff --git a/Analysis/plot_traj.py b/Analysis/plot_traj.py index cd6f8154ead0052033b0a8671d4a780b92ddb7b0..e7c1720c02658e55ee2f3f1f881ed651139423ba 100755 --- a/Analysis/plot_traj.py +++ b/Analysis/plot_traj.py @@ -24,35 +24,41 @@ def get_coord(n, k1, ishape_last, sr): return sr.shape(ishape).points[0] if __name__ == "__main__": + import itertools import shapefile import numpy as np from os import path - import csv + import report_graph import sys import matplotlib.pyplot as plt - + import networkx as nx + shp_tr_dir = sys.argv[1] + G = report_graph.read_eddy_graph() ishape_last = np.loadtxt(path.join(shp_tr_dir, "ishape_last.txt"), dtype = int) + pos = {} + + with shapefile.Reader(path.join(shp_tr_dir, "extremum")) as sr: + k1 = sr.record(0)["date_index"] + for n in G: pos[n] = get_coord(n, k1, ishape_last, sr) + + color_iter = itertools.cycle(('#1f77b4', '#aec7e8', '#ff7f0e', '#ffbb78', + '#2ca02c', '#98df8a', '#d62728', '#ff9896', + '#9467bd', '#c5b0d5', '#8c564b', '#c49c94', + '#e377c2', '#f7b6d2', '#7f7f7f', '#c7c7c7', + '#bcbd22', '#dbdb8d', '#17becf', '#9edae5')) plt.figure() - with shapefile.Reader(path.join(shp_tr_dir, "extremum")) as sr, \ - open("edgelist.csv", newline ="") as edgelist: - k1 = sr.record(0)["date_index"] - cr = csv.reader(edgelist, delimiter = " ", skipinitialspace = True) - - # Title lines: - next(cr) - next(cr) - - for row_in in cr: - k = int(row_in[0]) - i = int(row_in[1]) - extr_1 = get_coord((k, i), k1, ishape_last, sr) - k = int(row_in[2]) - i = int(row_in[3]) - extr_2 = get_coord((k, i), k1, ishape_last, sr) - plt.plot((extr_1[0], extr_2[0]), (extr_1[1], extr_2[1]), "-+", - color = "black") + for component, color in zip(nx.weakly_connected_components(G), color_iter): + nx.draw_networkx(G, pos, edgelist = G.edges(component), + nodelist = component, edge_color=color, node_size=10, + node_color=color, with_labels=False) + H= G.subgraph(component) + t = nx.topological_sort(H) + ancestor = next(t) + plt.annotate(str(ancestor), pos[ancestor], color = color, + xytext = (2, 2), textcoords = 'offset points') + plt.tick_params(reset = True) plt.show()