In [60]:
from PIL import Image
import numpy as np
import pandas as pd
import sklearn as sk
from sklearn import cluster
from sklearn import linear_model
from scipy.optimize import minimize
from sklearn.cluster import KMeans
import itertools
import os
from itertools import cycle
import scipy.stats as stats
from scipy.misc import factorial
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import bisect
import time
import pylab as pl
from IPython import display
import copy
import seaborn as sns
import cv2
import statsmodels.api as sm
%matplotlib inline
sns.set_style('whitegrid')
sns.set_context('poster')
In [61]:
threshold = 15 # 15 pixles out of place
# A log roughlt corresponds to 25 pixles
In [62]:
y_max = 1080 - 864
x_max = 384
def load_image( infilename ) :
img = Image.open(infilename)
img.load()
data = np.asarray( img, dtype="int32" )
return data
def save_image( npdata, outfilename ) :
img = Image.fromarray( np.asarray( np.clip(npdata,0,255), dtype="uint8"), "L" )
img.save( outfilename )
In [63]:
def extract_boundary(original,hsv_image, lower, upper, flag):
# need end points of the boundary too
mask = cv2.inRange(hsv_image, lower, upper)
# Bitwise-AND mask and original image
res = cv2.bitwise_and(original,original,mask= mask)
#boundaries in gray scale
gray = cv2.cvtColor(res,cv2.COLOR_BGR2GRAY)
# Otsu's thresholding and gaussian filtering to make the logs white and the background black for better detection
ret2,th2 = cv2.threshold(gray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
# Otsu's thresholding after Gaussian filtering
blur = cv2.GaussianBlur(gray,(5,5),0)
#logs will be white in th3
ret3,th3 = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
if(flag==1):
black, extLeft, extRight, cx,cy = find_contour(th3,original)
return black,extLeft,extRight,cx,cy
return th3
#just identify water flow path for drawing graphs
def detect_water(min_video_frame):
hsv = cv2.cvtColor(min_video_frame, cv2.COLOR_BGR2HSV)
# define range of green/yellow color in HSV
lower_green = np.array([29,86,6])
upper_green = np.array([64,255,255])
th3 = extract_boundary(min_video_frame,hsv,lower_green, upper_green,0)
store = th3
# morphing to get the skeletal structure/ medial line of the water flow
size = np.size(th3)
skel = np.zeros(th3.shape,np.uint8)
element = cv2.getStructuringElement(cv2.MORPH_CROSS,(3,3))
done = False
while(not done):
eroded = cv2.erode(th3,element)
temp = cv2.dilate(eroded,element)
temp = cv2.subtract(th3,temp)
skel = cv2.bitwise_or(skel,temp)
th3 = eroded.copy()
zeros = size - cv2.countNonZero(th3)
if zeros==size:
done = True
return store,skel
def detect_logs(min_video_frame):
hsv = cv2.cvtColor(min_video_frame, cv2.COLOR_BGR2HSV)
# define range of blue color in HSV
lower_blue = np.array([110,50,50])
upper_blue = np.array([130,255,255])
th3 = extract_boundary(min_video_frame,hsv,lower_blue, upper_blue,0)
#smooth the logs (current version very fat lines)
image ,contours, heirarchy = cv2.findContours(th3,1,2)#cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
# print(contours)
#Draw log contour + bonding rects
colored = cv2.cvtColor(image,cv2.COLOR_GRAY2BGR)
count =0
black = np.zeros(colored.shape)
return image,[]
In [64]:
cap = cv2.VideoCapture('./2017-09-19T10_01_14-0.mov')
count = 0
logs_minus1 = []
diffs = []
iter_num = 0
images = []
i = 0
while cap.isOpened():
i+=1
if i < 6:
continue
ret,frame = cap.read()
if ret==True:
logs,centers = detect_logs(frame[864:1080,0:384])
logs = logs / 255 # get a 0 or a 1
if len(logs_minus1) <= 0:
logs_minus1 = logs
diff = np.linalg.norm(logs_minus1 - logs)
if diff > threshold:
diffs.append(iter_num)
images.append(frame[864:1080,0:384])
logs_minus1 = logs
iter_num += 1
else:
cap.set(1, count-1)
print("something not working")
cv2.waitKey(1000)
if cv2.waitKey(10) == 27:
break
if cap.get(1) == cap.get(7):
break
if count >= 20:
break
cap.release()
In [65]:
print("There are {number_changes} in the file, each occured at:\n{change_points}".format(
number_changes = len(diffs),
change_points = ["%0.2f"%(d/60) for d in diffs]))
In [66]:
fig, ax = plt.subplots(1,1, figsize=(15,2))
data_diffs = {}
for i in diffs:
ax.axvline(x=i/60, alpha=0.5)
ax.yaxis.set_visible(False)
ax.set_xlabel('Events occuring in log movements over time (min)')
data_diffs['df_1'] = diffs
plt.show()
In [67]:
num_rows=min(len(diffs)//2, 8)
fig, axes = plt.subplots(num_rows, 2, figsize=(15, 8*num_rows))
axes = axes.flatten()
for i, ax in enumerate(axes):
ax.set_title('Timestamp: %i'%diffs[i])
ax.imshow(images[i])
ax.set_axis_off()
fig.tight_layout()
plt.show()
In [68]:
cap = cv2.VideoCapture('./2017-09-09T12_08_00-0.mov')
count = 0
log_centers = []
mu, cov = None, None
logs_minus1 = []
diffs = []
iter_num = 0
images = []
while cap.isOpened():
ret,frame = cap.read()
if ret==True:
# store, skel = detect_water(frame[864:1080,0:384])
logs,centers = detect_logs(frame[864:1080,0:384])
logs = logs / 255 # get a 0 or a 1
if len(logs_minus1) <= 0:
logs_minus1 = logs
diff = np.linalg.norm(logs_minus1 - logs)
# print(diff)
if diff > threshold:
diffs.append(iter_num)
images.append(frame[864:1080,0:384])
# plt.imshow(logs)
logs_minus1 = logs
iter_num += 1
# cv2.imwrite("frames/frame_log%04d.jpg" % count, logs)
# count = count + 1
# log_centers.append(centers)
else:
# The next frame is not ready, so we try to read it again
cap.set(1, count-1)
print("something not working")
# It is better to wait for a while for the next frame to be ready
cv2.waitKey(1000)
if cv2.waitKey(10) == 27:
break
if cap.get(1) == cap.get(7):
# If the number of captured frames is equal to the total number of frames,
# we stop
break
if count >= 20:
break
cap.release()
In [69]:
print("There are {number_changes} in the file, each occured at:\n{change_points}".format(
number_changes = len(diffs),
change_points = ["%0.2f"%(d/60) for d in diffs]))
In [70]:
fig, ax = plt.subplots(1,1, figsize=(15,2))
for i in diffs:
ax.axvline(x=i/60, alpha=0.1)
ax.yaxis.set_visible(False)
ax.set_xlabel('Events occuring in log movements over time (min)')
data_diffs['df_2'] = diffs
plt.show()
In [71]:
num_rows=min(len(diffs)//2, 8)
fig, axes = plt.subplots(num_rows, 2, figsize=(15, 8*num_rows))
axes = axes.flatten()
for i, ax in enumerate(axes):
ax.set_title('Timestamp: %i'%diffs[i])
ax.imshow(images[i])
ax.set_axis_off()
fig.tight_layout()
plt.show()
In [72]:
cap = cv2.VideoCapture('./2017-09-01T13_16_45-2.mov')
count = 0
logs_minus1 = []
diffs = []
iter_num = 0
images = []
i = 0
while cap.isOpened():
i+=1
if i < 6:
continue
ret,frame = cap.read()
if ret==True:
logs,centers = detect_logs(frame[864:1080,0:384])
logs = logs / 255 # get a 0 or a 1
if len(logs_minus1) <= 0:
logs_minus1 = logs
diff = np.linalg.norm(logs_minus1 - logs)
if diff > threshold:
diffs.append(iter_num)
images.append(frame[864:1080,0:384])
logs_minus1 = logs
iter_num += 1
else:
cap.set(1, count-1)
print("something not working")
cv2.waitKey(1000)
if cv2.waitKey(10) == 27:
break
if cap.get(1) == cap.get(7):
break
if count >= 20:
break
cap.release()
In [73]:
print("There are {number_changes} in the file, each occured at:\n{change_points}".format(
number_changes = len(diffs),
change_points = ["%0.2f"%(d/60) for d in diffs]))
In [74]:
fig, ax = plt.subplots(1,1, figsize=(15,2))
for i in diffs:
ax.axvline(x=i/60, alpha=0.1)
ax.yaxis.set_visible(False)
ax.set_xlabel('Events occuring in log movements over time (min)')
data_diffs['df_3'] = diffs
plt.show()
In [75]:
num_rows=min(len(diffs)//2, 8)
fig, axes = plt.subplots(num_rows, 2, figsize=(15, 8*num_rows))
axes = axes.flatten()
for i, ax in enumerate(axes):
ax.set_title('Timestamp: %i'%diffs[i])
ax.imshow(images[i])
ax.set_axis_off()
fig.tight_layout()
plt.show()
In [76]:
cap = cv2.VideoCapture('./2017-09-06T09_35_14-1.mov')
count = 0
logs_minus1 = []
diffs = []
iter_num = 0
images = []
i = 0
while cap.isOpened():
i+=1
if i < 6:
continue
ret,frame = cap.read()
if ret==True:
logs,centers = detect_logs(frame[864:1080,0:384])
logs = logs / 255 # get a 0 or a 1
if len(logs_minus1) <= 0:
logs_minus1 = logs
diff = np.linalg.norm(logs_minus1 - logs)
if diff > threshold:
diffs.append(iter_num)
images.append(frame[864:1080,0:384])
logs_minus1 = logs
iter_num += 1
else:
cap.set(1, count-1)
print("something not working")
cv2.waitKey(1000)
if cv2.waitKey(10) == 27:
break
if cap.get(1) == cap.get(7):
break
if count >= 20:
break
cap.release()
In [77]:
print("There are {number_changes} in the file, each occured at:\n{change_points}".format(
number_changes = len(diffs),
change_points = ["%0.2f"%(d/60) for d in diffs]))
data_diffs['df_4'] = diffs
In [20]:
data_csvs = {}
data_csvs['df_1'] = pd.read_csv("./2017-09-19T10_01_14-0.csv").ix[6:]
data_csvs['df_2'] = pd.read_csv("./2017-09-09T12_08_00-0.csv").ix[6:]
data_csvs['df_3'] = pd.read_csv("./2017-09-01T13_16_45-2.csv").ix[6:]
data_csvs['df_4'] = pd.read_csv("./2017-09-06T09_35_14-1.csv").ix[6:]
for k,v in data_csvs.items():
data_csvs[k].columns = v.columns.str.strip()
data_csvs['df_1'].head(2)
Out[20]:
In [21]:
data_csvs['df_1'].columns
Out[21]:
In [22]:
water_columns = ['Desert_Water', 'Plains_Water', 'Wetlands_Water', 'Jungle_Water', 'Reservoir_Water', 'MountainValley_Water',
'Waterfall_Water', 'Floor_Water']
In [23]:
import scipy
from scipy.signal import butter, filtfilt
time_series = data_csvs['df_1'][water_columns[1]]
b,a = butter(10, 0.05, btype='low')
time_series.plot()
plt.plot(time_series.index.values, filtfilt(b, a, time_series), 'r--', alpha=0.5)
Out[23]:
In [24]:
time_series = data_csvs['df_2'][water_columns[1]]
b,a = butter(5, 0.1, btype='low')
time_series.plot()
plt.plot(time_series.index.values, filtfilt(b, a, time_series), 'r--', alpha=0.5)
Out[24]:
In [25]:
signal_data = {}
for k,v in data_csvs.items():
signal_data[k] = {}
signal_data[k]['index'] = v.index.values
for i, tsname in enumerate(water_columns):
for k,v in signal_data.items():
signal_data[k][tsname] = filtfilt(b, a, data_csvs[k][water_columns[i]])
Let's perform stepwise greedy selection to choose the best wavelets
In [26]:
def construct_x_matrix(best_splits, stop):
start = 0
for i,split in enumerate(best_splits):
next_part = np.arange(start, split, 1)
constant = np.ones_like(next_part)
if i == 0:
X = np.concatenate([
np.concatenate([constant, np.zeros(stop - split)]).reshape(-1,1),
np.concatenate([next_part, np.zeros(stop - split)]).reshape(-1,1)], axis=1)
start = split
continue
X = np.concatenate([
X,
np.concatenate([np.zeros(start), constant, np.zeros(stop - split)]).reshape(-1,1),
np.concatenate([np.zeros(start), next_part, np.zeros(stop - split)]).reshape(-1,1),
], axis=1)
start = split
return X
In [27]:
step_size = 2
# start with the entire dataset
time_series = signal_data['df_1']['Plains_Water']
start, stop = 0, len(time_series)
best_splits = [stop]
best_bics = []
store_best_splits = []
for i in range(25):
# naive approach is to iterate through at window size
bics = []
splits = list(range(start+step_size, stop, step_size))
for i in splits:
temp_splits = copy.deepcopy(best_splits)
bisect.insort(temp_splits, i)
X = construct_x_matrix(temp_splits, stop)
y = time_series
result = sm.OLS(endog=y, exog=X).fit()
bics.append(result.bic)
arg_best_bic = np.argmin(bics)
best_split = splits[arg_best_bic]
bisect.insort(best_splits, best_split)
best_bics.append(bics[arg_best_bic])
store_best_splits.append(copy.deepcopy(best_splits))
In [28]:
plt.plot(best_bics)
Out[28]:
In [29]:
# first_part = np.arange(start, best_split, 1)
# second_part = np.arange(best_split, stop, 1)
# X = np.concatenate([
# np.concatenate([np.ones_like(first_part), np.zeros_like(second_part)]).reshape(-1,1),
# np.concatenate([first_part, np.zeros_like(second_part)]).reshape(-1,1),
# np.concatenate([np.zeros_like(first_part), np.ones_like(second_part)]).reshape(-1,1),
# np.concatenate([np.zeros_like(first_part), second_part]).reshape(-1,1)
# ], axis=1)
# result = sm.OLS(endog=y, exog=X).fit()
# I am trying to select the next best split based on the p-values but
# it might make more sense to just make a split based on the length of the
# training data
# column_to_break = np.argmax((X[::2] > 0).sum())
# column_to_break = np.argmax(result.pvalues[::2] / (X[::2] > 0).sum())
# column_to_break
In [30]:
X = construct_x_matrix(store_best_splits[7], len(signal_data['df_1']['index']))
result = sm.OLS(endog=y, exog=X).fit()
index = signal_data['df_1']['index']
plt.scatter(index, result.predict())
plt.plot(time_series)
plt.show()
In [31]:
def get_best_split_points_and_plot_BICs(axes, df):
split_points = []
for k, ts_name in enumerate(water_columns):
split_points.append([])
time_series = signal_data[df][ts_name]
y = time_series
index = signal_data[df]['index']
step_size = 10
# start with the entire dataset
start, stop = 0, len(time_series)
best_splits = [stop]
best_bics = []
store_best_splits = []
for i in range(15):
# naive approach is to iterate through at window size
bics = []
splits = list(range(start+step_size, stop, step_size))
for i in splits:
temp_splits = copy.deepcopy(best_splits)
bisect.insort(temp_splits, i)
X = construct_x_matrix(temp_splits, stop)
result = sm.OLS(endog=y, exog=X).fit()
bics.append(result.bic)
arg_best_bic = np.argmin(bics)
best_split = splits[arg_best_bic]
bisect.insort(best_splits, best_split)
best_bics.append(bics[arg_best_bic])
store_best_splits.append(copy.deepcopy(best_splits))
split_points[k].append(store_best_splits)
axes[k].plot(best_bics)
axes[k].set_title(ts_name)
return split_points
In [32]:
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
signal_data['df_1']['best_splits'] = get_best_split_points_and_plot_BICs(axes, 'df_1')
plt.show()
In [40]:
def generate_colour():
while True:
yield tuple(np.random.uniform(0,1,3))
def plot_best_splits(axes, best_individual_splitpoints, df):
events = []
for k, ts_name in enumerate(water_columns):
time_series = signal_data[df][ts_name]
y = time_series
index = signal_data[df]['index']
start, stop = 0, len(time_series)
splits = signal_data[df]['best_splits'][k][0][best_individual_splitpoints[k]]
events.append(splits)
X = construct_x_matrix(splits, stop)
result = sm.OLS(endog=y, exog=X).fit()
axes[k].scatter(index/60, result.predict())
axes[k].plot(index/60, time_series)
axes[k].set_title(ts_name)
axes[k].set_ylim([0,1.5])
axes[k].set_xlim([0,np.max(index/60)])
return events
In [34]:
# best_individual_splitpoints = [
# 1, #MV
# 2, #P
# 3, #WF
# 0, #F
# 2, #R
# 1, #J
# 3, #W
# 1] #D
best_individual_splitpoints = [
2, #MV
2, #P
2, #WF
2, #F
2, #R
2, #J
2, #W
2] #D
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
events = plot_best_splits(axes, best_individual_splitpoints, 'df_1')
plt.show()
Now study the events over time to see if these match up with the observed actions from the log images.
In [35]:
# flat = []
# [[flat.append(i) for i in e] for e in events]
# flat = list(np.sort(flat))
# difference = np.diff(np.sort(flat))
# c_map = generate_colour()
# colours = []
# for i, d in enumerate(difference):
# if i == 0:
# colours.append(next(c_map))
# if d <= 20:
# colours.append(colours[i])
# else:
# colours.append(next(c_map))
# colours
# fig, ax = plt.subplots(1,1, figsize=(15,9))
# i = 0
# for event_log in events[::-1]:
# ax.axhline(y=i, c='black', lw=0.5, alpha=0.5)
# for e in event_log[:-1]:
# ax.axvline(x=e/60, ymin=(i+1)/(len(events)), ymax=(i)/(len(events)), c=colours[flat.index(e)], lw=5)
# i += 1
# ax.set_ylim(0,i)
# ax.grid(False)
# ax.set_yticks(np.arange(0,len(water_columns))+0.5)
# ax.set_yticklabels(water_columns[::-1], position=(0,-1))
# for i in data_diffs['df_1']:
# ax.axvline(x=i/60, alpha=0.1)
# plt.show()
focussed_events = events[-5::-1]
def generate_colour():
while True:
yield tuple(np.random.uniform(0,1,3))
flat = []
[[flat.append(i) for i in e] for e in focussed_events]
flat = list(np.sort(flat))
difference = np.diff(np.sort(flat))
c_map = generate_colour()
colours = []
for i, d in enumerate(difference):
if i == 0:
colours.append(next(c_map))
if d <= 20:
colours.append(colours[i])
else:
colours.append(next(c_map))
colours
fig, ax = plt.subplots(1,1, figsize=(15,4))
i = 0
for event_log in focussed_events:
ax.axhline(y=i, c='black', lw=0.5, alpha=0.5)
for e in event_log[:-1]:
ax.axvline(x=e/60, ymin=(i+1)/(len(focussed_events)), ymax=(i)/(len(focussed_events)), c=colours[flat.index(e)], lw=5)
i += 1
ax.set_ylim(0,i)
ax.grid(False)
ax.set_yticks(np.arange(0,len(water_columns[:4]))+0.5)
ax.set_yticklabels(water_columns[:4][::-1], position=(0,-1))
for i in data_diffs['df1']:
ax.axvline(x=i/60, alpha=0.1)
plt.show()
In [36]:
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
signal_data['df_2']['best_splits'] = get_best_split_points_and_plot_BICs(axes, 'df_2')
plt.show()
In [46]:
best_individual_splitpoints = [
10, #D
7, #P
7, #W
10, #J
4, #R
2, #MV
3, #WF
2] #F
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
events = plot_best_splits(axes, best_individual_splitpoints, 'df_2')
plt.show()
In [51]:
focussed_events = events[-5::-1]
def generate_colour():
while True:
yield tuple(np.random.uniform(0,1,3))
flat = []
[[flat.append(i) for i in e] for e in focussed_events]
flat = list(np.sort(flat))
difference = np.diff(np.sort(flat))
c_map = generate_colour()
colours = []
for i, d in enumerate(difference):
if i == 0:
colours.append(next(c_map))
if d <= 20:
colours.append(colours[i])
else:
colours.append(next(c_map))
colours
fig, ax = plt.subplots(1,1, figsize=(15,4))
i = 0
for event_log in focussed_events:
ax.axhline(y=i, c='black', lw=0.5, alpha=0.5)
for e in event_log[:-1]:
ax.axvline(x=e/60, ymin=(i+1)/(len(focussed_events)), ymax=(i)/(len(focussed_events)), c=colours[flat.index(e)], lw=5,
zorder=100)
i += 1
ax.set_ylim(0,i)
ax.grid(False)
ax.set_yticks(np.arange(0,len(water_columns[:4]))+0.5)
ax.set_yticklabels(water_columns[:4][::-1], position=(0,-1))
for i in data_diffs['df_2']:
ax.axvline(x=i/60, alpha=0.1, zorder=10)
plt.show()
In [52]:
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
signal_data['df_3']['best_splits'] = get_best_split_points_and_plot_BICs(axes, 'df_3')
plt.show()
In [88]:
best_individual_splitpoints = [
8, #D
8, #P
8, #W
8, #J
8, #R
8, #MV
8, #WF
8] #F
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
events = plot_best_splits(axes, best_individual_splitpoints, 'df_3')
plt.show()
In [89]:
focussed_events = events[-5::-1]
def generate_colour():
while True:
yield tuple(np.random.uniform(0,1,3))
flat = []
[[flat.append(i) for i in e] for e in focussed_events]
flat = list(np.sort(flat))
difference = np.diff(np.sort(flat))
c_map = generate_colour()
colours = []
for i, d in enumerate(difference):
if i == 0:
colours.append(next(c_map))
if d <= 20:
colours.append(colours[i])
else:
colours.append(next(c_map))
colours
fig, ax = plt.subplots(1,1, figsize=(15,4))
i = 0
for event_log in focussed_events:
ax.axhline(y=i, c='black', lw=0.5, alpha=0.5)
for e in event_log[:-1]:
ax.axvline(x=e/60, ymin=(i+1)/(len(focussed_events)), ymax=(i)/(len(focussed_events)), c=colours[flat.index(e)], lw=5, zorder=100)
i += 1
ax.set_ylim(0,i)
ax.grid(False)
ax.set_yticks(np.arange(0,len(water_columns[:4]))+0.5)
ax.set_yticklabels(water_columns[:4][::-1], position=(0,-1))
for i in data_diffs['df_3']:
ax.axvline(x=i/60, alpha=0.1)
plt.show()
In [80]:
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
signal_data['df_4']['best_splits'] = get_best_split_points_and_plot_BICs(axes, 'df_4')
plt.show()
In [81]:
best_individual_splitpoints = [
0, #D
0, #P
0, #W
0, #J
0, #R
0, #MV
0, #WF
0] #F
fig, axes = plt.subplots(4,2, figsize=(15,16))
axes = axes.flatten()
events = plot_best_splits(axes, best_individual_splitpoints, 'df_4')
plt.show()
In [82]:
focussed_events = events[-5::-1]
def generate_colour():
while True:
yield tuple(np.random.uniform(0,1,3))
flat = []
[[flat.append(i) for i in e] for e in focussed_events]
flat = list(np.sort(flat))
difference = np.diff(np.sort(flat))
c_map = generate_colour()
colours = []
for i, d in enumerate(difference):
if i == 0:
colours.append(next(c_map))
if d <= 20:
colours.append(colours[i])
else:
colours.append(next(c_map))
colours
fig, ax = plt.subplots(1,1, figsize=(15,4))
i = 0
for event_log in focussed_events:
ax.axhline(y=i, c='black', lw=0.5, alpha=0.5)
for e in event_log[:-1]:
ax.axvline(x=e/60, ymin=(i+1)/(len(focussed_events)), ymax=(i)/(len(focussed_events)), c=colours[flat.index(e)], lw=5)
i += 1
ax.set_ylim(0,i)
ax.grid(False)
ax.set_yticks(np.arange(0,len(water_columns[:4]))+0.5)
ax.set_yticklabels(water_columns[:4][::-1], position=(0,-1))
for i in data_diffs['df_4']:
ax.axvline(x=i/60, alpha=0.1)
plt.show()