Traversing and Analyzing Coronal Hole Connectivity Graph¶

In [1]:

import pickle
import networkx as nx
import os 
import copy
import numpy as np 
import scipy
from scipy import stats 
import matplotlib.pyplot as plt
import matplotlib.colors as c
import matplotlib
import cv2
import datetime as dt

plt.rcParams['savefig.facecolor']='white'
matplotlib.rcParams.update({'font.size': 12})

In [2]:

os.chdir(os.path.dirname('/Users/osissan/PycharmProjects/CHMAP/chmap'))

from chmap.coronal_holes.tracking.tools.plots import plot_coronal_hole
import chmap.database.db_classes as db_class
import chmap.database.db_funs as db_funs
import chmap.maps.magnetic.flux.br_flux as br_flux

Set directory paths.¶

In [3]:

# folder with saved pkl files from running the tracking algorithm. 
res_dir = '/Users/osissan/desktop/CHT_RESULTS/2007to2020'
# folder with downloaded synchronic maps and magnetic flux maps. 
map_dir = '/Users/osissan/desktop/CH_DB'

Log in to the Database.¶

In [4]:

user = "opalissan"
password = ""     # encrypted.   
db_session = db_funs.init_db_conn_old(db_name='mysql-Q', chd_base=db_class.Base,
                                          user=user, password=password)

Attempting to connect to DB server  mysql://opalissan:****pwd****@q.predsci.com:3306/chd
Connection successful

Read in Coronal Hole Connectivity Graph¶

Saved as a pickle file in DropBox.

In [5]:

pickle_file = os.path.join(res_dir + '/connectivity_graph_2007-06-03-10-01-24.pkl')

In [6]:

graph = pickle.load(open(pickle_file, "rb"))

G = graph.G

Analyze Properties of the Connectivity Graph¶

In [7]:

print("Total number of nodes: ", G.number_of_nodes())
print("Total number of edges: ", G.number_of_edges())
print("Total number of subplots: ", len(list(nx.connected_components(G))))

Total number of nodes:  4127
Total number of edges:  4146
Total number of subplots:  81

In [8]:

# order subgraphs based on average node area
subgraph_ordered = graph.order_subgraphs_based_on_area()

In [9]:

# save the total number of nodes in a subgraph
num_nodes_list = []
for ii, g in enumerate(subgraph_ordered):
    subgraph = G.subgraph(g)
    num_nodes_list.append(subgraph.number_of_nodes())

In [10]:

fig, ax = plt.subplots(figsize=(20, 6))
_ = ax.scatter(np.arange(len(num_nodes_list)), num_nodes_list)
_ = ax.set_ylabel("Number of Nodes")
_ = ax.set_xlabel("Subgraph Index")
_ = ax.set_title("Connected Subgraphs Total Number of Nodes \n Ordered by Average Node Area")

plt.savefig(res_dir + '/figures/subgraph_node_size.png')

Note¶

From the plot above, it is visible that the graph has three main subgraphs. Assumption: One of the main subgraphs is connected to the north pole and the other to the south. Lets examine the first and second subgraph (with the largest number of nodes $\approx$ 1750)

In [11]:

main_subgraph_1 = G.subgraph(subgraph_ordered[0])
main_subgraph_2 = G.subgraph(subgraph_ordered[1])

In [12]:

# list with all weights in G. 
w_list = [w["weight"] for u,v,w in G.edges.data()]

# analyze the weight between nodes that have the same id and nodes that have a different id. 
adj_diff_list = []
adj_same_list = []

# loop over all edges in G. 
for u,v,w in G.edges.data():
    # different id. 
    if G.nodes[u]["id"] != G.nodes[v]["id"]:
        adj_diff_list.append(w["weight"])
    # same id. 
    else:
        adj_same_list.append(w["weight"])

In [13]:

print("Number of edges between nodes of the same ID = ", len(adj_same_list))
print("Number of edges between nodes with different ID = ", len(adj_diff_list))

Number of edges between nodes of the same ID =  4012
Number of edges between nodes with different ID =  134

In [14]:

fig, ax = plt.subplots(ncols=3, nrows=2, sharex=True, figsize=(12, 6))
mu, sigma = scipy.stats.norm.fit(w_list)
n , bins, _ = ax[0][0].hist(w_list, 25, density=True)
best_fit_line = stats.norm.pdf(bins, mu, sigma)
_ = ax[0][0].plot(bins, best_fit_line, color="red")
_ = ax[0][0].set_title("All Edge Weights \n $\mu$ = %.2f, median = %.2f, $\sigma$ = %.2f" % (mu, np.median(w_list), sigma))
_ = ax[0][0].set_ylabel("Probability density")


_ = ax[1][0].hist(w_list, 25)
_ = ax[1][0].set_xlabel("Weight")
_ = ax[1][0].set_ylabel("Count")


mu2, sigma2 = scipy.stats.norm.fit(adj_diff_list)
n , bins2, _ = ax[0][1].hist(adj_diff_list, 25, density=True, color="orange")
best_fit_line2 = stats.norm.pdf(bins2, mu2, sigma2)
_ = ax[0][1].plot(bins2, best_fit_line2, color="red")
_ = ax[0][1].set_title("Adjacent different IDs \n $\mu$ = %.2f, median = %.2f, $\sigma$ = %.2f" % (mu2, np.median(adj_diff_list), sigma2))
_ = ax[0][1].set_ylabel("Probability Density")


_ = ax[1][1].hist(adj_diff_list, 25, color="orange")
_ = ax[1][1].set_xlabel("Weight")
_ = ax[1][1].set_ylabel("Count")



mu3, sigma3 = scipy.stats.norm.fit(adj_same_list)
n , bins3, _ = ax[0][2].hist(adj_same_list, 25, density=True, color="green")
best_fit_line3 = stats.norm.pdf(bins3, mu3, sigma3)
_ = ax[0][2].plot(bins3, best_fit_line3, color="red")
_ = ax[0][2].set_title("Adjacent same ID \n $\mu$ = %.2f, median = %.2f, $\sigma$ = %.2f" % (mu3, np.median(adj_same_list), sigma3))
_ = ax[0][2].set_ylabel("Probability Density")


_ = ax[1][2].hist(adj_same_list, 25, color="green")
_ = ax[1][2].set_xlabel("Weight")
_ = ax[1][2].set_ylabel("Count")


fig.suptitle("Main Graph")
plt.tight_layout()

plt.savefig(res_dir + '/figures/edge_distribution_G.png')

How to Traverse the Graph to explore all the nodes related to the North Pole Coronal Hole?¶

Step 1: First find the coronal holes that are in similar regions in the first and last frame.

In [15]:

for node in G.nodes:
    if G.nodes[node]["id"] == 4 and G.nodes[node]["frame_num"] == 1000:
        print(G.nodes[node])
        print(node)
        
    if G.nodes[node]["id"] == 4 and G.nodes[node]["frame_num"] == 1:
        print(G.nodes[node])
        print(node)

{'area': 0.1777568, 'id': 4, 'frame_num': 1, 'frame_timestamp': Timestamp('2007-03-01 00:01:48'), 'count': 0, 'color': [244, 230, 94], 'net_flux': 0.5546990722427789, 'abs_flux': 0.559596492926508}
1_4_0
{'area': 0.41803846, 'id': 4, 'frame_num': 1000, 'frame_timestamp': Timestamp('2007-06-03 10:01:24'), 'count': 0, 'color': [244, 230, 94], 'net_flux': 1.3488836106873416, 'abs_flux': 1.3749790336408365}
1000_4_0

Define the start and end nodes.

In [16]:

start = "1_4_0"
end = "1000_4_0"

Step 2: Create a directed version of our main Graph

In [17]:

DiG = G.to_directed()

Step 3: Find the shortest path by using the Dijkstra's algorithm

In [18]:

# define the cost function
def func(n1, n2, d):
    if G.nodes[n1]["id"] == DiG.nodes[n2]["id"]:
        return 0
    else:
        weight = DiG.get_edge_data(n1, n2)["weight"]
        return 1 - weight

