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
%matplotlib inline
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))
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>''')
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')
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)
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)
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)
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)
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')
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)
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)
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')
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')
f0['Beaching_Time'][:,:].sum()
7373
f1['Beaching_Time'][:,:].sum()
7032
f2['Beaching_Time'][:,:].sum()
7033
f3['Beaching_Time'][:,:].sum()
6968
control['Beaching_Time'][:,:].sum()
6924
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
problessthanone = (f0lessthanone+f1lessthanone+f2lessthanone+f3lessthanone+controllessthanone)/5
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
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
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)
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')
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')
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
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')
n=0
for item in header:
print(str(n), item)
n=n+1
0 Seconds 1 YY 2 MM 3 DD 4 hh 5 mm 6 ss 7 MassOil 8 VolOilBeached 9 VolumeBeached 10 VolumeOil 11 Volume 12 Area 13 TeoricalArea 14 Thickness 15 MEvaporated 16 VEvaporated 17 FMEvaporated 18 MDispersed 19 VDispersed 20 FMDispersed 21 MSedimented 22 VSedimented 23 FMSedimented 24 MDissolved 25 VDissolved 26 FMDissolved 27 MChemDisp 28 VChemDisp 29 FMChemDisp 30 MOilRecovered 31 VOilRecovered 32 FMOilRecovered 33 MWaterContent 34 VWaterContent 35 Density 36 Viscosity 37 MBio 38 VBio 39 FMBio 40 CharacteristicDiameter 41 P_Star 42 AnalyteMass1 43 AnalyteMass2 44 AnalyteMass3 45 AnalyteMass4 46 AnalyteMass5 47 AnalyteBio1 48 AnalyteBio2 49 AnalyteBio3 50 AnalyteBio4 51 AnalyteBio5
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]);
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]);
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]);
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]);
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)
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)
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')
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)))
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');
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')
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)
/home/vdo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:2: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
filelist = glob.glob("/ocean/vdo/MIDOSS/spatial_test/L*")
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)
/home/vdo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). This is separate from the ipykernel package so we can avoid doing imports until
total = np.zeros((896, 396))
for file in filelist:
f = nc.Dataset(file)
total = total + np.sum(f['OilConcentration_3D'][1,:,:], axis=0)
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)
<matplotlib.colorbar.Colorbar at 0x7ff6f9503550>
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)
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)