In [1]:
from matplotlib.path import Path
%matplotlib inline
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import matplotlib.patches as patches
from salishsea_tools import viz_tools, geo_tools, tidetools
from bathy_helpers import *
In [2]:
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')

547 particles were released at various depthts at time=0 in the following domain marked by the black rectangle.

In [41]:
fig, ax = plt.subplots(1,1, figsize=(16,12))
ax.pcolormesh(grid.variables['Bathymetry'][:])
ax.set_ylim(650,800)
ax.set_xlim(100,150)
viz_tools.set_aspect(ax)
ax.add_patch(patches.Rectangle((120,735), 10, 15, fill=False, linewidth=3))
ax.grid('on')
In [39]:
fig, ax = plt.subplots(1,1, figsize=(16,12))
ax.pcolormesh(grid.variables['Bathymetry'][:])
#ax.set_ylim(650,800)
#ax.set_xlim(100,150)
viz_tools.set_aspect(ax)
ax.add_patch(patches.Rectangle((45,750), 90, 898-750, fill=False, linewidth=3))
ax.add_patch(patches.Rectangle((45+90,735), 180, 898-735, fill=False, linewidth=3))
ax.add_patch(patches.Rectangle((100,200), 298, 535, fill=False, linewidth=3))
ax.grid('on')
#plt.savefig('Boxes.png', bbox_inches='tight')
In [43]:
fig, ax = plt.subplots(1,1, figsize=(16,12))
ax.pcolormesh(grid.variables['Bathymetry'][:])
ax.set_ylim(590,670)
ax.set_xlim(115,135)
viz_tools.set_aspect(ax)
ax.add_patch(patches.Rectangle((118,598), 15, 59, fill=False, linewidth=3))
ax.add_patch(patches.Rectangle((120,598), 13, 10, fill=False, linewidth=3))
Out[43]:
<matplotlib.patches.Rectangle at 0x7f8da834ca58>
In [40]:
fig,ax= plt.subplots(1,1)
l = 743
ax.plot(Jul_lont[l,:], Jul_latt[l,:], 'r.')
viz_tools.plot_land_mask(ax, grid, coords='map')
#plt.savefig('Final_destination.png', bbox_inches='tight')
In [24]:
fig,ax= plt.subplots(1,1)
l = 743
ax.plot(Oct_lont[l,:], Oct_latt[l,:], 'r.')
viz_tools.plot_land_mask(ax, grid, coords='map')
Out[24]:
<matplotlib.contour.QuadContourSet at 0x7f8da850a9e8>
In [25]:
fig,ax= plt.subplots(1,1)
l = 743
ax.plot(Jan_lont[l,:], Jan_latt[l,:], 'r.')
viz_tools.plot_land_mask(ax, grid, coords='map')
Out[25]:
<matplotlib.contour.QuadContourSet at 0x7f8da844eba8>
In [27]:
fig,ax= plt.subplots(1,1)
l = 719
ax.plot(Apr_lont[l,:], Apr_latt[l,:], 'r.')
viz_tools.plot_land_mask(ax, grid, coords='map')
Out[27]:
<matplotlib.contour.QuadContourSet at 0x7f8da8614550>
In [70]:
def still_inside(il, ir, jl, jr, lont, latt):
    p = makebox(glamfe,gphife,il,ir,jl,jr)
    poly = path.Path(p, closed=True) 
    time, particles = lont.shape
    number_of_particles = np.zeros(time)
    for l in range(time):
        cx = lont[l,:]
        cy = latt[l,:]
        pts = np.array([cx,cy]).T
        test = poly.contains_points(pts)
        number_of_particles[l]= sum(test) / particles
    return number_of_particles
def deep_particles(lont,latt,init_z, il, ir, jl, jr):
    mask = lont[:].mask
    p = makebox(glamfe,gphife,il,ir,jl,jr)
    poly = path.Path(p, closed=True) 
    time, particles = lont.shape
    number_of_particles = np.zeros(time)
    index_deep_particles=[]
    for n in range(particles):
        if init_z[n] > 6:
            index_deep_particles.append(n)
    length_of_deep_particles = len(index_deep_particles)
    number_of_deep_particles = np.zeros(time)
    for n in range(time):
        for m in index_deep_particles:
            if (mask[n,m]) == False: 
                tf = poly.contains_point(np.array((lont[n,m], latt[n,m])).T)
                if tf == True:
                    number_of_deep_particles[n] = number_of_deep_particles[n] + 1
    return (number_of_deep_particles / length_of_deep_particles)