In [19]:

# apply dijkstra's algorithm
path1 = nx.algorithms.shortest_paths.weighted.dijkstra_path(DiG, source=start, target=end, weight=func)

In [20]:

set_of_classes = set()

In [21]:

for node in path1:
    #print(str(node) + ", id: " + str(G.nodes[node]["id"]) + ", frame: "+ str(G.nodes[node]["frame_num"]) )
    set_of_classes.add(G.nodes[node]["id"])

In [22]:

list(set_of_classes)

Out[22]:

[4]

In [23]:

path_g = G.subgraph(path1)

In [24]:

unique_id_class = []

In [25]:

for id_class in list(set_of_classes):
    for node in path_g.nodes:
        if path_g.nodes[node]["id"] == id_class:
            unique_id_class.append(node)
            break

In [26]:

unique_id_class

Out[26]:

['766_4_0']

In [27]:

# set time array for plotting purposes
timearray = np.arange('2007-03-01T00', '2007-06-03T10',np.timedelta64(2*68,'m'), dtype='datetime64')

In [28]:

len(timearray)

Out[28]:

1000

In [29]:

fig, ax = plt.subplots(figsize=(10, 10))
# draw graph, nodes positions are based on their count and frame_num.
# labels are the coronal hole id number.
pos, labels = graph.get_plot_features(sub_graph=path_g)

edge_weights = nx.get_edge_attributes(G=path_g, name='weight')
edges, weights = zip(*edge_weights.items())

# plot nodes and labels.
nx.draw(path_g, pos=pos, font_weight='bold', ax=ax, node_size=700,
        node_color=[c.to_rgba(np.array(path_g.nodes[ch]["color"]) / 255)
                    for ch in path_g.nodes], edgelist=[])

for ch in unique_id_class:
    ax.scatter([],[], c=[c.to_rgba(np.array(path_g.nodes[ch]["color"]) / 255)], label='ID {}'.format(path_g.nodes[ch]["id"]))

nx.draw_networkx_edges(path_g, pos=pos, edge_color=weights, edgelist=edges,
                                        edge_cmap=plt.cm.get_cmap('Greys'), edge_vmin=0,
                                        edge_vmax=1, width=3, ax=ax)

nx.draw_networkx_edge_labels(G=path_g, pos=pos, edge_labels=edge_weights, ax=ax, alpha=1, font_size=10)

# Only show ticks on the left and bottom spines
ax.yaxis.set_ticks_position('left')
ax.set_xlim(tuple(sum(i) for i in zip(ax.get_xlim(), (-0.5, 0.5))))


# set y ticks
ax.yaxis.get_major_locator().set_params(integer=True)
ax.set_yticks(np.linspace(0, len(timearray), len(timearray))[::150])
ax.set_yticklabels(timearray[::150], fontsize=12)
ax.tick_params(axis='y', rotation=40)

ax.axis('on')
_ = ax.set_title("Coronal Hole Connectivity North")
_ = plt.gca().legend()

plt.savefig(res_dir + '/figures/north_ch_shortest_path_res_4to4.png')

Spatio-temporal Analysis of the North pole - Area and Flux vs. Time¶

First, lets examine the area over time of the north coronal hole.

In [30]:

frame_array = np.arange(1, 1000)

In [31]:

fig, ax = plt.subplots(figsize=(15, 5))


for frame in frame_array:
    holder = []
    for node in path_g:
        if path_g.nodes[node]["frame_num"] == frame:
            holder.append(node)
            
    area = 0 
    for node in holder:
        area += path_g.nodes[node]["area"]
    
    if area != 0:
        ax.scatter(frame, area, c=[c.to_rgba(np.array(path_g.nodes[holder[0]]["color"]) / 255)])
        
for ch in unique_id_class:
    ax.scatter([],[], c=[c.to_rgba(np.array(path_g.nodes[ch]["color"]) / 255)], label='ID {}'.format(path_g.nodes[ch]["id"]))

# set x ticks to be timestamps
ax.set_xticks(np.linspace(0, len(timearray), len(timearray))[::100])
ax.set_xticklabels(timearray[::100], fontsize=12)
ax.tick_params(axis='x', rotation=15)

# label axis
_ = ax.set_ylabel("Area ($R_{\odot}^2$)")
_ = ax.set_title("North Pole Coronal Hole Shortest Path Spatiotemporal Properties")
_ = plt.legend()
plt.savefig(res_dir + '/figures/north_ch_shortest_path_area.png')

In [32]:

fig, ax = plt.subplots(figsize=(15, 5))

for frame in frame_array:
    holder = []
    for node in path_g:
        if path_g.nodes[node]["frame_num"] == frame:
            holder.append(node)
            
    # initialize net flux and absolute flux.         
    nf, af = 0, 0
    for node in holder:
        nf += path_g.nodes[node]["net_flux"]
        af += path_g.nodes[node]["abs_flux"]
        
    if nf !=0:
        ax.scatter(frame, nf, c="b")
    if af !=0:
        ax.scatter(frame, af, c="r")

# add label to net and absolute flux. 
ax.scatter([],[], c="b", label="Net-Flux")
ax.scatter([], [], c="r", label="Abs-Flux")

# set x ticks to be timestamps
ax.set_xticks(np.linspace(0, len(timearray), len(timearray))[::100])
ax.set_xticklabels(timearray[::100])
ax.tick_params(axis='x', rotation=15)

# label axis
_ = ax.set_ylabel("$\Phi_{B}$")
_ = plt.legend()
_ = ax.set_title("North Pole Coronal Hole Shortest Path Flux")
plt.savefig(res_dir + '/figures/north_ch_shortest_path_net_flux.png')

North Shortest Path and Neighbors¶

Shorstest path described above and neighoring nodes that their edge connection is greater than a certain threshold.

In [33]:

THRESH = 0.5

In [34]:

north_and_ne = path_g.copy()

In [35]:

# add all nodes that are adjacent to nodes in shortest path
for node in path_g:
    for adj_node in G[node]:
        # if the edge weight is greater than threshold, then add the neighbor node. 
        if G[node][adj_node]['weight'] >= THRESH:
            # check if the node already exists in shortest path
            if adj_node not in north_and_ne.nodes:
                # add the node.
                north_and_ne.add_node(str(adj_node), 
                            area=G.nodes[adj_node]["area"],
                            id=G.nodes[adj_node]["id"],
                            frame_num=G.nodes[adj_node]["frame_num"],
                            frame_timestamp=G.nodes[adj_node]["frame_timestamp"],
                            count=G.nodes[adj_node]["count"],
                            color=G.nodes[adj_node]["color"], 
                            net_flux=G.nodes[adj_node]["net_flux"], 
                            abs_flux=G.nodes[adj_node]["abs_flux"])

In [36]:

for u,v,w in G.edges.data():
    # if the two nodes exist in the shorstest path and neighbors then add the edge. 
    if u in north_and_ne and v in north_and_ne:
        if not north_and_ne.has_edge(u, v):
            north_and_ne.add_edge(u, v, weight=w)

In [37]:

ne_classes = set()

for node in north_and_ne:
    ne_classes.add(G.nodes[node]["id"])

In [38]:

str(ne_classes)

Out[38]:

'{49, 4, 54}'

In [39]:

fig, ax = plt.subplots(figsize=(15, 5))


for frame in frame_array:
    holder = []
    for node in north_and_ne:
        if north_and_ne.nodes[node]["frame_num"] == frame:
            holder.append(node)
            
    area = 0 
    for node in holder:
        area += north_and_ne.nodes[node]["area"]
    
    if area != 0:
        ax.scatter(frame, area, c=['r' if len(holder) > 1 else 'b'])
        
# set color label        
ax.scatter([], [], c="r", label="npf $>$ 1")
ax.scatter([], [], c="b", label="npf $\leq$ 1")


# set x ticks to be timestamps
ax.set_xticks(np.linspace(0, len(timearray), len(timearray))[::150])
ax.set_xticklabels(timearray[::150])
ax.tick_params(axis='x', rotation=15)

_ = ax.set_ylabel("Area ($R_{\odot}^2$)")
# _ = ax.set_title("North Pole Coronal Hole Shortest Path and Neighbors Spatiotemporal Properties \n Edge Threshold is " + str(THRESH) + "\n " + str(ne_classes))
_ = ax.set_title("North Pole Coronal Hole Shortest Path and Neighbors Spatiotemporal Properties \n Edge Threshold is " + str(THRESH))
_ = plt.legend()
plt.savefig(res_dir + '/figures/north_ch_shortest_path_and_ne_area' + str(THRESH) + '.png', bbox_inches='tight')

In [40]:

fig, ax = plt.subplots(figsize=(15, 5))

for frame in frame_array:
    holder = []
    for node in north_and_ne:
        if north_and_ne.nodes[node]["frame_num"] == frame:
            holder.append(node)
            
    # initialize net flux and absolute flux.         
    nf, af = 0, 0
    for node in holder:
        nf += north_and_ne.nodes[node]["net_flux"]
        af += north_and_ne.nodes[node]["abs_flux"]
        
    if nf !=0:
        ax.scatter(frame, nf, c="b")
    if af !=0:
        ax.scatter(frame, af, c="r")

# add label to net and absolute flux. 
ax.scatter([],[], c="b", label="Net-Flux")
ax.scatter([], [], c="r", label="Abs-Flux")

# set x ticks to be timestamps
ax.set_xticks(np.linspace(0, len(timearray), len(timearray))[::150])
ax.set_xticklabels(timearray[::150])
ax.tick_params(axis='x', rotation=15)

# label axis
_ = ax.set_ylabel("$\Phi_{B}$")
_ = plt.legend()
_ = ax.set_title("North Pole Coronal Hole Shortest Path and Neighbors Flux \n Edge Threshold is " + str(THRESH))
plt.savefig(res_dir + '/figures/north_ch_shortest_path_and_ne_flux' + str(THRESH) + '.png', bbox_inches='tight')

North shortest path and neighbors classes¶

Shorstest path described above and neighoring nodes and all nodes with the same ID that their edge connection is greater than a certain threshold.

In [41]:

north_and_ne_c = north_and_ne.copy()

In [42]:

# add all nodes that are adjacent to nodes in shortest path
for node in G:
    if node not in north_and_ne_c.nodes and G.nodes[node]["id"] in ne_classes:
        # add the node.
        north_and_ne_c.add_node(str(node), 
                    area=G.nodes[node]["area"],
                    id=G.nodes[node]["id"],
                    frame_num=G.nodes[node]["frame_num"],
                    frame_timestamp=G.nodes[node]["frame_timestamp"],
                    count=G.nodes[node]["count"],
                    color=G.nodes[node]["color"], 
                    net_flux=G.nodes[node]["net_flux"], 
                    abs_flux=G.nodes[node]["abs_flux"])

In [43]:

for u,v,w in G.edges.data():
    # if the two nodes exist in the shorstest path and neighbors then add the edge. 
    if u in north_and_ne_c and v in north_and_ne_c:
        if not north_and_ne_c.has_edge(u, v):
            north_and_ne_c.add_edge(u, v, weight=w)

In [44]:

fig, ax = plt.subplots(figsize=(15, 5))


for frame in frame_array:
    holder = []
    for node in north_and_ne_c:
        if north_and_ne_c.nodes[node]["frame_num"] == frame:
            holder.append(node)
            
    area = 0 
    for node in holder:
        area += north_and_ne_c.nodes[node]["area"]
        
    #if len(holder) == 0:
    #    print("no nodes in this frame", frame)
    
    if area != 0:
        ax.scatter(frame, area, c=['r' if len(holder) > 1 else 'b'])

ax.scatter([], [], c="r", label="npf $>$ 1")
ax.scatter([], [], c="b", label="npf $\leq$ 1")

# set x ticks to be timestamps
ax.set_xticks(np.linspace(0, len(timearray), len(timearray))[::150])
ax.set_xticklabels(timearray[::150])
ax.tick_params(axis='x', rotation=15)

_ = ax.set_ylabel("Area ($R_{\odot}^2$)")
_ = ax.set_title("North Pole Coronal Hole Shortest Path and Neighbors Classes Spatiotemporal Properties \n Edge Threshold is " + str(THRESH))
_ = plt.legend()
plt.savefig(res_dir + '/figures/north_ch_shortest_path_and_ne_c_area' + str(THRESH) + '.png', bbox_inches='tight')

In [45]:

fig, ax = plt.subplots(figsize=(15, 5))

for frame in frame_array:
    holder = []
    for node in north_and_ne_c:
        if north_and_ne_c.nodes[node]["frame_num"] == frame:
            holder.append(node)
            
    # initialize net flux and absolute flux.         
    nf, af = 0, 0
    for node in holder:
        nf += north_and_ne_c.nodes[node]["net_flux"]
        af += north_and_ne_c.nodes[node]["abs_flux"]
        
    if nf !=0:
        ax.scatter(frame, nf, c="b")
    if af !=0:
        ax.scatter(frame, af, c="r")

# add label to net and absolute flux. 
ax.scatter([],[], c="b", label="Net-Flux")
ax.scatter([], [], c="r", label="Abs-Flux")

# set x ticks to be timestamps
ax.set_xticks(np.linspace(0, len(timearray), len(timearray))[::150])
ax.set_xticklabels(timearray[::150])
ax.tick_params(axis='x', rotation=15)

# label axis
_ = ax.set_ylabel("$\Phi_{B}$")
_ = plt.legend()
_ = ax.set_title("North Pole Coronal Hole Shortest Path and Neighbors Classes Flux \n Edge Threshold is " + str(THRESH))
plt.savefig(res_dir + '/figures/north_ch_shortest_path_and_ne_c_flux' + str(THRESH) + '.png', bbox_inches='tight')

In [46]:

print("Total number of nodes shortest path + neighbors: ", north_and_ne.number_of_nodes())
print("Total number of edges: ", north_and_ne.number_of_edges())


print("Total number of nodes shortest path + neighbors and classes: ", north_and_ne_c.number_of_nodes())
print("Total number of edges: ", north_and_ne_c.number_of_edges())

Total number of nodes shortest path + neighbors:  1003
Total number of edges:  1003
Total number of nodes shortest path + neighbors and classes:  1027
Total number of edges:  1038

Create Video¶

In [47]:

# choose codec according to format needed
fourcc = cv2.VideoWriter_fourcc(*'mp4v')

video = cv2.VideoWriter(os.path.join(res_dir, "north_pole_ne_c" + str(THRESH) + ".mov"), fourcc, 40, (1200, 800))

for frame in frame_array:
    
    holder = []
    
    for node in north_and_ne_c:
        if north_and_ne_c.nodes[node]["frame_num"] == frame:
            holder.append(node)
            
    else:
        timestamp = None
        for node in G:
            if G.nodes[node]["frame_num"] == frame:
                timestamp = G.nodes[node]["frame_timestamp"]
                break
        
        if timestamp is not None:
            pickle_file = str(timestamp).replace(':', '-') + ".pkl"
            pickle_file = pickle_file.replace(" ", "-")
            frame_read = pickle.load(open(os.path.join(res_dir +  "/pkl/" + pickle_file), "rb"))

            ch_list = frame_read.contour_list
            
            for ch in ch_list:
                ch_id = str(frame) + "_" + str(ch.id) + "_" + str(ch.count)
                if ch_id not in list(north_and_ne_c.nodes):
                    ch.color = (169,169,169)

            plot_coronal_hole(ch_list, 160, 400, "North Coronal Hole \n Time: " + str(frame_read.timestamp), 
                              filename=os.path.join(res_dir, "frames_north.png"),
                              plot_rect=False, plot_circle=True, circle_radius=50, 
                              thickness_circle=1, thickness_rect=2, fontscale=0.3, 
                              origin=None, dpi=200)
            
            
            
            image_file_name = os.path.join(res_dir, "frames_north.png")
            img = cv2.imread(image_file_name)
            video.write(img)
            

cv2.destroyAllWindows()
video.release()