This notebook will be used to show some functions in ellipse.py and maximize the accessibility of this module.
The velocities have to be prepared and set up in a certain way before they can be put through the tidal analysis. They need to be orthogonal, rotated to align with N/S, unstaggered and masked. The functions in ellipse.py handle these calculations.
This notebook outlines how to calculate tidal ellipses for :
import datetime
import matplotlib.pylab as plt
import numpy as np
import netCDF4 as nc
from matplotlib.patches import Ellipse
from salishsea_tools import (viz_tools, tidetools as tt, ellipse, nc_tools)
from salishsea_tools.nowcast import (analyze, research_VENUS)
%matplotlib inline
The functions were originally developed to analyze tidal currents from nowcast results. Given a date range and grid coordinate, a single function can be used to analyze the time series and generate ellipse parameters.
By default, this function will calculate tidal ellipses at each model depth. An example is shown below.
First, define the date range and grid coordinates of interest.
path = '/data/dlatorne/MEOPAR/SalishSea/nowcast/'
#The longer the date range the longer the analysis takes but the more tidal constituents it can resolve
to=datetime.datetime(2015,7,10)
tf=datetime.datetime(2015,7,17)
i = [280] # x-index
j = [400] # y-index
nconst = 2 # number of constituents to be analyzed
Note: The i, j coordinates must be represented as a list or array. If multiple grid points are to be analyzed, they must cover an area. This means that the values in i (and j) must be separated by one. Further, since analyzing currents requires unstaggering, both i and j must be greater than 0.
We also have to define the number of constituents to be analyzed. The possibilities are:
Note that a long time series is required to analyze many constituents.
Next, a simlple call to get_params_nowcast() returns the ellipse parameters in a dictionary object (params) and the depths associated with the ellipse parameters (dep).
params, dep = ellipse.get_params_nowcast(to, tf, i, j, path, nconst)
A quick demonstration of plotting the semi-major axis over depth is shown below.
for key in params.iterkeys():
plt.plot(np.squeeze(params[key]['Semi-Major Axis']), dep[:],label=key)
plt.axis([0,0.7,150,0])
plt.legend(loc=0)
plt.xlabel('Semi-Major Axis (m/s)')
plt.ylabel('Depth (m)')
<matplotlib.text.Text at 0x7ff816be7b50>
If more constituents were analyzed then params would have additional keys representing those constituents.
Each constituent in params also contains keys for the Semi-Minor Axis, Phase, and Inclination.
We can also analyze the tidal currents over an area and at a single depth.
i = np.arange(280,300) # x-indices
j = np.arange(400,440) # y-indices
nconst = 2 # number of constituents to be analyzed
depth = 1 # depth of analysis in meters
params, dep = ellipse.get_params_nowcast(to, tf, i, j, path, nconst,depthrange=depth)
print dep
0.5
Note: the ellipse is analyzed at the closest model level that is less than depth.
It is straightforward to plot the phase.
grid = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
bathy = grid.variables['Bathymetry']
lons = grid.variables['nav_lon']
lats = grid.variables['nav_lat']
fig,ax=plt.subplots(1,1)
mesh=ax.pcolormesh(lons[j,i],lats[j,i],params['M2']['Phase'])
plt.colorbar(mesh,ax=ax)
viz_tools.plot_coastline(ax,grid,coords='map')
ax.set_xlim([-123.5,-123])
ax.set_ylim([48.5,49.5])
ax.set_title('M2 Phase')
viz_tools.set_aspect(ax)
1.1363636363636362
Finally, we can also calculate the ellipses for the depth averaged currents over a depth range. The depav argument should be set to True.
i = [280] # x-index
j = [400] # y-index
depthrange = [0,100] # depth range in meters.
nconst = 4 # number of constituents to be analyzed
params, dep = ellipse.get_params_nowcast(to, tf, i, j, path, nconst,depthrange=depthrange, depav=True)
print 'S2 ellipse parameters'
print params['S2']
print 'Range of depth averaging:'
print dep
S2 ellipse parameters {'Phase': array([[ 354.9102976]]), 'Semi-Major Axis': array([[ 0.04274396]]), 'Semi-Minor Axis': array([[-0.00943675]]), 'Inclination': array([[ 131.62047723]])} Range of depth averaging: [ 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]
And an example for plotting the S2 ellipse.
fig, ax = plt.subplots(1,1)
major = params['S2']['Semi-Major Axis'][0]
minor = params['S2']['Semi-Minor Axis'][0]
inc = params['S2']['Inclination'][0]
#color is blue if counterclockise (minor<0) and red if clockwise
color = 'b'
if minor >= 0:
color = 'r'
ell = Ellipse(xy=(lons[j,i],lats[j,i]),width = major,
height =minor, angle= inc,
color=color)
ax.add_artist(ell)
viz_tools.plot_coastline(ax,grid,coords='map')
ax.set_xlim([-123.5,-123])
ax.set_ylim([48.6,49.2])
viz_tools.set_aspect(ax)
ax.set_title('S2 ellipse depth averaged over {0:.3}-{1:.3} m'.format(dep[0],dep[-1]))
<matplotlib.text.Text at 0x7ff81693e750>
We might want to analyze other simulations outside of the nowcast system.
For example, we have a time series of U/V at the VENUS east node from a simulation saved at a 15 minute frequency. To analyze this time series we need to
These steps are shown in each of cells below.
Load data:
grid = '/ocean/mdunn/MEOPAR/jun_VENUS_east_gridded.nc'
G = nc.Dataset(grid)
u = G.variables['vozocrtx'][:]
v = G.variables['vomecrty'][:]
dep = G.variables['depthu'][:]
#Set up the time array that corresponds to the times of the u and v arrays in hours
t = G.variables['time_counter'][:]/3600.
Rorate and unstagger velocities:
u_rot,v_rot = ellipse.prepare_vel(u, v, depav=True, depth=dep)
Define nodal corrections:
NodalCorr = {
'K1': {
'ft': 0.763545,
'uvt': 344.740346},
'M2': {
'ft': 1.035390,
'uvt': 346.114490}}
nconst = 2
Calculate ellipse parameters:
params = ellipse.get_params(u_rot, v_rot, t, nconst, tidecorr = NodalCorr)
print params['K1']
{'Phase': array([[ 284.96515857]]), 'Semi-Major Axis': array([[ 0.34321001]]), 'Semi-Minor Axis': array([[-0.03599698]]), 'Inclination': array([[ 128.21677008]])}
The above result is a depth averaged analysis. Like earlier, we also have the flexibilty to analyze a range of depths or a single depth by slicing u and v at the desired model levels before passing to prepare_vel().
If you would like to analyze your own tidal currents you need to consider the following:
$u_{sub} = U[..., 399:401,299:301]$ and $v_{sub} = V[..., 399:401,299:301]$.
The result is that $u_{sub}$ and $v_{sub}$ are both of shape (..., 2,2). If you subset this way, then $u_{sub}$ and $v_{sub}$ can be passed into prepare_vel() to perform the unstaggering and rotation.
In these calculations, the phase is defined as the time that the current is aligned in the inclination direction with respect to GMT. The inclination is the direction of the maximum currents. Since this happens twice in an ellipse, the convention is to take the inclination to be between 0 and 180 degrees. As such, if the inclination in a region is close to 0 or 180, then the nearby ellipses may abruptly change inclination and it may appear as though the phase has a large discontinuity. The phase should be interpretted with caution.
While we have shown some plotting basics in this notebook, there are many other functions useful for plotting in the notebooks listed below:
At the moment, the harmonic analysis can only be performed for up to 8 constituents and in predetermined pairings, but, it may be necessary to expand this functionality. Some recommendations are