%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import netCDF4 as NC
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from salishsea_tools import tidetools
from salishsea_tools import (nc_tools,viz_tools)
import matplotlib.cm as cm
from matplotlib import animation
from numpy import *
from pylab import *
grid = NC.Dataset('/ocean/imachuca/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
bathy, X, Y = tidetools.get_bathy_data(grid)
lats = grid.variables['nav_lat']
lons = grid.variables['nav_lon']
bath = grid.variables['Bathymetry']
#nc_tools.show_variables(grid)
nlats = lats[350:480,250:340]
nlons = lons[350:480,250:340]
nbath = bath[350:480,250:340]
#np.min(nlons)
ngrid=[nlons,nlats,nbath]
To get an idea of what indices to use. An initial position on land might be causing errors in Ariane.
# Bathymetry (2D)
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
viz_tools.set_aspect(ax)
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(bathy, cmap=cmap)
fig.colorbar(mesh)
# Bathymetry (Close-Up)
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
viz_tools.set_aspect(ax)
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(bathy, cmap=cmap)
fig.colorbar(mesh)
plt.axis((300, 340, 390, 430))
# Bathymetry (3D)
fig = plt.figure(figsize=(20,5))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(lons[300:500,200:350],lats[300:500,200:350],-bath[300:500,200:350], cmap=cm.BrBG, alpha='0.75')
ax.view_init(elev=50, azim=-90)
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.set_ylim3d(np.min(lats[300:500,200:350]),np.max(lats[300:500,200:350]))
ax.set_xlim3d(np.min(lons[300:500,200:350]),np.max(lons[300:500,200:350]))
(-124.02245330810547, -122.58061981201172)
2006 December 14 - 15
f = NC.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/dec2006/ariane_trajectories_qualitative.nc','r');
nc_tools.show_variables(f)
[u'init_x', u'init_y', u'init_z', u'init_t', u'init_age', u'init_transp', u'final_x', u'final_y', u'final_z', u'final_t', u'final_age', u'final_transp', u'traj_lon', u'traj_lat', u'traj_depth', u'traj_time']
lont=f.variables['traj_lon']
latt=f.variables['traj_lat']
dept=f.variables['traj_depth']
xs=f.variables['init_x']
ys=f.variables['init_y']
lont.shape
(48, 8)
latt.shape
(48, 8)
Condition 2: delta_t × frequency × nb_output < tunit × ntfic × (lmt + 0.5 - max(fl))
3600x1x47<3600x1x(48.5-1)===47<48.5-1===47<47.5
rd = NC.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/riverparticles/depths/ariane_trajectories_qualitative.nc','r');
rd_lont=rd.variables['traj_lon']
rd_latt=rd.variables['traj_lat']
rd_dept=rd.variables['traj_depth']
rd_xs=rd.variables['init_x']
rd_ys=rd.variables['init_y']
rd_lont.shape
(48, 6)
rt = NC.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/riverparticles/times/ariane_trajectories_qualitative.nc','r');
rt_lont=rt.variables['traj_lon']
rt_latt=rt.variables['traj_lat']
rt_dept=rt.variables['traj_depth']
rt_xs=rt.variables['init_x']
rt_ys=rt.variables['init_y']
rt_time=rt.variables['init_t']
rt_lont.shape
(25, 6)
tracers = NC.Dataset('/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/dec2006/all_forcing/1hour/SalishSea_1h_20061214_20061215_grid_T.nc')
sal = tracers.variables['vosaline']
timesteps = tracers.variables['time_counter']
s_depth = tracers.variables['deptht']
s_lon = tracers.variables['nav_lon']
s_lat = tracers.variables['nav_lat']
nc_tools.show_variables(tracers)
shape0 = sal.shape
shape0
[u'nav_lon', u'nav_lat', u'deptht', u'time_counter', u'time_counter_bnds', u'sossheig', u'votemper', u'vosaline', u'rain_rate', u'snow_rate']
(48, 40, 898, 398)
Tracking multiple particles (7 sharing same x-index and 1 coming out of Fraser). No background, just land mask. All particles are at the surface.
n = np.arange(8)
colors = cm.rainbow(np.linspace(0, 1, len(n)))
# 2D
fig, ax = plt.subplots(1,1,figsize=(10,6.6))
for N,c in zip(n,colors):
ax.scatter(lont[1:,N],latt[1:,N],color=c)
ax.scatter(lont[0,N],latt[0,N],color='0.30',marker='s')
viz_tools.plot_land_mask(ax,grid,coords='map')
ax.set_xlim([-123.8,-122.8])
ax.set_ylim([48.5,49.5])
ax.set_title('Surface particles')
ax.set_xlabel('lon')
ax.set_ylabel('lat')
# 3D
fig = plt.figure(figsize=(20,10))
ax = fig.gca(projection='3d')
for N,c in zip(n,colors):
ax.scatter(lont[1:,N],latt[1:,N],dept[1:,N],color=c)
ax.scatter(lont[0,N],latt[0,N],dept[0,N],color='0.30',marker='s')
ax.set_xlabel('lon')
ax.set_ylabel('lat')
ax.set_zlabel('depth')
ax.set_title('Surface particles')
<matplotlib.text.Text at 0x7f9dd4b4a7d0>
The river particle is one of the particles in the multiple particle simulation (above).
# 2D
fig, ax = plt.subplots(1,1,figsize=(10,6.6))
NN=0
ax.scatter(lont[1:,NN],latt[1:,NN],c='purple')
ax.scatter(lont[0,NN],latt[0,NN],color='0.30',marker='s')
viz_tools.plot_land_mask(ax,grid,coords='map')
ax.set_xlim([-123.5,-122.9])
ax.set_ylim([48.9,49.3])
ax.set_title('River particle at the surface')
ax.set_xlabel('lon')
ax.set_ylabel('lat')
# 3D
fig = plt.figure(figsize=(20,10))
ax = fig.gca(projection='3d')
NN=0
ax.scatter(lont[1:,NN],latt[1:,NN],dept[1:,NN],c='purple')
ax.scatter(lont[0,NN],latt[0,NN],dept[0,NN],color='0.30',marker='s')
ax.set_xlabel('lon')
ax.set_ylabel('lat')
ax.set_zlabel('depth')
ax.set_title('Surface particles')
<matplotlib.text.Text at 0x7f9dd41626d0>
Animation of trajectory of river particle shown above.
# 2D
#Empty map
fig, ax = plt.subplots(1,1,figsize=(10,6.6))
NN=0
viz_tools.plot_land_mask(ax,grid,coords='map')
ax.set_xlim([-123.5,-122.9])
ax.set_ylim([48.9,49.3])
ax.set_title('River particle at the surface')
ax.set_xlabel('lon')
ax.set_ylabel('lat')
#Initial image i.e. initial position only
def init():
ax.scatter(lont[0,NN],latt[0,NN],color='0.30',marker='s')
#Progression of particles - frames max = 47
def animate(p):
ax.scatter(lont[(p+1),NN],latt[(p+1),NN],c='purple')
#The animation function
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=3, interval=300, blit=True, repeat=False)
#A line that makes it all work
mywriter = animation.FFMpegWriter()
#Save in current folder
anim.save('RiverParticle_Animation.mp4',writer=mywriter,fps=1)
#Show as a pop-up window
#plt.show()
# ********************************************************
# 3D
#Empty map
fig = plt.figure(figsize=(20,10))
ax = fig.gca(projection='3d')
ax.set_xlabel('lon')
ax.set_ylabel('lat')
ax.set_zlabel('depth')
ax.set_title('River particle at the surface')
NN=0
ax.set_xlim([-123.5,-122.9])
ax.set_ylim([48.9,49.3])
ax.set_zlim([-0.01,0.01])
#Initial image i.e. initial position only
def init():
ax.scatter(lont[0,NN],latt[0,NN],dept[0,NN],color='0.30',marker='s')
#Progression of particles - frames max = 47
def animate(p):
ax.scatter(lont[(p+1),NN],latt[(p+1),NN],c='purple')
#The animation function
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=3, interval=300, blit=True, repeat=False)
#A line that makes it all work
mywriter = animation.FFMpegWriter()
#Save in current folder
anim.save('RiverParticle3D_Animation.mp4',writer=mywriter,fps=1)
#Show as a pop-up window
#plt.show()
Trajectories differ between particles that have the same x-y coordinates but different depths. In the z direction, the grid is 1 metre - so, particles within the same metre of depth will have the same trajectory (they have the same velocity information). Also, note that the lower particles seem to get stuck before they leave the mouth of the river - this is probably due to a salt wedge - see section 8.
# 3D
n = np.arange(6)
colors = cm.rainbow(np.linspace(0, 1, len(n)))
fig = plt.figure(figsize=(20,10))
ax = fig.gca(projection='3d')
for N,c in zip(n,colors):
ax.scatter(rd_lont[1:,N],rd_latt[1:,N],rd_dept[1:,N],color=c)
ax.scatter(rd_lont[0,N],rd_latt[0,N],rd_dept[0,N],color='0.30',marker='s')
ax.set_xlabel('lon')
ax.set_ylabel('lat')
ax.set_zlabel('depth')
ax.set_title('Particles at varying depths')
ax.view_init(elev=0, azim=-90)
#***********************************************************
# 3D Animation
n = np.arange(6)
colors = cm.rainbow(np.linspace(0, 1, len(n)))
#Empty map
fig = plt.figure(figsize=(20,15))
ax1 = fig.add_subplot(2, 2, 1)
viz_tools.plot_land_mask(ax1,grid,coords='map')
ax1.set_xlim([-123.5,-122.9])
ax1.set_ylim([48.9,49.3])
ax1.set_xlabel('lon')
ax1.set_ylabel('lat')
ax1.set_title('Particles at varying depths')
ax2 = fig.add_subplot(2, 2, 2, projection='3d')
ax2.view_init(elev=0, azim=-90)
ax2.set_xlim([-123.5,-122.9])
ax2.set_ylim([48.9,49.3])
ax2.set_xlabel('lon')
ax2.set_ylabel('lat')
ax2.set_zlabel('depth')
ax2.set_title('Particles at varying depths')
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
ax3.view_init(elev=30, azim=-90)
ax3.set_xlim([-123.5,-122.9])
ax3.set_ylim([48.9,49.3])
ax3.set_xlabel('lon')
ax3.set_ylabel('lat')
ax3.set_zlabel('depth')
ax3.set_title('Particles at varying depths')
ax4 = fig.add_subplot(2, 2, 4, projection='3d')
ax4.view_init(elev=40, azim=45)
ax4.set_xlim([-123.5,-122.9])
ax4.set_ylim([48.9,49.3])
ax4.set_xlabel('lon')
ax4.set_ylabel('lat')
ax4.set_zlabel('depth')
ax4.set_title('Particles at varying depths')
#Initial image i.e. initial position only
def init():
for N,c in zip(n,colors):
ax1.scatter(rd_lont[0,N],rd_latt[0,N],color='0.30',marker='s')
ax2.scatter(rd_lont[0,N],rd_latt[0,N],rd_dept[0,N],color='0.30',marker='s')
ax3.scatter(rd_lont[0,N],rd_latt[0,N],rd_dept[0,N],color='0.30',marker='s')
ax4.scatter(rd_lont[0,N],rd_latt[0,N],rd_dept[0,N],color='0.30',marker='s')
#Progression of particles - frames max = 47
def animate(p):
for N,c in zip(n,colors):
ax1.scatter(rd_lont[(p+1),N],rd_latt[(p+1),N],color=c)
ax2.scatter(rd_lont[(p+1),N],rd_latt[(p+1),N],rd_dept[(p+1),N],color=c)
ax3.scatter(rd_lont[(p+1),N],rd_latt[(p+1),N],rd_dept[(p+1),N],color=c)
ax4.scatter(rd_lont[(p+1),N],rd_latt[(p+1),N],rd_dept[(p+1),N],color=c)
#The animation function
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=47, interval=300, blit=True, repeat=False)
#A line that makes it all work
mywriter = animation.FFMpegWriter()
#Save in current folder
anim.save('RiverParticleDepths_Animation.mp4',writer=mywriter,fps=1)
#Show as a pop-up window
#plt.show()
Now, we'll look at the same river particle we have been working with but we'll see what the trajectories are like when they start at different times. Initial times are 1st, 2nd, 6th, 12th, 18th, 24th hours.
fig, axs = plt.subplots(1, 3, figsize=(16, 8))
part = np.arange(3)
times = rt_time[part]
for ax, p, t in zip(axs, part,times):
ax.scatter(rt_lont[0,p],rt_latt[0,p],color='0.30',marker='s')
ax.scatter(rt_lont[1:,p],rt_latt[1:,p],c='purple')
viz_tools.plot_land_mask(ax,grid,coords='map')
ax.set_xlim([-123.5,-122.9])
ax.set_ylim([48.9,49.3])
ax.set_xlabel('lon')
ax.set_ylabel('lat')
ax.set_title('t = {:.1f}h'.format(t))
fig, axs = plt.subplots(1, 3, figsize=(16, 8))
part = np.arange(3)+3
times = rt_time[part]
for ax, p, t in zip(axs, part,times):
ax.scatter(rt_lont[0,p],rt_latt[0,p],color='0.30',marker='s')
ax.scatter(rt_lont[1:,p],rt_latt[1:,p],c='purple')
viz_tools.plot_land_mask(ax,grid,coords='map')
ax.set_xlim([-123.5,-122.9])
ax.set_ylim([48.9,49.3])
ax.set_xlabel('lon')
ax.set_ylabel('lat')
ax.set_title('River particle at t = {:.1f}h'.format(t))