Fig_20.ipynb

Kernel: Python 3 (system-wide)

import os
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.transforms as mtransforms

from matplotlib import colors

from scipy import signal
from scipy.interpolate import InterpolatedUnivariateSpline

plt.style.use('../jfm.mplstyle')

!find ../Experimental_Dataset -type f -name "elevation.csv.gz" | xargs dirname | sort | grep "undular" | tee "data_directories.txt" > /dev/null

with open("data_directories.txt","r") as data_directories_file:
    data_directories = data_directories_file.read().splitlines()

x_loc = "-100.0cm"
#x_loc = "-50.00cm"

data_dir = sorted([d for d in data_directories if x_loc in d])[0]

in_to_mm = 1/25.4

linestyles = ['-','--','-.',':']
line_list = []

time  = np.loadtxt(f"{data_dir}/time.csv",delimiter=',',dtype=int)
space = np.loadtxt(f"{data_dir}/space.csv",delimiter=',',dtype=int)
elev  = np.loadtxt(f"{data_dir}/elevation.csv",delimiter=',')

c_t1   = np.loadtxt(f"{data_dir}/crests_1st.csv",delimiter=',',dtype=int,usecols=(0,))
c_t2   = np.loadtxt(f"{data_dir}/crests_2nd.csv",delimiter=',',dtype=int,usecols=(0,))
c_t3   = np.loadtxt(f"{data_dir}/crests_3rd.csv",delimiter=',',dtype=int,usecols=(0,))
c_t4   = np.loadtxt(f"{data_dir}/crests_4th.csv",delimiter=',',dtype=int,usecols=(0,))

c_v1   = np.loadtxt(f"{data_dir}/crests_1st.csv",delimiter=',',usecols=(1,))
c_v2   = np.loadtxt(f"{data_dir}/crests_2nd.csv",delimiter=',',usecols=(1,))
c_v3   = np.loadtxt(f"{data_dir}/crests_3rd.csv",delimiter=',',usecols=(1,))
c_v4   = np.loadtxt(f"{data_dir}/crests_4th.csv",delimiter=',',usecols=(1,))

h0 = 50
a0 = 5
g  = 9810
c0 = np.sqrt(g*h0)
t0 = h0/c0

eta = elev/a0
y = space*0.3/h0
t = time*0.008/t0

c1 = c_t1*0.008/t0
c2 = c_t2*0.008/t0
c3 = c_t3*0.008/t0
c4 = c_t4*0.008/t0

c1v = c_v1/a0 #- 1
c2v = c_v2/a0 #- 1
c3v = c_v3/a0 #- 1
c4v = c_v4/a0 #- 1

T,Y = np.meshgrid(t,y,indexing='ij')

def smooth(lag):
    # low pass filtering
    w = signal.get_window('bartlett',30)
    w /= np.sum(w)
    lag = signal.convolve(lag,w,'same')
    for _ in range(3):
        lag = signal.convolve(lag,w,'same')

    return lag

mask = np.empty(eta.shape,dtype='uint8')
mask[eta > 0.9] = 1

# note that the connected components is slightly inconsistent because
# of random seed points to start grouping. Fixed starting points could
# add consistency. Easier to just try again.
num_groups, labels, stats, centroids = \
    cv.connectedComponentsWithStats(mask,connectivity=8)

# ignore the background image, note that last is largest
idx_candidates = np.argsort(stats[:,cv.CC_STAT_AREA])[:-1]
# sort by arrival
ic_sort_order = np.argsort(stats[idx_candidates,cv.CC_STAT_TOP])

nums_sorted = idx_candidates[ic_sort_order]    

arr   = np.zeros(eta.shape,dtype=int)
s = y.shape + (len(nums_sorted),)
c_i = np.zeros(s,dtype=int  )
c_t = np.zeros(s,dtype=float)
c_v = np.zeros(s,dtype=float)
for i,n in enumerate(nums_sorted):
    group = labels == n
    arr[group] = i + 1
    
    invalid_locs = np.sum(group,axis=0) < 1
    
    c_i[:,i] = np.argmax(eta*group,axis=0)
    c_i[invalid_locs,i] = 0
    
    c_t[:,i] = t[c_i[:,i]]
    c_t[invalid_locs,i] = np.nan

    crest = np.empty(c_v[:,i].shape)
    for y_idx, t_idx in enumerate(c_i[:,i]):
        crest[y_idx] = eta[t_idx,y_idx]
    
    c_v[:,i] = crest
    c_v[invalid_locs,i] = np.nan


fig, axes = plt.subplots(nrows=2,figsize=(5.0*0.8,3.0),layout="constrained",sharex='all')

ax = axes[0]

cmap = ax.pcolormesh(Y,T,eta,vmin=0,vmax=2,cmap='coolwarm')
cbar = fig.colorbar(cmap,label=r'$\eta/a_0$',
                   location='top',orientation='horizontal')

ls  = ['-','--','-.',':',(0,(5,1,1,1,1,1))]
lbl = ['1st','2nd','3rd','4th','5th']
rng = range(len(ls))
for i,linestyle,label in zip(rng,ls,lbl):
    ax.plot(y,smooth(c_t[:,i]),'k',linestyle=linestyle,label=label)
    
ax.set_xlim(-2.5,12)
ax.set_ylim(45,85)

ax.set_ylabel(r'$t$')

ax = axes[1]
ax.axhline(1,color='k',linewidth=0.5)

for i,linestyle,label in zip(rng,ls,lbl):
    ax.plot(y,smooth(c_v[:,i]),'k',linestyle=linestyle)

#for label, ax in zip(['a)', 'b)'], axes.flat):
#    trans = mtransforms.ScaledTranslation(-29/72, 7/72, fig.dpi_scale_trans)
#    ax.text(0.0, 1.0, label, transform=ax.transAxes + trans,
#            fontsize='9.0', verticalalignment='top', fontfamily='Times New Roman',
#            bbox=dict(facecolor='1.0', edgecolor='none', pad=2.0))

ax.set_ylim(0.7,2.0)
ax.set_xlim(-2.5,12)

ax.set_ylabel(r'$a/a_0$')
ax.set_xlabel(r'$y$')

fig.legend(bbox_to_anchor=(1, 0.5), loc="center left", title='crest')

plt.savefig('Fig_20.jpg', bbox_inches="tight")