def shallow_particles(lont,latt,init_z, il, ir, jl, jr):
    mask = lont[:].mask
    p = makebox(glamfe,gphife,il,ir,jl,jr)
    poly = path.Path(p, closed=True) 
    time, particles = lont.shape
    number_of_particles = np.zeros(time)
    index_shallow_particles=[]
    for n in range(particles):
        if init_z[n] < 6:
            index_shallow_particles.append(n)
    length_of_shallow_particles = len(index_shallow_particles)
    number_of_shallow_particles = np.zeros(time)
    for n in range(time):
        for m in index_shallow_particles:
            if (mask[n,m]) == False: 
                tf = poly.contains_point(np.array((lont[n,m], latt[n,m])).T)
                if tf == True:
                    number_of_shallow_particles[n] = number_of_shallow_particles[n] + 1
    return (number_of_shallow_particles / length_of_shallow_particles)

Starting Location

In [98]:
fig, ax = plt.subplots(1,1, figsize=(16,12))
ax.pcolormesh(grid.variables['Bathymetry'][:])
ax.set_ylim(650,800)
ax.set_xlim(100,150)
viz_tools.set_aspect(ax)
ax.add_patch(patches.Rectangle((120,735), 10, 15, fill=False, linewidth=3))
ax.grid('on')
In [92]:
Jul_CR_all = still_inside(120, 130, 735, 750, Jul_lont, Jul_latt)
Oct_CR_all = still_inside(120, 130, 735, 750, Oct_lont, Oct_latt)
Jan_CR_all = still_inside(120, 130, 735, 750, Jan_lont, Jan_latt)
Apr_CR_all = still_inside(120, 130, 735, 750, Apr_lont, Apr_latt)
In [93]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_CR_all, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_CR_all, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_CR_all, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_CR_all, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[93]:
<matplotlib.text.Text at 0x7f8d9b166048>
In [94]:
Jul_CR_deep = deep_particles(Jul_lont, Jul_latt, Jul.variables['init_z'], 120, 130, 735, 750)
Oct_CR_deep = deep_particles(Oct_lont, Oct_latt, Oct.variables['init_z'], 120, 130, 735, 750)
Jan_CR_deep = deep_particles(Jan_lont, Jan_latt, Jan.variables['init_z'], 120, 130, 735, 750)
Apr_CR_deep = deep_particles(Apr_lont, Apr_latt, Apr.variables['init_z'], 120, 130, 735, 750)

Deep Particles

In [95]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_CR_deep, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_CR_deep, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_CR_deep, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_CR_deep, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[95]:
<matplotlib.text.Text at 0x7f8d9aecb208>
In [96]:
Jul_CR_shallow = shallow_particles(Jul_lont, Jul_latt, Jul.variables['init_z'], 120, 130, 735, 750)
Oct_CR_shallow = shallow_particles(Oct_lont, Oct_latt, Oct.variables['init_z'], 120, 130, 735, 750)
Jan_CR_shallow = shallow_particles(Jan_lont, Jan_latt, Jan.variables['init_z'], 120, 130, 735, 750)
Apr_CR_shallow = shallow_particles(Apr_lont, Apr_latt, Apr.variables['init_z'], 120, 130, 735, 750)

Shallow Particles

In [97]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_CR_shallow, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_CR_shallow, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_CR_shallow, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_CR_shallow, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[97]:
<matplotlib.text.Text at 0x7f8d9ace1198>

Northwest

In [45]:
fig, ax = plt.subplots(1,1, figsize=(16,12))
ax.pcolormesh(grid.variables['Bathymetry'][:])
viz_tools.set_aspect(ax)
ax.add_patch(patches.Rectangle((45,750), 90, 898-750, fill=False, linewidth=3))
Out[45]:
<matplotlib.patches.Rectangle at 0x7f8da06aa710>
In [58]:
Jul_NW_all = still_inside(45, 135, 750, 897, Jul_lont, Jul_latt)
Oct_NW_all = still_inside(45, 135, 750, 897, Oct_lont, Oct_latt)
Jan_NW_all = still_inside(45, 135, 750, 897, Jan_lont, Jan_latt)
Apr_NW_all = still_inside(45, 135, 750, 897, Apr_lont, Apr_latt)
In [59]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_NW_all, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_NW_all, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_NW_all, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_NW_all, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[59]:
<matplotlib.text.Text at 0x7f8da1979240>
In [71]:
Jul_NW_deep = deep_particles(Jul_lont, Jul_latt, Jul.variables['init_z'], 45, 135, 750, 897)
Jul_NW_deep[:5]
Out[71]:
array([ 0.       ,  0.7881137,  1.       ,  1.       ,  1.       ])
In [72]:
Oct_NW_deep = deep_particles(Oct_lont, Oct_latt, Oct.variables['init_z'], 45, 135, 750, 897)
Jan_NW_deep = deep_particles(Jan_lont, Jan_latt, Jan.variables['init_z'], 45, 135, 750, 897)
Apr_NW_deep = deep_particles(Apr_lont, Apr_latt, Apr.variables['init_z'], 45, 135, 750, 897)

Deep Particles

In [73]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_NW_deep, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_NW_deep, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_NW_deep, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_NW_deep, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[73]:
<matplotlib.text.Text at 0x7f8da0cda8d0>
In [74]:
Jul_NW_shallow = shallow_particles(Jul_lont, Jul_latt, Jul.variables['init_z'], 45, 135, 750, 897)
Oct_NW_shallow = shallow_particles(Oct_lont, Oct_latt, Oct.variables['init_z'], 45, 135, 750, 897)
Jan_NW_shallow = shallow_particles(Jan_lont, Jan_latt, Jan.variables['init_z'], 45, 135, 750, 897)
Apr_NW_shallow = shallow_particles(Apr_lont, Apr_latt, Apr.variables['init_z'], 45, 135, 750, 897)

Shallow particles

In [75]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_NW_shallow, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_NW_shallow, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_NW_shallow, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_NW_shallow, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[75]:
<matplotlib.text.Text at 0x7f8da0a66198>

Northeast

In [99]:
fig, ax = plt.subplots(1,1, figsize=(16,12))
ax.pcolormesh(grid.variables['Bathymetry'][:])
viz_tools.set_aspect(ax)
ax.add_patch(patches.Rectangle((45+90,735), 180, 898-735, fill=False, linewidth=3))
ax.grid('on')
In [60]:
Jul_NE_all = still_inside(135, 315, 735, 897, Jul_lont, Jul_latt)
Oct_NE_all = still_inside(135, 315, 735, 897, Oct_lont, Oct_latt)
Jan_NE_all = still_inside(135, 315, 735, 897, Jan_lont, Jan_latt)
Apr_NE_all = still_inside(135, 315, 735, 897, Apr_lont, Apr_latt)
In [61]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_NE_all, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_NE_all, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_NE_all, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_NE_all, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[61]:
<matplotlib.text.Text at 0x7f8da16f62b0>
In [76]:
Jul_NE_deep = deep_particles(Jul_lont, Jul_latt, Jul.variables['init_z'], 135, 315, 735, 897)
Oct_NE_deep = deep_particles(Oct_lont, Oct_latt, Oct.variables['init_z'], 135, 315, 735, 897)
Jan_NE_deep = deep_particles(Jan_lont, Jan_latt, Jan.variables['init_z'], 135, 315, 735, 897)
Apr_NE_deep = deep_particles(Apr_lont, Apr_latt, Apr.variables['init_z'], 135, 315, 735, 897)

Deep particles

In [77]:
fig,ax = plt.subplots(2,2,figsize=(12,12))
ax[0,0].plot(range(744), Jul_NE_deep, 'r.')
ax[0,0].grid('on')
ax[0,0].set_title('July 2016')
ax[0,1].plot(range(744), Oct_NE_deep, 'r.')
ax[0,1].grid('on')
ax[0,1].set_title('October 2016')
ax[1,0].plot(range(744), Jan_NE_deep, 'r.')
ax[1,0].grid('on')
ax[1,0].set_title('January 2017')
ax[1,1].plot(range(720), Apr_NE_deep, 'r.')
ax[1,1].grid('on')
ax[1,1].set_title('April 2017')
Out[77]:
<matplotlib.text.Text at 0x7f8da04f4390>