In [3]:

    

import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
from scipy.interpolate import interp2d
import sys
import argparse
import os
import numpy as np
from numpy.random import uniform
import pandas as pd
from itertools import product
import datetime
import glob
import re
from numpy import ma
%matplotlib inline
# %matplotlib notebook
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
# plt.rcParams['figure.figsize'] = (10,5) 
# plt.rcParams['figure.figsize'] = (10,5.625)   # 16:9
plt.rcParams['figure.figsize'] = (10,6.180)    #golden ratio
# plt.rcParams['figure.figsize'] = (10*2,6.180*2)    #golden ratio


    

In [4]:

    

# https://bougui505.github.io/2016/08/31/compute_the_shortest_path_on_a_grid_using_python.html
def dijkstra(V):
    mask = V.mask
    visit_mask = mask.copy() # mask visited cells
    m = np.ones_like(V) * np.inf
    connectivity = [(i,j) for i in [-1, 0, 1] for j in [-1, 0, 1] if (not (i == j == 0))]
    cc = np.unravel_index(V.argmin(), m.shape) # current_cell
    m[cc] = 0
    P = {}  # dictionary of predecessors 
    #while (~visit_mask).sum() > 0:
    for _ in range(V.size):
        #print cc
        neighbors = [tuple(e) for e in np.asarray(cc) - connectivity 
                     if e[0] > 0 and e[1] > 0 and e[0] < V.shape[0] and e[1] < V.shape[1]]
        neighbors = [ e for e in neighbors if not visit_mask[e] ]
        tentative_distance = np.asarray([V[e]-V[cc] for e in neighbors])
        for i,e in enumerate(neighbors):
            d = tentative_distance[i] + m[cc]
            if d < m[e]:
                m[e] = d
                P[e] = cc
        visit_mask[cc] = True
        m_mask = ma.masked_array(m, visit_mask)
        cc = np.unravel_index(m_mask.argmin(), m.shape)
    return m, P


    

In [5]:

    

def shortestPath(start, end, P):
    Path = []
    step = end
    while 1:
        Path.append(step)
        if step == start: break
        step = P[step]
    Path.reverse()
    return np.asarray(Path)


    

In [6]:

    

all_data = pd.read_feather("/Users/weilu/Research/data/pulling/19_Feb_data_1.feather")


    

In [7]:

    

test = all_data.query("temp=='400'").query("perturbation=='original'").query("mode=='2d_z_qw'")


    

In [8]:

    

x=1
y=2
z=3
zmin = 0
zmax=20
data = test[["index", "x","y","f"]].values
#         print(data)
data = data[~np.isnan(data).any(axis=1)] # remove rows with nan
data = data[~(data[:,z] > zmax)] # remove rows of data for z not in [zmin zmax]
data = data[~(data[:,z] < zmin)]
size = 20
xi = np.linspace(min(data[:,x]), max(data[:,x]), size)
yi = np.linspace(min(data[:,y]), max(data[:,y]), size)
zi = griddata((data[:,x], data[:,y]), data[:,z], (xi[None,:], yi[:,None]), method='linear')


    

In [13]:

    

P


    

    
Out[13]:





{(1, 2): (2, 3),
 (1, 3): (2, 4),
 (1, 4): (2, 4),
 (1, 5): (2, 4),
 (1, 6): (2, 5),
 (1, 7): (2, 8),
 (1, 8): (2, 8),
 (2, 1): (2, 2),
 (2, 2): (2, 3),
 (2, 3): (3, 4),
 (2, 4): (3, 4),
 (2, 5): (3, 4),
 (2, 6): (3, 5),
 (2, 7): (3, 6),
 (2, 8): (3, 7),
 (2, 9): (2, 8),
 (2, 10): (3, 9),
 (2, 11): (3, 10),
 (2, 12): (3, 11),
 (2, 13): (2, 12),
 (3, 1): (2, 2),
 (3, 2): (2, 3),
 (3, 3): (4, 4),
 (3, 4): (4, 4),
 (3, 5): (4, 4),
 (3, 6): (3, 5),
 (3, 7): (3, 6),
 (3, 8): (3, 7),
 (3, 9): (2, 8),
 (3, 10): (3, 9),
 (3, 11): (4, 10),
 (3, 12): (4, 13),
 (3, 13): (4, 14),
 (3, 14): (4, 15),
 (3, 15): (4, 15),
 (3, 16): (4, 15),
 (3, 17): (4, 16),
 (4, 2): (3, 3),
 (4, 3): (4, 4),
 (4, 4): (5, 5),
 (4, 5): (5, 5),
 (4, 6): (5, 5),
 (4, 7): (3, 6),
 (4, 8): (3, 7),
 (4, 9): (3, 9),
 (4, 10): (3, 9),
 (4, 11): (5, 12),
 (4, 12): (5, 13),
 (4, 13): (5, 14),
 (4, 14): (5, 13),
 (4, 15): (5, 16),
 (4, 16): (5, 15),
 (4, 17): (5, 16),
 (4, 18): (5, 17),
 (4, 19): (5, 18),
 (5, 2): (4, 3),
 (5, 3): (4, 4),
 (5, 4): (6, 5),
 (5, 5): (6, 5),
 (5, 6): (6, 5),
 (5, 7): (6, 6),
 (5, 8): (5, 9),
 (5, 9): (4, 10),
 (5, 10): (6, 11),
 (5, 11): (6, 12),
 (5, 12): (6, 13),
 (5, 13): (6, 14),
 (5, 14): (6, 14),
 (5, 15): (6, 14),
 (5, 16): (6, 15),
 (5, 17): (5, 16),
 (5, 18): (6, 17),
 (5, 19): (5, 18),
 (6, 3): (5, 4),
 (6, 4): (6, 5),
 (6, 5): (7, 6),
 (6, 6): (7, 6),
 (6, 7): (7, 6),
 (6, 8): (5, 9),
 (6, 9): (5, 10),
 (6, 10): (6, 11),
 (6, 11): (7, 12),
 (6, 12): (7, 13),
 (6, 13): (7, 14),
 (6, 14): (7, 14),
 (6, 15): (7, 14),
 (6, 16): (7, 15),
 (6, 17): (7, 16),
 (6, 18): (7, 17),
 (6, 19): (7, 18),
 (7, 4): (6, 5),
 (7, 5): (8, 6),
 (7, 6): (8, 6),
 (7, 7): (8, 6),
 (7, 8): (8, 7),
 (7, 9): (8, 10),
 (7, 10): (8, 11),
 (7, 11): (8, 12),
 (7, 12): (8, 13),
 (7, 13): (8, 13),
 (7, 14): (8, 13),
 (7, 15): (7, 14),
 (7, 16): (7, 15),
 (7, 17): (7, 16),
 (7, 18): (7, 17),
 (7, 19): (7, 18),
 (8, 4): (7, 5),
 (8, 5): (9, 6),
 (8, 6): (9, 6),
 (8, 7): (9, 6),
 (8, 8): (8, 7),
 (8, 9): (9, 10),
 (8, 10): (8, 11),
 (8, 11): (8, 12),
 (8, 12): (8, 13),
 (8, 13): (9, 14),
 (8, 14): (9, 14),
 (8, 15): (9, 14),
 (8, 16): (7, 15),
 (8, 17): (7, 16),
 (8, 18): (7, 17),
 (8, 19): (7, 18),
 (9, 5): (9, 6),
 (9, 6): (10, 7),
 (9, 7): (10, 7),
 (9, 8): (10, 7),
 (9, 9): (10, 10),
 (9, 10): (8, 11),
 (9, 11): (10, 12),
 (9, 12): (10, 13),
 (9, 13): (10, 13),
 (9, 14): (10, 13),
 (9, 15): (9, 14),
 (9, 16): (9, 15),
 (9, 17): (8, 16),
 (9, 18): (8, 17),
 (9, 19): (9, 18),
 (10, 6): (11, 7),
 (10, 7): (11, 7),
 (10, 8): (11, 7),
 (10, 9): (11, 10),
 (10, 10): (11, 11),
 (10, 11): (11, 12),
 (10, 12): (11, 13),
 (10, 13): (11, 13),
 (10, 14): (11, 13),
 (10, 15): (11, 14),
 (10, 16): (11, 15),
 (10, 17): (11, 17),
 (10, 18): (11, 17),
 (10, 19): (11, 18),
 (11, 6): (11, 7),
 (11, 7): (12, 8),
 (11, 8): (12, 8),
 (11, 9): (12, 10),
 (11, 10): (12, 10),
 (11, 11): (12, 12),
 (11, 12): (12, 13),
 (11, 13): (12, 13),
 (11, 14): (12, 13),
 (11, 15): (12, 14),
 (11, 16): (12, 15),
 (11, 17): (12, 16),
 (11, 18): (12, 17),
 (11, 19): (12, 18),
 (12, 7): (12, 8),
 (12, 8): (13, 9),
 (12, 9): (13, 10),
 (12, 10): (13, 11),
 (12, 11): (13, 12),
 (12, 12): (13, 13),
 (12, 13): (13, 13),
 (12, 14): (13, 13),
 (12, 15): (13, 14),
 (12, 16): (13, 15),
 (12, 17): (13, 16),
 (12, 18): (13, 17),
 (12, 19): (13, 18),
 (13, 8): (14, 9),
 (13, 9): (14, 10),
 (13, 10): (14, 11),
 (13, 11): (14, 12),
 (13, 12): (14, 13),
 (13, 13): (14, 13),
 (13, 14): (14, 13),
 (13, 15): (14, 14),
 (13, 16): (14, 15),
 (13, 17): (14, 16),
 (13, 18): (13, 17),
 (13, 19): (14, 18),
 (14, 8): (14, 9),
 (14, 9): (15, 10),
 (14, 10): (15, 11),
 (14, 11): (15, 12),
 (14, 12): (15, 13),
 (14, 13): (15, 13),
 (14, 14): (15, 13),
 (14, 15): (15, 14),
 (14, 16): (15, 15),
 (14, 17): (15, 16),
 (14, 18): (15, 17),
 (14, 19): (15, 18),
 (15, 9): (15, 10),
 (15, 10): (15, 11),
 (15, 11): (15, 12),
 (15, 12): (15, 13),
 (15, 14): (15, 13),
 (15, 15): (15, 14),
 (15, 16): (15, 15),
 (15, 17): (15, 16),
 (15, 18): (15, 17),
 (15, 19): (14, 18),
 (16, 10): (15, 11),
 (16, 11): (15, 12),
 (16, 12): (15, 13),
 (16, 13): (15, 13),
 (16, 14): (15, 13),
 (16, 15): (15, 14),
 (16, 16): (15, 15),
 (16, 17): (15, 16),
 (16, 18): (15, 17),
 (16, 19): (16, 18),
 (17, 10): (16, 11),
 (17, 11): (18, 12),
 (17, 12): (16, 13),
 (17, 13): (16, 13),
 (17, 14): (16, 13),
 (17, 15): (16, 14),
 (17, 16): (18, 15),
 (17, 17): (17, 16),
 (18, 11): (17, 12),
 (18, 12): (17, 13),
 (18, 13): (17, 13),
 (18, 14): (17, 13),
 (18, 15): (17, 14),
 (18, 16): (17, 15),
 (19, 13): (18, 13),
 (19, 14): (18, 13)}

In [10]:

    

zi = np.where(np.isnan(zi), 50, zi)
V = ma.masked_array(zi, zi>40)
D, P = dijkstra(V)
source = np.unravel_index(V.argmin(), V.shape)
# source = (17,35)
path = shortestPath(source, (1,2), P)
plt.contourf(xi, yi, V, 30, cmap='jet')
# plt.plot(path[:,1], path[:,0], 'r.-')
plt.plot(xi[path[:,1]], yi[path[:,0]], 'r.-')


    

    
Out[10]:





[<matplotlib.lines.Line2D at 0x10c8a8518>]






    

In [12]:

    

f_on_path = [zi[tuple(p)] for p in path]
plt.plot(f_on_path)


    

    
Out[12]:





[<matplotlib.lines.Line2D at 0x181885fe10>]






    

In [ ]: