In [1]:

import cmocean.cm as cm
import h5py
import matplotlib as mpl
import matplotlib.colorbar as colorbar
import matplotlib.colors as mplcolours
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

%matplotlib inline

In [2]:

mesh = xr.open_dataset('~/MEOPAR/grid/mesh_mask201702.nc')

In [4]:

oil0d = pd.read_csv('/data/sallen/results/MIDOSS/ParticleNoTests/SoG_2000_bunkerc/resOilOutput.sro', sep='\s+', skiprows=4)
oil0d = oil0d.drop([0], axis=0)
length =len(oil0d)
oil0d = oil0d.drop([length-3, length-2, length-1, length], axis=0)

In [5]:

fig, ax = plt.subplots(1, 1, figsize=(15,5))
oil0d.Area.plot(ax=ax);
plt.grid();

In [6]:

fig, ax = plt.subplots(1, 1, figsize=(15,5))
oil0d.MDissolved.plot();
oil0d.MEvaporated.plot();
oil0d.MassOil.plot();
oil0d.MDispersed.plot();
oil0d.MBio.plot()
massbeached = (oil0d.VolOilBeached*oil0d.Density/(1-oil0d.VWaterContent)
    *(1-oil0d.MWaterContent))

massbeached.plot(label="Beached")

plt.plot(oil0d.MDissolved + oil0d.MEvaporated + oil0d.MassOil 
         + massbeached + oil0d.MDispersed);


plt.legend();
plt.grid();

In [7]:

(oil0d.Density/(1-oil0d.VWaterContent)
    *(1-oil0d.MWaterContent)).plot();

In [8]:

oil0d.MDispersed.plot();

In [9]:

fig, ax = plt.subplots(1, 1, figsize=(15, 3.5))
ax.plot(oil0d.MDissolved + oil0d.MEvaporated + oil0d.MassOil 
         + massbeached + oil0d.MDispersed + oil0d.MBio);
ax.grid();

In [10]:

massbeached.plot(label="Beached");
plt.grid()

In [11]:

(oil0d.MDispersed/(oil0d.MDispersed + oil0d.MassOil))[:-2].plot()
plt.grid()

In [12]:

oil0d.VolumeOil.plot();
plt.legend();
plt.grid();

In [13]:

oil0d.AnalyteMass1.plot()
oil0d.AnalyteMass2.plot()
oil0d.AnalyteMass3.plot()
oil0d.AnalyteMass4.plot()
oil0d.AnalyteMass5.plot();

In [14]:

oil0d.Thickness.plot();

In [15]:

oilLag = xr.open_dataset('/data/sallen/results/MIDOSS/ParticleNoTests/SoG_2000_bunkerc/Lagrangian_SoG_2000_bunkerc_SoG_2000_bunkerc.nc')

In [16]:

# 460:480, 240:260
imin, imax = 460, 560
jmin, jmax = 180, 260

In [17]:

fig, axs = plt.subplots(1, 2, figsize=(15, 5))
it = 20

field = oilLag.OilWaterColumnOilVol_3D[it]
(field[:, imin:imax , jmin:jmax].sum(axis=0)).plot(ax=axs[0], cmap='copper')
print (oilLag.OilWaterColumnOilVol_3D[it, :, imin:imax, jmin:jmax].sum(axis=0).sum(axis=0).sum(axis=0))

axs[1].plot(oilLag.time[0:72], oilLag.OilWaterColumnOilVol_3D[0:72].sum(axis=1).sum(axis=1).sum(axis=1)
           )

<xarray.DataArray 'OilWaterColumnOilVol_3D' ()>
array(9998.7025)
Coordinates:
    time     datetime64[ns] 2018-01-15T20:30:00

Out[17]:

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

In [18]:

oil0d.VolumeOil[0:72].plot()

Out[18]:

<matplotlib.axes._subplots.AxesSubplot at 0x7f10d09c4a10>

In [19]:

tmax = 20
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
imin, imax = 465, 490
jmin, jmax = 230, 255
oilLag.OilWaterColumnOilVol_3D[1:tmax+1, 39, imin:imax , jmin:jmax].sum(axis=0).plot(ax=axs[0]);
oilLag.OilWaterColumnOilVol_3D[1:tmax+1, 38, imin:imax , jmin:jmax].sum(axis=0).plot(ax=axs[1]);
oilLag.OilWaterColumnOilVol_3D[1:tmax+1, 0:38, imin:imax , jmin:jmax].sum(axis=0).sum(axis=0).plot(ax=axs[2]);

In [20]:

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
imin, imax = 465, 490
jmin, jmax = 230, 255
oilLag.OilWaterColumnOilVol_3D[:, 39, imin:imax, jmin:jmax].sum(axis=0).plot(ax=ax, cmap='gist_ncar');

In [21]:

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
imin, imax = 460, 560
jmin, jmax = 180, 260
oilLag.OilWaterColumnOilVol_3D[:, 39, imin:imax, jmin:jmax].sum(axis=0).plot(ax=ax, cmap='gist_ncar');

In [22]:

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
imin, imax = 540, 800
jmin, jmax = 100, 220
oilLag.Beaching_Volume[imin:imax , jmin:jmax].plot(ax=ax, cmap='copper')
print (oilLag.Beaching_Volume.sum(axis=0).sum(axis=0))
print (oilLag.Beaching_Volume[imin:imax, jmin:jmax].sum(axis=0).sum(axis=0))

<xarray.DataArray 'Beaching_Volume' ()>
array(10011.67851866)
<xarray.DataArray 'Beaching_Volume' ()>
array(10006.67136475)

In [23]:

bt = np.array(oilLag.Beaching_Time[imin:imax, jmin:jmax] - oilLag.Beaching_Time[imin:imax, jmin:jmax].min())
plt.plot(bt/1e9/3600, 'ro');

In [25]:

watercolour = 'lightskyblue'
landcolour = 'papayawhip'
waterland_cmap = mplcolours.LinearSegmentedColormap.from_list('mycmap', [(0, watercolour),  (0.85, watercolour), 
                                                    (0.850001, landcolour), (1, landcolour)])
xs = range(jmin, jmax)
ys = range(imin, imax)
xx, yy = np.meshgrid(xs, ys)
mymask = oilLag.Beaching_Time[imin:imax, jmin:jmax] == oilLag.Beaching_Time[imin:imax, jmin:jmax].min()
days = np.ma.array(bt/1e9/86400., mask=mymask)

In [26]:

cnorm = mplcolours.Normalize(vmin=0, vmax=7)
mycmap = mpl.cm.copper
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.set_aspect(500/440.)
ax.tick_params(axis='both', which='both', bottom=False, left=False, labelbottom=False, labelleft=False)
ax.pcolormesh(1-mesh.tmask[0, 0, imin:imax, jmin:jmax], cmap=waterland_cmap)
ax.scatter(xx-jmin+1, yy-imin+1, s=10, c=days, cmap='copper', vmin=0, vmax=7, marker='o', edgecolors='face');
ax_cbar = fig.add_axes([0.9, 0.12, 0.05, 0.76])
cb = colorbar.ColorbarBase(ax=ax_cbar, cmap=mycmap, norm=cnorm);
cb.set_label('Days to Beaching')
fig.suptitle('10,000,000 $\ell$ Oil Spilled on Jan 15, 2018 in Strait of Georgia');

In [27]:

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
np.log(oilLag.Beaching_Volume[imin:imax , jmin:jmax]).plot(ax=ax, cmap='copper_r', 
                                                               vmax=10, vmin=0);

/home/sallen/anaconda/envs/py37/lib/python3.7/site-packages/xarray/core/computation.py:604: RuntimeWarning: divide by zero encountered in log
  result_data = func(*input_data)

In [ ]: