#!/usr/bin/env python
# coding: utf-8

# In[100]:


import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
import scipy.io as sio
import pandas as pd
import datetime
import xarray as xr
from salishsea_tools import tidetools, geo_tools, viz_tools
from matplotlib.colors import LinearSegmentedColormap
import os
import h5py
import glob
get_ipython().run_line_magic('matplotlib', 'inline')

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))


# In[41]:


from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>

<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')


# In[ ]:





# # In one SalishSeaCast grid

# In[16]:


f0 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-0_InOneSSGrid-0.nc')
f1 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-1_InOneSSGrid-1.nc')
f2 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-2_InOneSSGrid-2.nc')
f3 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-3_InOneSSGrid-3.nc')
f4 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-4_InOneSSGrid-4.nc')


# In[20]:


n = 1
fig, ax = plt.subplots(1,5, figsize = (15,5))
ax[0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
ax[1].pcolormesh(np.sum(f1['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
ax[2].pcolormesh(np.sum(f2['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
ax[3].pcolormesh(np.sum(f3['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
ax[4].pcolormesh(np.sum(f4['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
for a in ax:
    a.set_ylim(320,350)
    a.set_xlim(240,260)


# ## surface

# ## top left - bottom right

# In[4]:


for n in np.arange(1,21):
    fig, ax = plt.subplots(1,3, figsize = (15,5))
    ax[0].set_title('h = ' +str(n))
    ax[0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,250:360, 200:270], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1].pcolormesh(np.sum(f4['OilConcentration_3D'][n,:,250:360, 200:270], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    mesh = ax[2].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,250:360, 200:270], axis=0) - 
                         np.sum(f4['OilConcentration_3D'][n,:,250:360, 200:270], axis=0),
                         cmap = 'bwr', vmin = -0.3, vmax = 0.3)
    fig.colorbar(mesh, ax=ax)


# ## top right - bottom left

# In[5]:


for n in np.arange(1,21):
    fig, ax = plt.subplots(1,3, figsize = (15,5))
    ax[0].set_title('h = ' +str(n))
    ax[0].pcolormesh(np.sum(f1['OilConcentration_3D'][n,:,250:360, 200:270], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1].pcolormesh(np.sum(f3['OilConcentration_3D'][n,:,250:360, 200:270], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    mesh = ax[2].pcolormesh(np.sum(f1['OilConcentration_3D'][n,:,250:360, 200:270], axis=0) - 
                         np.sum(f3['OilConcentration_3D'][n,:,250:360, 200:270], axis=0),
                         cmap = 'bwr', vmin = -0.3, vmax = 0.3)
    fig.colorbar(mesh, ax=ax)


# ## top left - center

# In[6]:


for n in np.arange(1,21):
    fig, ax = plt.subplots(1,3, figsize = (15,5))
    ax[0].set_title('h = ' +str(n))
    ax[0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,250:360, 200:270], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1].pcolormesh(np.sum(f2['OilConcentration_3D'][n,:,250:360, 200:270], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    mesh = ax[2].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,250:360, 200:270], axis=0) - 
                         np.sum(f2['OilConcentration_3D'][n,:,250:360, 200:270], axis=0),
                         cmap = 'bwr', vmin = -0.3, vmax = 0.3)
    fig.colorbar(mesh, ax=ax)


# # In one Geotiff grid

# In[21]:


f0 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-0_InOneGeotiffGrid-0.nc')
f1 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-1_InOneGeotiffGrid-1.nc')
f2 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-2_InOneGeotiffGrid-2.nc')


# In[22]:


n = 1
fig, ax = plt.subplots(1,3, figsize = (15,5))
ax[0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
ax[1].pcolormesh(np.sum(f1['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
ax[2].pcolormesh(np.sum(f2['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
for a in ax:
    a.set_ylim(320,350)
    a.set_xlim(240,260)


# In[26]:


for n in np.arange(1,21):
    fig, ax = plt.subplots(1,3, figsize = (15,5))
    ax[0].set_title('h = ' +str(n))
    ax[0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1].pcolormesh(np.sum(f1['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    mesh = ax[2].pcolormesh(np.sum(f2['OilConcentration_3D'][n,:,:], axis=0),
                         cmap = 'Blues_r', vmin = 0, vmax = 0.5)
    fig.colorbar(mesh, ax=ax)
    for a in ax:
        a.set_ylim(280,360)
        a.set_xlim(200,270)


# # Beaching tests

# In[119]:


f0 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-0_beachingtests-0.nc')
f1 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-1_beachingtests-1.nc')
control = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude_control_beaching.nc')


# In[120]:


f2 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-2_beachingtests-2.nc')
f3 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-3_beachingtests-3.nc')


# # control: prob=0.5, limit = 5
# # f0: prob=0.5, limit = 250
# # f1: prob = 1, limit = 250
# # f2: prob = 0.5, limit = 500
# # f3: prob = 1, limit = 500

# ## total beaching hours for f0

# In[95]:


f0['Beaching_Time'][:,:].sum()


# ## total beaching hours for f1

# In[96]:


f1['Beaching_Time'][:,:].sum()


# ## total beaching hours for f2

# In[97]:


f2['Beaching_Time'][:,:].sum()


# ## total beaching hours for f3

# In[98]:


f3['Beaching_Time'][:,:].sum()


# ## total beaching hours for control

# In[99]:


control['Beaching_Time'][:,:].sum()


# In[ ]:





# In[ ]:





# In[299]:


f0lessthanone = f0['Beaching_Time'][:,:]
f0lessthanone[(f0lessthanone <= 24) & (f0lessthanone > 0)] = 1
f0lessthanone[(f0lessthanone > 24) | (f0lessthanone == 0)] = 0

f1lessthanone = f1['Beaching_Time'][:,:]
f1lessthanone[(f1lessthanone <= 24) & (f1lessthanone > 0)] = 1
f1lessthanone[(f1lessthanone > 24) | (f1lessthanone == 0)] = 0

f2lessthanone = f2['Beaching_Time'][:,:]
f2lessthanone[(f2lessthanone <= 24) & (f2lessthanone > 0)] = 1
f2lessthanone[(f2lessthanone > 24) | (f2lessthanone == 0)] = 0

f3lessthanone = f0['Beaching_Time'][:,:]
f3lessthanone[(f0lessthanone <= 24) & (f3lessthanone > 0)] = 1
f3lessthanone[(f0lessthanone > 24) | (f3lessthanone == 0)] = 0

controllessthanone = control['Beaching_Time'][:,:]
controllessthanone[(controllessthanone <= 24) & (controllessthanone > 0)] = 1
controllessthanone[(controllessthanone > 24) | (controllessthanone == 0)] = 0


# In[300]:


problessthanone = (f0lessthanone+f1lessthanone+f2lessthanone+f3lessthanone+controllessthanone)/5


# In[306]:


f0lessthanthree = f0['Beaching_Time'][:,:]
f0lessthanthree[(f0lessthanthree <= 72) & (f0lessthanthree > 0)] = 1
f0lessthanthree[(f0lessthanthree > 72) | (f0lessthanthree == 0)] = 0

f1lessthanthree = f1['Beaching_Time'][:,:]
f1lessthanthree[(f1lessthanthree <= 72) & (f1lessthanthree > 0)] = 1
f1lessthanthree[(f1lessthanthree > 72) | (f1lessthanthree == 0)] = 0

f2lessthanthree = f2['Beaching_Time'][:,:]
f2lessthanthree[(f2lessthanthree <= 72) & (f2lessthanthree > 0)] = 1
f2lessthanthree[(f2lessthanthree > 72) | (f2lessthanthree == 0)] = 0

f3lessthanthree = f3['Beaching_Time'][:,:]
f3lessthanthree[(f3lessthanthree <= 72) & (f3lessthanthree > 0)] = 1
f3lessthanthree[(f3lessthanthree > 72) | (f3lessthanthree == 0)] = 0

controllessthanthree = control['Beaching_Time'][:,:]
controllessthanthree[(controllessthanthree <= 72) & (controllessthanthree > 0)] = 1
controllessthanthree[(controllessthanthree > 72) | (controllessthanthree == 0)] = 0

problessthanthree = (f0lessthanthree+f1lessthanthree+f2lessthanthree+f3lessthanthree+controllessthanthree)/5


# In[307]:


f0lessthanseven = f0['Beaching_Time'][:,:]
f0lessthanseven[(f0lessthanseven <= 168) & (f0lessthanseven > 0)] = 1
f0lessthanseven[(f0lessthanseven > 168) | (f0lessthanseven == 0)] = 0

f1lessthanseven = f1['Beaching_Time'][:,:]
f1lessthanseven[(f1lessthanseven <= 168) & (f1lessthanseven > 0)] = 1
f1lessthanseven[(f1lessthanseven > 168) | (f1lessthanseven == 0)] = 0

f2lessthanseven = f2['Beaching_Time'][:,:]
f2lessthanseven[(f2lessthanseven <= 168) & (f2lessthanseven > 0)] = 1
f2lessthanseven[(f2lessthanseven > 168) | (f2lessthanseven == 0)] = 0

f3lessthanseven = f3['Beaching_Time'][:,:]
f3lessthanseven[(f3lessthanseven <= 168) & (f3lessthanseven > 0)] = 1
f3lessthanseven[(f3lessthanseven > 168) | (f3lessthanseven == 0)] = 0

controllessthanseven = control['Beaching_Time'][:,:]
controllessthanseven[(controllessthanseven <= 168) & (controllessthanseven > 0)] = 1
controllessthanseven[(controllessthanseven > 168) | (controllessthanseven == 0)] = 0

problessthanseven = (f0lessthanseven+f1lessthanseven+f2lessthanseven+f3lessthanseven+controllessthanseven)/5


# # Probability maps

# In[315]:


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

ax[0].pcolormesh(np.ma.masked_equal(problessthanone,0), vmin = 0, vmax = 1)
ax[1].pcolormesh(np.ma.masked_equal(problessthanthree,0), vmin = 0, vmax = 1)
mesh = ax[2].pcolormesh(np.ma.masked_equal(problessthanseven,0), vmin = 0, vmax = 1)
fig.colorbar(mesh, ax=ax)

ax[0].set_title('< 1 day')
ax[1].set_title('< 3 day')
ax[2].set_title('< 7 day')
for a in ax:
    a.set_ylim(280,360)
    a.set_xlim(120,280)


# In[ ]:





# In[ ]:





# In[ ]:





# In[324]:


bins = np.arange(0,8)
fig, ax = plt.subplots(2,3, figsize = (15,10))
ax[0,0].set_title("beaching time, prob = 0.5, limit = 250m")
ax[0,1].set_title("prob = 1, lim = 250m")
ax[0,2].set_title("control, prob = 0.5, lim = 5m")
ax[1,0].set_title("beaching time, prob = 0.5, limit = 500m")
ax[1,1].set_title("prob = 1, lim = 500m")
ax[0,0].hist(f0['Beaching_Time'][:,:].flatten()[f0['Beaching_Time'][:,:].flatten() !=0]/24, bins=bins)
ax[0,1].hist(f1['Beaching_Time'][:,:].flatten()[f1['Beaching_Time'][:,:].flatten() !=0]/24, bins=bins)
ax[0,2].hist(control['Beaching_Time'][:,:].flatten()[control['Beaching_Time'][:,:].flatten() !=0]/24, 
             bins=bins)
ax[1,0].hist(f2['Beaching_Time'][:,:].flatten()[f2['Beaching_Time'][:,:].flatten() !=0]/24, bins=bins)
ax[1,1].hist(f3['Beaching_Time'][:,:].flatten()[f3['Beaching_Time'][:,:].flatten() !=0]/24, 
             bins=bins)
for a in ax.flatten():
    a.set_ylim(0,60)
    a.set_xlabel('days')


# In[327]:


bins = np.arange(0,8)
fig, ax = plt.subplots(2,3, figsize = (15,10))
ax[0,0].set_title("beaching time, prob = 0.5, limit = 250m")
ax[0,1].set_title("prob = 1, lim = 250m")
ax[0,2].set_title("control, prob = 0.5, lim = 5m")
ax[1,0].set_title("beaching time, prob = 0.5, limit = 500m")
ax[1,1].set_title("prob = 1, lim = 500m")
ax[0,0].hist(f0['Beaching_Time'][:,:].flatten()[f0['Beaching_Time'][:,:].flatten() !=0]/24, 
             bins=bins, cumulative = True)
ax[0,1].hist(f1['Beaching_Time'][:,:].flatten()[f1['Beaching_Time'][:,:].flatten() !=0]/24, 
             bins=bins, cumulative = True)
ax[0,2].hist(control['Beaching_Time'][:,:].flatten()[control['Beaching_Time'][:,:].flatten() !=0]/24, 
             bins=bins, cumulative = True)
ax[1,0].hist(f2['Beaching_Time'][:,:].flatten()[f2['Beaching_Time'][:,:].flatten() !=0]/24, 
             bins=bins, cumulative = True)
ax[1,1].hist(f3['Beaching_Time'][:,:].flatten()[f3['Beaching_Time'][:,:].flatten() !=0]/24, 
             bins=bins, cumulative = True)
for a in ax.flatten():
    a.set_ylim(0,160)
    a.set_xlabel('days')


# In[147]:


def load_sro(filepath):
    with open(filepath, 'r') as the_file:
        all_data = [line.strip() for line in the_file.readlines()]
        header = all_data[4]  

    # Order header into list array by splitting up string
    header_arr = []
    header_arr = header.split(' ')
    # Remove emtpy entries from list
    header_arr = np.asarray([x for x in header_arr if x != ''])
    
    data2D = np.genfromtxt(filepath, skip_header=6, skip_footer=4)
    
    return header_arr, data2D


# In[151]:


header, f0data = load_sro('/ocean/vdo/MIDOSS/results/beaching-0.sro')
header, f1data = load_sro('/ocean/vdo/MIDOSS/results/beaching-1.sro')
header, f2data = load_sro('/ocean/vdo/MIDOSS/results/beaching-2.sro')
header, f3data = load_sro('/ocean/vdo/MIDOSS/results/beaching-3.sro')
header, controldata = load_sro('/ocean/vdo/MIDOSS/results/beaching-control.sro')


# In[150]:


n=0
for item in header:
    print(str(n), item)
    n=n+1


# In[164]:


fig, ax = plt.subplots(figsize = (8,8))
ax.plot(f0data[::24,8], label = 'f0')
ax.plot(f1data[::24,8], label = 'f1')
ax.plot(f2data[::24,8], label = 'f2')
ax.plot(f3data[::24,8], label = 'f3')
ax.plot(controldata[::24,8], label = 'control')
ax.legend()
ax.set_xlabel('days')
ax.set_ylabel(header[8]);


# In[320]:


fig, ax = plt.subplots(figsize = (8,8))
ax.plot(f0data[::24,9], label = 'f0')
ax.plot(f1data[::24,9], label = 'f1')
ax.plot(f2data[::24,9], label = 'f2')
ax.plot(f3data[::24,9], label = 'f3')
ax.plot(controldata[::24,9], label = 'control')
ax.legend()
ax.set_xlabel('days')
ax.set_ylabel(header[9]);


# In[322]:


fig, ax = plt.subplots(figsize = (8,8))
ax.plot(f0data[:,8] + f0data[:,34], label = 'VolOilBeached + VWaterContent')
ax.plot(f0data[:,9], label = 'VolumeBeached')
ax.legend()
ax.set_xlabel('days');
#ax.set_ylabel(header[9]);


# In[323]:


fig, ax = plt.subplots(figsize = (8,8))
ax.plot(f0data[:,34] + f0data[:,10], label = 'VolumeOil + VWaterContent')
ax.plot(f0data[:,11], label = 'Volume')
ax.legend()
ax.set_xlabel('days');
#ax.set_ylabel(header[9]);


# In[ ]:





# In[ ]:





# In[ ]:





# In[319]:


fig, ax = plt.subplots(2,3, figsize = (15,10))
ax[0,0].set_title("beaching time, prob = 0.5, limit = 250m")
ax[0,1].set_title("prob = 1, lim = 250m")
ax[0,2].set_title("control, prob = 0.5, lim = 5m")
ax[1,0].set_title("beaching time, prob = 0.5, limit = 500m")
ax[1,1].set_title("prob = 1, lim = 500m")
ax[0,0].pcolormesh(np.ma.masked_equal(f0['Beaching_Time'][:,:], 0),
                 vmin=0, vmax=160)
ax[0,1].pcolormesh(np.ma.masked_equal(f1['Beaching_Time'][:,:], 0),
                 vmin=0, vmax=160)
ax[0,2].pcolormesh(np.ma.masked_equal(control['Beaching_Time'][:,:], 0),
                 vmin=0, vmax=160)
ax[1,0].pcolormesh(np.ma.masked_equal(f2['Beaching_Time'][:,:], 0),
                 vmin=0, vmax=160)
mesh = ax[1,1].pcolormesh(np.ma.masked_equal(f3['Beaching_Time'][:,:], 0),
                 vmin=0, vmax=160)
fig.colorbar(mesh, ax=ax)
for a in ax.flatten():
    a.set_ylim(280,370)
    a.set_xlim(120,270)


# In[94]:


fig, ax = plt.subplots(2,3, figsize = (15,10))
ax[0,0].set_title("f0-f1")
ax[0,1].set_title("f0-f2")
ax[0,2].set_title("f0-f3")
ax[1,0].set_title("f1-f2")
ax[1,1].set_title("f1-f3")
ax[1,2].set_title("f2-f3")
ax[0,0].pcolormesh(f0['Beaching_Time'][:,:] - 
                     f1['Beaching_Time'][:,:],
                     cmap = 'seismic', vmin = -100, vmax = 100)
ax[0,1].pcolormesh(f0['Beaching_Time'][:,:] - 
                     f2['Beaching_Time'][:,:],
                     cmap = 'seismic', vmin = -100, vmax = 100)
ax[0,2].pcolormesh(f0['Beaching_Time'][:,:] - 
                     f3['Beaching_Time'][:,:],
                     cmap = 'seismic', vmin = -100, vmax = 100)
ax[1,0].pcolormesh(f1['Beaching_Time'][:,:] - 
                     f2['Beaching_Time'][:,:],
                     cmap = 'seismic', vmin = -100, vmax = 100)
ax[1,1].pcolormesh(f1['Beaching_Time'][:,:] - 
                     f3['Beaching_Time'][:,:],
                     cmap = 'seismic', vmin = -100, vmax = 100)
mesh = ax[1,2].pcolormesh(f2['Beaching_Time'][:,:] - 
                     f3['Beaching_Time'][:,:],
                     cmap = 'seismic', vmin = -100, vmax = 100)
fig.colorbar(mesh, ax=ax)
for a in ax.flatten():
    a.set_ylim(280,360)
    a.set_xlim(120,270)


# In[222]:


t = h5py.File('/results2/MIDOSS/forcing/SalishSeaCast/MF0/21nov17-28nov17/t.hdf5', 'r')
winds = h5py.File('/results2/MIDOSS/forcing/SalishSeaCast/MF0/21nov17-28nov17/winds.hdf5', 'r')
currents = h5py.File('/results2/MIDOSS/forcing/SalishSeaCast/MF0/21nov17-28nov17/currents.hdf5', 'r')
waves = h5py.File('/results2/MIDOSS/forcing/SalishSeaCast/MF0/21nov17-28nov17/waves.hdf5', 'r')


# In[237]:


f0tem = np.array([])
f0sal = np.array([])
f0wcc = np.array([])
f0ssh = np.array([])
f0wnd = np.array([])
f0cur = np.array([])

f1tem = np.array([])
f1sal = np.array([])
f1wcc = np.array([])
f1ssh = np.array([])
f1wnd = np.array([])
f1cur = np.array([])

f2tem = np.array([])
f2sal = np.array([])
f2wcc = np.array([])
f2ssh = np.array([])
f2wnd = np.array([])
f2cur = np.array([])

f3tem = np.array([])
f3sal = np.array([])
f3wcc = np.array([])
f3ssh = np.array([])
f3wnd = np.array([])
f3cur = np.array([])

control_tem = np.array([])
control_sal = np.array([])
control_wcc = np.array([])
control_ssh = np.array([])
control_wnd = np.array([])
control_cur = np.array([])

for hr in range(1,168):
    
    cur_speed = np.sqrt(currents['Results']['velocity U']['velocity U_00'+ "{0:03}".format(hr)][-1,:]**2 + 
                        currents['Results']['velocity V']['velocity V_00'+ "{0:03}".format(hr)][-1,:]**2)
    wnd_speed = np.sqrt(winds['Results']['wind velocity X']['wind velocity X_00'+ "{0:03}".format(hr)][:]**2 + 
                        winds['Results']['wind velocity Y']['wind velocity Y_00'+ "{0:03}".format(hr)][:]**2)
    
    
    
    f0mask = np.ma.masked_equal(f0['OilConcentration_3D'][hr,-1,:].T, 0).mask
    f0tem = np.append(f0tem, 
                    np.ma.mean(np.ma.masked_array(t['Results']['temperature']['temperature_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f0mask)))
    f0sal = np.append(f0sal, 
                    np.ma.mean(np.ma.masked_array(t['Results']['salinity']['salinity_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f0mask)))
    f0ssh = np.append(f0ssh, 
                    np.ma.mean(np.ma.masked_array(t['Results']['water level']['water level_00' + "{0:03}".format(hr)][:], 
                                                  mask = f0mask)))
    f0wcc = np.append(f0wcc, 
                    np.ma.mean(np.ma.masked_array(waves['Results']['whitecap coverage']['whitecap coverage_00' + "{0:03}".format(hr)][:], 
                                                  mask = f0mask)))
    f0cur = np.append(f0cur, np.ma.mean(np.ma.masked_array(cur_speed, mask=f0mask)))
    f0wnd = np.append(f0wnd, np.ma.mean(np.ma.masked_array(wnd_speed, mask=f0mask)))
    
    
    
    
    f1mask = np.ma.masked_equal(f1['OilConcentration_3D'][hr,-1,:].T, 0).mask
    f1tem = np.append(f1tem, 
                    np.ma.mean(np.ma.masked_array(t['Results']['temperature']['temperature_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f1mask)))
    f1sal = np.append(f1sal, 
                    np.ma.mean(np.ma.masked_array(t['Results']['salinity']['salinity_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f1mask)))
    f1ssh = np.append(f1ssh, 
                    np.ma.mean(np.ma.masked_array(t['Results']['water level']['water level_00' + "{0:03}".format(hr)][:], 
                                                  mask = f1mask)))
    f1wcc = np.append(f1wcc, 
                    np.ma.mean(np.ma.masked_array(waves['Results']['whitecap coverage']['whitecap coverage_00' + "{0:03}".format(hr)][:], 
                                                  mask = f1mask)))
    f1cur = np.append(f1cur, np.ma.mean(np.ma.masked_array(cur_speed, mask=f1mask)))
    f1wnd = np.append(f1wnd, np.ma.mean(np.ma.masked_array(wnd_speed, mask=f1mask)))
    
    f2mask = np.ma.masked_equal(f2['OilConcentration_3D'][hr,-1,:].T, 0).mask
    f2tem = np.append(f2tem, 
                    np.ma.mean(np.ma.masked_array(t['Results']['temperature']['temperature_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f2mask)))
    f2sal = np.append(f2sal, 
                    np.ma.mean(np.ma.masked_array(t['Results']['salinity']['salinity_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f2mask)))
    f2ssh = np.append(f2ssh, 
                    np.ma.mean(np.ma.masked_array(t['Results']['water level']['water level_00' + "{0:03}".format(hr)][:], 
                                                  mask = f2mask)))
    f2wcc = np.append(f2wcc, 
                    np.ma.mean(np.ma.masked_array(waves['Results']['whitecap coverage']['whitecap coverage_00' + "{0:03}".format(hr)][:], 
                                                  mask = f2mask)))
    f2cur = np.append(f2cur, np.ma.mean(np.ma.masked_array(cur_speed, mask=f2mask)))
    f2wnd = np.append(f2wnd, np.ma.mean(np.ma.masked_array(wnd_speed, mask=f2mask)))
    
    
    
    f3mask = np.ma.masked_equal(f3['OilConcentration_3D'][hr,-1,:].T, 0).mask
    f3tem = np.append(f3tem, 
                    np.ma.mean(np.ma.masked_array(t['Results']['temperature']['temperature_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f3mask)))
    f3sal = np.append(f3sal, 
                    np.ma.mean(np.ma.masked_array(t['Results']['salinity']['salinity_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = f3mask)))
    f3ssh = np.append(f3ssh, 
                    np.ma.mean(np.ma.masked_array(t['Results']['water level']['water level_00' + "{0:03}".format(hr)][:], 
                                                  mask = f3mask)))
    f3wcc = np.append(f3wcc, 
                    np.ma.mean(np.ma.masked_array(waves['Results']['whitecap coverage']['whitecap coverage_00' + "{0:03}".format(hr)][:], 
                                                  mask = f3mask)))
    f3cur = np.append(f3cur, np.ma.mean(np.ma.masked_array(cur_speed, mask=f3mask)))
    f3wnd = np.append(f3wnd, np.ma.mean(np.ma.masked_array(wnd_speed, mask=f3mask)))
    
    
    
    
    control_mask = np.ma.masked_equal(control['OilConcentration_3D'][hr,-1,:].T, 0).mask
    control_tem = np.append(control_tem, 
                    np.ma.mean(np.ma.masked_array(t['Results']['temperature']['temperature_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = control_mask)))
    control_sal = np.append(control_sal, 
                    np.ma.mean(np.ma.masked_array(t['Results']['salinity']['salinity_00' + "{0:03}".format(hr)][-1,:], 
                                                  mask = control_mask)))
    control_ssh = np.append(control_ssh, 
                    np.ma.mean(np.ma.masked_array(t['Results']['water level']['water level_00' + "{0:03}".format(hr)][:], 
                                                  mask = control_mask)))
    control_wcc = np.append(control_wcc, 
                    np.ma.mean(np.ma.masked_array(waves['Results']['whitecap coverage']['whitecap coverage_00' + "{0:03}".format(hr)][:], 
                                                  mask = control_mask)))
    control_cur = np.append(control_cur, np.ma.mean(np.ma.masked_array(cur_speed, mask=control_mask)))
    control_wnd = np.append(control_wnd, np.ma.mean(np.ma.masked_array(wnd_speed, mask=control_mask)))


# In[244]:


fig, ax = plt.subplots(2,3, figsize = (20,10))
ax[0,0].set_title('Temperature')
ax[0,1].set_title('Salinity')
ax[0,2].set_title('Whitecap coverage')
ax[1,0].set_title('SSH')
ax[1,1].set_title('Wind speed')
ax[1,2].set_title('Current speed')

ax[0,0].plot(f0tem, label = 'f0')
ax[0,0].plot(f1tem,  label = 'f1')
ax[0,0].plot(f2tem,  label = 'f2')
ax[0,0].plot(f3tem,  label = 'f3')
ax[0,0].plot(control_tem,  label = 'control')
ax[0,0].legend()

ax[0,1].plot(f0sal,  label = 'f0')
ax[0,1].plot(f1sal,  label = 'f1')
ax[0,1].plot(f2sal,  label = 'f2')
ax[0,1].plot(f3sal,  label = 'f3')
ax[0,1].plot(control_sal,  label = 'control')

ax[0,2].plot(f0wcc,  label = 'f0')
ax[0,2].plot(f1wcc,  label = 'f1')
ax[0,2].plot(f2wcc,  label = 'f2')
ax[0,2].plot(f3wcc,  label = 'f3')
ax[0,2].plot(control_wcc,  label = 'control')

ax[1,0].plot(f0ssh,  label = 'f0')
ax[1,0].plot(f1ssh,  label = 'f1')
ax[1,0].plot(f2ssh,  label = 'f2')
ax[1,0].plot(f3ssh,  label = 'f3')
ax[1,0].plot(control_ssh,  label = 'control')

ax[1,1].plot(f0wnd,  label = 'f0')
ax[1,1].plot(f1wnd,  label = 'f1')
ax[1,1].plot(f2wnd,  label = 'f2')
ax[1,1].plot(f3wnd, label = 'f3')
ax[1,1].plot(control_wnd,  label = 'control')

ax[1,2].plot(f0cur,  label = 'f0')
ax[1,2].plot(f1cur,  label = 'f1')
ax[1,2].plot(f2cur,  label = 'f2')
ax[1,2].plot(f3cur,  label = 'f3')
ax[1,2].plot(control_cur, label = 'control')

for a in ax.flatten():
    a.set_xlabel('hours');


# # Delayed start

# In[33]:


f0 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-0_delayed_start-0.nc')
f1 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-1_delayed_start-1.nc')
f2 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-2_delayed_start-2.nc')
f3 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-3_delayed_start-3.nc')
f4 = nc.Dataset('/ocean/vdo/MIDOSS/Lagrangian_AKNS_crude-4_delayed_start-4.nc')


# In[40]:


for n in np.arange(1,25):
    fig, ax = plt.subplots(2,3, figsize = (15,10))
    ax[0,0].set_title('h = ' +str(n))
    ax[0,1].set_title('delayed start: +3h' )
    ax[0,2].set_title('delayed start: +6h' )
    ax[1,0].set_title('delayed start: +24h' )
    ax[1,1].set_title('delayed start: +72h' )
    ax[0,0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[0,1].pcolormesh(np.sum(f1['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[0,2].pcolormesh(np.sum(f2['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1,0].pcolormesh(np.sum(f3['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1,1].pcolormesh(np.sum(f4['OilConcentration_3D'][n,:,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    for a in ax.flatten():
        a.set_ylim(280,360)
        a.set_xlim(200,270)


# # 5x5 SalishSeaCast grid

# In[103]:


filelist = glob.glob("/ocean/vdo/MIDOSS/spatial_test/L*")


# In[105]:


for file in filelist:
    f = nc.Dataset(file)
    fig, ax = plt.subplots()
    ax.pcolormesh(np.sum(f['OilConcentration_3D'][1,:,:], axis=0), vmin=0, vmax=0.5, cmap = "Blues_r")
    ax.set_title(file)
    ax.set_ylim(320,360)
    ax.set_xlim(240,270)


# In[107]:


total = np.zeros((896, 396))
for file in filelist:
    f = nc.Dataset(file)
    total = total + np.sum(f['OilConcentration_3D'][1,:,:], axis=0)


# In[112]:


fig, ax = plt.subplots()
mesh = ax.pcolormesh(total, cmap = "Greens_r")
ax.set_ylim(335,355)
ax.set_xlim(245,260)
fig.colorbar(mesh, ax=ax)


# In[117]:


f0 = nc.Dataset('/ocean/vdo/MIDOSS/spatial_test/Lagrangian_AKNS_crude-0_AKNS-spatial-0.nc')
f21 = nc.Dataset('/ocean/vdo/MIDOSS/spatial_test/Lagrangian_AKNS_crude-21_AKNS-spatial-21.nc')
for n in np.arange(1,21):
    fig, ax = plt.subplots(1,3, figsize = (15,5))
    ax[0].set_title('hour = ' +str(n))
    ax[0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1].pcolormesh(np.sum(f21['OilConcentration_3D'][n,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    mesh = ax[2].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:], axis=0) - 
                         np.sum(f21['OilConcentration_3D'][n,:,], axis=0),
                         cmap = 'bwr', vmin = -0.3, vmax = 0.3)
    fig.colorbar(mesh, ax=ax)
    for a in ax:
        a.set_xlim(220,300)
        a.set_ylim(300,400)


# In[118]:


f0 = nc.Dataset('/ocean/vdo/MIDOSS/spatial_test/Lagrangian_AKNS_crude-0_AKNS-spatial-0.nc')
f1 = nc.Dataset('/ocean/vdo/MIDOSS/spatial_test/Lagrangian_AKNS_crude-1_AKNS-spatial-1.nc')
for n in np.arange(1,21):
    fig, ax = plt.subplots(1,3, figsize = (15,5))
    ax[0].set_title('hour = ' +str(n))
    ax[0].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    ax[1].pcolormesh(np.sum(f1['OilConcentration_3D'][n,:], axis=0), 
                     vmin=0, vmax=0.5, cmap = "Blues_r")
    mesh = ax[2].pcolormesh(np.sum(f0['OilConcentration_3D'][n,:], axis=0) - 
                         np.sum(f1['OilConcentration_3D'][n,:,], axis=0),
                         cmap = 'bwr', vmin = -0.3, vmax = 0.3)
    fig.colorbar(mesh, ax=ax)
    for a in ax:
        a.set_xlim(220,300)
        a.set_ylim(300,400)


# In[ ]: