Understanding particle trajectories near islands along Haro Strait and along the thalweg¶
References¶
1. Imports¶
In [1]:
from __future__ import division
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
from salishsea_tools import (nc_tools,viz_tools)
from IPython.core.display import Image
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
from salishsea_tools import tidetools
%matplotlib inline
2. Data¶
September 18-27, 2003
Salinity used in thalweg¶
In [2]:
tracers = nc.Dataset('/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/18sep27sep/SalishSea_1d_20030918_20030927_grid_T.nc')
sal = tracers.variables['vosaline']
zlevels = tracers.variables['deptht']
ssal = sal.shape
ssal
Out[2]:
(10, 40, 898, 398)
Vertical velocity entered into Ariane namelist¶
In [3]:
tracersW = nc.Dataset('/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/18sep27sep/SalishSea_1d_20030918_20030927_grid_W.nc')
nc_tools.show_variables(tracers)
velW=tracersW.variables['vovecrtz']
svelW = velW.shape
svelW
[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']
Out[3]:
(10, 40, 898, 398)
Grid for bathymetry plot¶
In [4]:
grid = nc.Dataset('/ocean/imachuca/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
nc_tools.show_variables(grid)
lats = grid.variables['nav_lat']
lons = grid.variables['nav_lon']
bathy = grid.variables['Bathymetry']
bathy.shape
bath, X, Y = tidetools.get_bathy_data(grid)
[u'nav_lon', u'nav_lat', u'Bathymetry']
Thalweg data¶
In [5]:
!head thalweg.txt
thalweg = np.loadtxt('thalweg.txt', dtype=int, unpack=True)
405 2 404 3 404 4 404 5 404 6 404 7 403 8 403 9 402 10 401 11
Ariane Namelist¶
Unconfined particles¶
Namelist input concerning data read: nmax = 7, tunit = 86400 (one day), ntfic = 1 (each time sample covers one day), lmt = 10 (time steps in input data = 10 days).
I want to output trajectory points every 8 hours so I could ideally plot 3 points per day = 30 points total (read below).
delta_t = 3600 (hour unit), frequency = 8 (successive positions are eight hours apart), nb_out = 28 (total number of points). nb_out is 28 instead of 30 in order to satisfy Conditions 1 and 2. Refer to http://salishsea-meopar-docs.readthedocs.org/en/latest/particles/ariane.html
No interpolation¶
delta_t=86400, frequency=1, nb_out=9
Particle Trajectories¶
Along thalweg: unconfined to depths
In [6]:
PT_thal = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/thalweg_islands/along_thalweg/unconfined/ariane_trajectories_qualitative.nc','r')
lonT=PT_thal.variables['traj_lon']
latT=PT_thal.variables['traj_lat']
depT=PT_thal.variables['traj_depth']
xT=PT_thal.variables['init_x']
yT=PT_thal.variables['init_y']
tT=PT_thal.variables['traj_time']
nc_tools.show_variables(PT_thal)
[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']
Islands: unconfined to depths
In [7]:
PT_isl = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/thalweg_islands/islands/unconfined/ariane_trajectories_qualitative.nc','r')
lonI=PT_isl.variables['traj_lon']
latI=PT_isl.variables['traj_lat']
depI=PT_isl.variables['traj_depth']
xI=PT_isl.variables['init_x']
yI=PT_isl.variables['init_y']
tI=PT_thal.variables['traj_time']
Along thalweg: no interpolation
In [8]:
PTN_thal = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/thalweg_islands/along_thalweg/nointerp/ariane_trajectories_qualitative.nc','r')
lonTN=PTN_thal.variables['traj_lon']
latTN=PTN_thal.variables['traj_lat']
depTN=PTN_thal.variables['traj_depth']
xTN=PTN_thal.variables['init_x']
yTN=PTN_thal.variables['init_y']
tTN=PTN_thal.variables['traj_time']
nc_tools.show_variables(PTN_thal), latTN.shape
[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']
Out[8]:
(None, (10, 7))
Islands: no interpolation
In [9]:
PTN_isl = nc.Dataset('/ocean/imachuca/MEOPAR/Ariane/results/thalweg_islands/islands/nointerp/ariane_trajectories_qualitative.nc','r')
lonIN=PTN_isl.variables['traj_lon']
latIN=PTN_isl.variables['traj_lat']
depIN=PTN_isl.variables['traj_depth']
xIN=PTN_isl.variables['init_x']
yIN=PTN_isl.variables['init_y']
tIN=PTN_thal.variables['traj_time']
nc_tools.show_variables(PTN_isl), latIN.shape
[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']
Out[9]:
(None, (10, 7))
3. Thalweg, Salinity, Bathymetry¶
Images from notebooks in Reference section above.
In [10]:
Image(filename='ThalSal.png')
Out[10]:
In [11]:
Image(filename='ThalDepth.png')
Out[11]:
In [12]:
Image(filename='ThalMap.png')
Out[12]:
4. Initial Positions¶
In [13]:
# Set up the figure and axes
fig, axr= plt.subplots(1,1, figsize=(30, 15))
land_colour = 'indigo'
axr.set_axis_bgcolor(land_colour)
axr.set_position((0.83, 0.125, 0.2, 0.775))
# Plot thalweg points on bathymetry map
viz_tools.set_aspect(axr)
cmap = plt.get_cmap('winter_r')
cmap.set_bad(land_colour)
bathy = grid.variables['Bathymetry']
x_slice = np.arange(bathy.shape[1])
y_slice = np.arange(200, 800)
axr.pcolormesh(x_slice, y_slice, bathy[y_slice, x_slice], cmap=cmap)
axr.plot(thalweg[1], thalweg[0], linestyle='-', marker='+', color='wheat', label='Thalweg Points')
legend = axr.legend(loc='best', fancybox=True, framealpha=0.25)
axr.set_xlabel('x Index')
axr.set_ylabel('y Index')
axr.grid()
axr.set_xlim([200,355])
axr.set_ylim([240,450])
xk = [268,272,285,300,289,262,239]
yk = [397,385,365,350,342,345,330]
xw = [320,293,266,324,292,265,302]
yw = [350,331,317,316,295,287,274]
lab = np.arange(1,8)
for xxk, yyk, xxw, yyw, labb in zip (xk,yk,xw,yw,lab):
axr.scatter(xxk,yyk,color='black',marker='s')
axr.scatter(xxw,yyw,color='white',marker='s')
axr.annotate(labb,(xxk,yyk),fontsize=15,color='black')
axr.annotate(labb,(xxw,yyw),fontsize=15,color='white')
5. Depths at initial positions¶
These will be the initial depths for the particles. We are choosing the particles to be at half of the total depth at their initial positions.
In [14]:
bathyk = np.zeros(7)
bathyw = np.zeros(7)
for i in range (7):
bathyk[i] = 0.5*bathy[yk[i],xk[i]]
bathyw[i] = 0.5*bathy[yw[i],xw[i]]
bathyk,bathyw
Out[14]:
(array([ 98.5, 90.5, 111. , 110.5, 102. , 134. , 109. ]), array([ 43.5, 75. , 92. , 48.5, 13. , 34.5, 29. ]))
Ariane takes grid indices not depths. I'm choosing the grid to be that whose corresponding depth is closest to the depths calculated above.
In [15]:
zlevels[:]
Out[15]:
array([ 0.5000003 , 1.5000031 , 2.50001144, 3.50003052,
4.50007057, 5.50015068, 6.50031042, 7.50062323,
8.50123596, 9.50243282, 10.50476551, 11.50931168,
12.51816654, 13.53541183, 14.56898212, 15.63428783,
16.76117325, 18.00713539, 19.48178482, 21.38997841,
24.10025597, 28.22991562, 34.68575668, 44.51772308,
58.48433304, 76.58558655, 98.06295776, 121.86651611,
147.08946228, 173.11448669, 199.57304382, 226.26029968,
253.06663513, 279.93453979, 306.834198 , 333.75018311,
360.67453003, 387.60321045, 414.53408813, 441.46609497], dtype=float32)
Bathyk ~ 98 98 121 121 98 147 98 ~ 27 27 28 28 27 29 27
Bathyw ~ 44 76 98 44 13 34 28 ~ 24 26 27 24 14 23 22
6. Bathymetry: Index-Coords Conversion¶
In [16]:
iyTN = np.zeros(shape=(10,7))
ixTN = np.zeros(shape=(10,7))
iyIN = np.zeros(shape=(10,7))
ixIN = np.zeros(shape=(10,7))
depthTNS=np.empty_like(iyTN)
depthINS=np.empty_like(iyIN)
partcc = np.arange(7)
dyss = np.arange(10)
for partc in zip(partcc):
for dys in zip(dyss):
iyTN,ixTN= tidetools.find_closest_model_point(lonTN[dys,partc], latTN[dys,partc], X, Y, bath, allow_land=False)
iyIN, ixIN = tidetools.find_closest_model_point(round(lonIN[dys,partc],4), round(latIN[dys,partc],4), X, Y, bath, allow_land=False)
depthTNS[dys,partc]=bathy[int(float(iyTN)),int(float(ixTN))]
depthINS[dys,partc]=bathy[int(float(iyIN)),int(float(ixIN))]
7. Along Thalweg¶
In [17]:
(fig,axs)=plt.subplots(7,2,figsize=(15,70))
n = np.arange(7)
for ax,N in zip(axs[:,0],n):
aspect = viz_tools.set_aspect(ax, coords='map', lats=lats)
cmap = plt.get_cmap('Blues')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(lons[:], lats[:], bathy[:], cmap=cmap)
ax.set_xlim([-123.8,-122.7])
ax.set_ylim([48,49.6])
if (ax==axs[6,0]):
ax.set_xlim([-124.6,-122.7])
ax.set_ylim([48,49.6])
ax.scatter(lonT[0,N],latT[0,N],color='red',marker='s', label='Initial Position')
cm = plt.cm.get_cmap('RdYlBu')
ax.scatter(lonT[1:,N],latT[1:,N],c=tT, cmap=cm,marker='o', label='Interp. Trajectory', alpha=0.7)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('Trajectory for particle along thalweg')
ax.legend(scatterpoints = 1, loc='lower right')
for ax,N in zip(axs[:,1],n):
cm = plt.cm.get_cmap('RdYlBu')
ax.plot(tT[:,N]*10, depT[:,N]*-1, '-.', label = 'Interpolated', color='grey')
ax.plot(tTN[:,N]*10, depTN[:,N]*-1, '-o', label = 'Not Interpolated', color='green')
ax.plot(tTN[:,N]*10,depthTNS[:,N],'-o', color='blue', label='Bathymetry')
ax.set_ylim(400, 0)
ax.set_xticks(np.arange(0, 10, 1))
ax.set_xticklabels(np.arange(0, 10, 1))
ax.set_xlabel('Time [days]')
ax.set_ylabel('Depth [m]')
ax.set_title('Depth of particle along trajectory')
ax.legend(scatterpoints = 1, loc='lower left')
8. Around Islands¶
In [18]:
(fig,axs)=plt.subplots(7,2,figsize=(15,70))
n = np.arange(7)
for ax,N in zip(axs[:,0],n):
aspect = viz_tools.set_aspect(ax, coords='map', lats=lats)
cmap = plt.get_cmap('Blues')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(lons[:], lats[:], bathy[:], cmap=cmap)
ax.set_xlim([-123.5,-122.5])
ax.set_ylim([48.2,49])
ax.scatter(lonIN[0,N],latIN[0,N],color='red',marker='s', label='Initial Position')
cm = plt.cm.get_cmap('RdYlBu')
ax.scatter(lonI[1:,N],latI[1:,N],c=tI, cmap=cm,marker='o', label='Interp. Trajectory', alpha=0.7)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('Trajectory for particle around islands')
ax.legend(scatterpoints = 1, loc='lower left')
for ax,N in zip(axs[:,1],n):
cm = plt.cm.get_cmap('RdYlBu')
ax.plot(tI[:,N]*10, depI[:,N]*-1, '-.', label = 'Interpolated', color='grey')
ax.plot(tIN[:,N]*10, depIN[:,N]*-1, '-o', label = 'Not Interpolated', color='green')
ax.plot(tIN[:,N]*10,depthINS[:,N],'-o', color='blue', label='Bathymetry')
ax.set_ylim(400, 0)
ax.set_xticks(np.arange(0, 10, 1))
ax.set_xticklabels(np.arange(0, 10, 1))
ax.set_xlabel('Time [days]')
ax.set_ylabel('Depth [m]')
ax.set_title('Depth of particle along trajectory')
ax.legend(scatterpoints = 1, loc='lower left')
In [18]: