Find Tidal Ellipses in CODAR region and VENUS ADCP vertical profiles
# imports
%matplotlib inline
import os
import matplotlib.pylab as plt
import matplotlib.ticker as ticker
import numpy as np
import netCDF4 as NC
from scipy.optimize import curve_fit
from salishsea_tools import tidetools
from salishsea_tools import viz_tools
from matplotlib.patches import Ellipse
# get bathymetric and grid information
bathy, X, Y = tidetools.get_SS_bathy_data()
grid = NC.Dataset('../../NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
grid_lats = grid.variables["nav_lat"][:]
grid_lons = grid.variables["nav_lon"][:]
grid_bathy = grid.variables["Bathymetry"]
The time series data (hourly) is written out in files named CODAR. Because they do not cover the whole region, they are not "rebuilt" and so here we are handling the individual node outputs.
v = np.zeros((480,40,898,398)); u =np.zeros((480,40,898,398))
DIRECTORY = '../../myResults/CBase_FTK1b/'
filenames = [f for f in os.listdir(DIRECTORY) if f.startswith('CODAR')]
for filename in sorted(filenames):
fn = NC.Dataset(os.path.join(DIRECTORY, filename),'r')
# read the velocity components
vvalues = fn.variables["vomecrty"][:]
uvalues = fn.variables["vozocrtx"][:]
# get the latitude and longitudes so we can place this on the bigger grid
lats = fn.variables["nav_lat"][:]
lons = fn.variables["nav_lon"][:]
depths = fn.variables["deptht"][:]
time = fn.variables["time_counter"][:]/3600. # want hours not seconds
x1, y1 = tidetools.find_closest_model_point(lons[0,0],lats[0,0],X,Y,bathy,allow_land=True)
v[:,:,x1:x1+lats.shape[0],y1:y1+lats.shape[1]] = vvalues
u[:,:,x1:x1+lats.shape[0],y1:y1+lats.shape[1]] = uvalues
# Just a quick plot to check things are approximately reasonable
imin = 379; imax = 461; jmin = 238; jmax=321
fig, ax = plt.subplots(1,1)
mesh = ax.pcolormesh(v[0,0,imin:imax,jmin:jmax])
fig.colorbar(mesh)
thesize = v[0,0].shape
print thesize
(898, 398)
# Initialize the Tide Information
# frequencies
M2freq = 28.984106 # degrees per hour
M2freq = M2freq*np.pi/180. # radians per hour
K1freq = 15.041069*np.pi/180.
# initialize fit
vM2amp = np.zeros(thesize); vM2pha = np.zeros(thesize)
vK1amp = np.zeros(thesize); vK1pha = np.zeros(thesize)
uM2amp = np.zeros(thesize); uM2pha = np.zeros(thesize)
uK1amp = np.zeros(thesize); uK1pha = np.zeros(thesize)
mean = np.zeros(thesize)
# function for fit, this assumes two tides only
def double(x, M2amp, M2pha, K1amp, K1pha, mean):
return (mean + M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
K1amp*np.cos(K1freq*x-K1pha*np.pi/180.))
need to unstagger the velocities, before I calculate the components so v is north/south and that's the first index so 2 vus = v[i,j]+v[i-1,j] to put it on the T-grid 2 uus = u[i,j]+u[i,j-1] to put it on the T-grid
print v.shape
ts = 240 # last half of series
length = v.shape[0]-ts
id = 1 # depth
uaus = np.zeros((length,thesize[0],thesize[1]))
vaus = np.zeros((length,thesize[0],thesize[1]))
print vaus.shape
uaus[:,imin:imax,jmin:jmax] = 0.5*(u[ts:,id,imin:imax,jmin+1:jmax+1]+u[ts:,id,imin:imax,jmin:jmax])
vaus[:,imin:imax,jmin:jmax] = 0.5*(v[ts:,id,imin+1:imax+1,jmin:jmax]+v[ts:,id,imin:imax,jmin:jmax])
plt.plot(uaus[:,imin+50,jmin+50],'b',vaus[:,imin+50,jmin+50],'r')
(480, 40, 898, 398) (240, 898, 398)
[<matplotlib.lines.Line2D at 0x3d669d0>, <matplotlib.lines.Line2D at 0x5265a10>]
# function to do the fitting and actually find the tidal components
def fittit (uaus, vaus, time, imin, imax, jmin, jmax, dj=1):
for i in np.arange(imin,imax):
for j in np.arange(jmin,jmax,dj):
if vaus[:,i,j].any() != 0.:
fitted, cov = curve_fit(double,time[ts:],vaus[:,i,j])
if fitted[0] < 0:
fitted[0] = -fitted[0]
fitted[1] = fitted[1]+180.
if fitted[1] > 180:
fitted[1] = fitted[1] - 360.
elif fitted[1] < -180-360:
fitted[1] = fitted[1] + 720.
elif fitted [1] < -180:
fitted[1] = fitted[1] + 360.
if fitted[2] < 0:
fitted[2] = -fitted[2]
fitted[3] = fitted[3]+180.
vM2amp[i,j] = fitted[0]
vM2pha[i,j] = fitted[1]
vK1amp[i,j] = fitted[2]
vK1pha[i,j] = fitted[3]
print 'vaus done'
for i in np.arange(imin,imax):
for j in np.arange(jmin,jmax,dj):
if uaus[:,i,j].any() !=0.:
fitted, cov = curve_fit(double,time[ts:],uaus[:,i,j])
if fitted[0] < 0:
fitted[0] = -fitted[0]
fitted[1] = fitted[1]+180.
if fitted[1] > 180+360:
fitted[1] = fitted[1] -720
elif fitted[1] > 180:
fitted[1] = fitted[1] - 360.
elif fitted[1] < -180-360:
fitted[1] = fitted[1] + 720.
elif fitted [1] < -180:
fitted[1] = fitted[1] + 360.
if fitted[2] < 0:
fitted[2] = -fitted[2]
fitted[3] = fitted[3]+180.
uM2amp[i,j] = fitted[0]
uM2pha[i,j] = fitted[1]
uK1amp[i,j] = fitted[2]
uK1pha[i,j] = fitted[3]
return vM2amp, vM2pha, vK1amp, vK1pha, uM2amp, uM2pha, uK1amp, uK1pha
Find the tidal amplitudes and phases. Print one out to make sure its working.
vM2amp, vM2pha, vK1amp, vK1pha, uM2amp, uM2pha, uK1amp, uK1pha = fittit (uaus, vaus,
time, imin, imax, jmin, jmax)
print uM2amp[420,280],uM2pha[420,280],vM2amp[420,280],vM2pha[420,280]
vaus done 0.146836064382 -178.916482371 0.261093909426 69.0354522705
Plot the amplitude and phase of the two velocities
fig = plt.figure(figsize=(2*5.5+1,2*5))
cmaprot = plt.get_cmap('hsv')
plt.subplot(2,2,1)
plt.pcolormesh(grid_lons[imin:imax,jmin:jmax],grid_lats[imin:imax,jmin:jmax],vM2amp[imin:imax,jmin:jmax])
plt.colorbar()
plt.title('v M2 amp')
plt.ylim(grid_lats[imin,jmin],grid_lats[imax,jmax])
plt.xlim(grid_lons[imax,jmin],grid_lons[imin,jmax])
plt.subplot(2,2,2)
plt.pcolormesh(grid_lons[imin:imax,jmin:jmax],grid_lats[imin:imax,jmin:jmax],vM2pha[imin:imax,jmin:jmax],cmap=cmaprot)
plt.colorbar()
plt.title('v M2 pha')
plt.ylim(grid_lats[imin,jmin],grid_lats[imax,jmax])
plt.xlim(grid_lons[imax,jmin],grid_lons[imin,jmax])
plt.subplot(2,2,3)
plt.pcolormesh(grid_lons[imin:imax,jmin:jmax],grid_lats[imin:imax,jmin:jmax],uM2amp[imin:imax,jmin:jmax])
plt.colorbar()
plt.title('u M2 amp')
plt.ylim(grid_lats[imin,jmin],grid_lats[imax,jmax])
plt.xlim(grid_lons[imax,jmin],grid_lons[imin,jmax])
plt.subplot(2,2,4)
plt.pcolormesh(grid_lons[imin:imax,jmin:jmax],grid_lats[imin:imax,jmin:jmax],uM2pha[imin:imax,jmin:jmax],cmap=cmaprot)
plt.colorbar()
plt.title('u M2 pha')
plt.ylim(grid_lats[imin,jmin],grid_lats[imax,jmax])
plt.xlim(grid_lons[imax,jmin],grid_lons[imin,jmax]);
# this function calculates ellipse parameters based on the amplitude and phase. From Foreman (1977)
def ellipse_params (uM2amp, uM2pha, vM2amp, vM2pha):
CX = uM2amp*np.cos(np.pi*uM2pha/180.)
SX = uM2amp*np.sin(np.pi*uM2pha/180.)
CY = vM2amp*np.cos(np.pi*vM2pha/180.)
SY = vM2amp*np.sin(np.pi*vM2pha/180.)
ap = np.sqrt((CX+SY)**2+(CY-SX)**2)/2. # amplitude of positively rotating component
am = np.sqrt((CX-SY)**2+(CY+SX)**2)/2. # amplitude of negatively rotating component
ep = np.arctan2(CY-SX,CX+SY)+np.pi # phase of positive
em = np.arctan2(CY+SX,CX-SY)+np.pi # phase of negative
major = ap+am # major axis
minor = ap-am # signed minor axis
theta = (ep+em)/2.*180./np.pi # axis tilt
theta %= 180 # by convention theta is between 0 and 180.
return CX, SX, CY, SY, ap, am, ep, em, major, minor, theta
Calculate ellipse parameters and plot angles of them.
CX, SX, CY, SY, ap, am, ep, em, major, minor, theta = ellipse_params (uM2amp, uM2pha, vM2amp, vM2pha)
fig=plt.figure(figsize=(2*5.5+1,2*5))
plt.subplot(221)
plt.pcolormesh(ep[imin:imax,jmin:jmax])
plt.colorbar()
plt.title('ep: phase of pos')
plt.ylim(0,imax-imin)
plt.xlim(0,jmax-jmin)
plt.subplot(222)
plt.pcolormesh(em[imin:imax,jmin:jmax])
plt.colorbar()
plt.title('em: phse of neg')
plt.ylim(0,imax-imin)
plt.xlim(0,jmax-jmin)
plt.subplot(223)
plt.pcolormesh(90-theta[imin:imax,jmin:jmax],cmap=cmaprot)
plt.colorbar()
plt.title('phi: flipped theta')
plt.ylim(0,imax-imin)
plt.xlim(0,jmax-jmin)
plt.subplot(224)
plt.pcolormesh(theta[imin:imax,jmin:jmax], cmap=cmaprot)
plt.colorbar()
plt.title('theta: ellipse tilt')
plt.ylim(0,imax-imin)
plt.xlim(0,jmax-jmin);
# Failed attempt to plot on lat/long.... this one distorts the ellipses, note that the 45 degree
# ellipse is not at 45 degrees.
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
viz_tools.set_aspect(ax, coords='map', lats=grid_lats)
scale = 0.04
interval = 4
for i in np.arange(imin-1,imax,interval):
for j in np.arange(jmin,jmax,interval):
if ap[i,j] > am[i,j]:
thec = 'b'
else:
thec = 'r'
ells = Ellipse(xy=(grid_lons[i,j],grid_lats[i,j]), width=scale*major[i,j],
height=scale*minor[i,j], angle=theta[i,j]+29.,
color=thec)
ax.add_artist(ells)
ells.set_facecolor(thec)
ax.set_ylim(grid_lats[imin,jmin], grid_lats[imax,jmax])
ax.set_xlim(grid_lons[imax,jmin],grid_lons[imin,jmax]);
ell = Ellipse(xy=(grid_lons[390,255],grid_lats[390,255]),width = scale*5, height = scale*0.1, angle = 45,
color='g')
ax.add_artist(ell)
viz_tools.plot_land_mask(ax, '../../NEMO-forcing/grid/bathy_meter_SalishSea2.nc', coords='map')
viz_tools.plot_coastline(ax, '../../NEMO-forcing/grid/bathy_meter_SalishSea2.nc',coords='map', isobath = 5)
ax.set_ylabel('Latitude')
ax.set_xlabel('Longitude')
ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%3.1f'))
ax.set_title("M2 tidal ellipses at 1.5 m depth")
print "Green ellipse = 0.5 m/s, red is clockwise"
Green ellipse = 0.5 m/s, red is clockwise
# rotate grid, then label axes in lats/longs
phi = 29. # the twist in the grid
k = np.zeros((898,398)); m = np.zeros((898,398))
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
viz_tools.set_aspect(ax)
ex = 50
for i in np.arange(imin-ex,imax+ex):
for j in np.arange(jmin-ex,jmax+ex):
k[i,j] = i*np.cos(phi*np.pi/180.)+j*np.sin(phi*np.pi/180.)
m[i,j] = -i*np.sin(phi*np.pi/180.)+j*np.cos(phi*np.pi/180.)
scale = 8.
interval = 4
for i in np.arange(imin-1,imax,interval):
for j in np.arange(jmin,jmax,interval):
if ap[i,j] > am[i,j]:
thec = 'b'
else:
thec = 'r'
ells = Ellipse(xy=(m[i,j],k[i,j]), width=scale*major[i,j],
height=scale*minor[i,j], angle=theta[i,j]+29.,
color=thec)
ax.add_artist(ells)
ells.set_facecolor(thec)
# y-axis in k, but labelled in latitude
ax.set_ylim(445,560)
slope = (grid_lats[imax,jmin]-grid_lats[imin,jmin])/(k[imax,jmin]-k[imin,jmin])
mylist = (k[imin,jmin]+(np.arange(48.8, 49.35, 0.1)-
grid_lats[imin,jmin])/slope).tolist()
labels = ['48.8', '48.9', '49', '49.1', '49.2', '49.3']
ax.set_yticks(mylist); ax.set_yticklabels(labels)
ax.set_ylabel('Latitude (degrees N)')
# x-axis in m, but labelled in longitude
ax.set_xlim(-14,94)
slope = (grid_lons[imin,jmax]-grid_lons[imin,jmin])/(m[imin,jmax]-m[imin,jmin])
mylist = (m[imin,jmin]+(np.arange(-123.7,-123.05,0.1)-
grid_lons[imin,jmin])/slope).tolist()
labels = ['123.7','123.6','123.5','123.4','123.3','123.2','123.1','123']
ax.set_xticks(mylist); ax.set_xticklabels(labels)
ax.set_xlabel('Longitude (degrees W)')
# scale ellipse
ell = Ellipse(xy=(35,465),width = scale*0.5, height = scale*0.1, angle = 45,
color='g')
ax.add_artist(ell)
# land, and 5 m contour
contour_interval = [-0.01, 0.01]
ax.contourf(m[imin-ex:imax+ex,jmin-ex:jmax+ex],k[imin-ex:imax+ex,jmin-ex:jmax+ex],
bathy.data[imin-ex:imax+ex,jmin-ex:jmax+ex],contour_interval,colors='black')
ax.contour(m[imin-ex:imax+ex,jmin-ex:jmax+ex],k[imin-ex:imax+ex,jmin-ex:jmax+ex],
bathy.data[imin-ex:imax+ex,jmin-ex:jmax+ex],[5],colors='black')
ax.set_title("M2 tidal ellipses at 1.5 m depth")
print "Green ellipse = 0.5 m/s, red is clockwise"
Green ellipse = 0.5 m/s, red is clockwise
# figure out transformation for latitude/longitude angles above
plt.subplot(2,2,1)
plt.plot(k[imin:imax,jmin],grid_lats[imin:imax,jmin])
slope = (grid_lats[imax,jmin]-grid_lats[imin,jmin])/(k[imax,jmin]-k[imin,jmin])
plt.plot(k[imin:imax,jmin],grid_lats[imin,jmin]+(k[imin:imax,jmin]-k[imin,jmin])*slope)
plt.subplot(2,2,2)
plt.plot(grid_lats[imin:imax,jmin],k[imin:imax,jmin])
plt.plot(grid_lats[imin:imax,jmin],k[imin,jmin]+(grid_lats[imin:imax,jmin]-grid_lats[imin,jmin])/slope)
plt.subplot(2,2,3)
plt.plot(m[imin,jmin:jmax],grid_lons[imin,jmin:jmax])
slope = (grid_lons[imin,jmax]-grid_lons[imin,jmin])/(m[imin,jmax]-m[imin,jmin])
plt.plot(m[imin,jmin:jmax],grid_lons[imin,jmin]+(m[imin,jmin:jmax]-m[imin,jmin])*slope)
plt.subplot(2,2,4)
plt.plot(grid_lons[imin,jmin:jmax],m[imin,jmin:jmax])
plt.plot(grid_lons[imin,jmin:jmax],m[imin,jmin]+(grid_lons[imin,jmin:jmax]-grid_lons[imin,jmin])/slope)
j=384;i=308
plt.figure(figsize=(5,5))
plt.plot(uaus[:,384,308],vaus[:,384,308])
plt.xlim(-1.4,1.4)
plt.ylim(-1.4,1.4)
print 360-theta[384,308], major[384,308], minor[384,308]
print CX[384,308], CY[384,308], SX[384,308], SY[384,308]
print uM2amp[384,308], uM2pha[384,308], vM2amp[384,308],vM2pha[384,308]
x=np.arange(ts,ts+length)
plt.plot(CX[j,i]*np.cos(M2freq*x)+SX[j,i]*np.sin(M2freq*x),
CY[j,i]*np.cos(M2freq*x)+SY[j,i]*np.sin(M2freq*x));
228.471888748 0.61444943943 0.0403693666508 -0.183195714243 0.259114523203 -0.365109351471 0.381014376364 0.408491748079 -116.645472016 0.460773579029 55.7817532098
j=415; i=285
plt.figure(figsize=(5,5))
plt.plot(uaus[:,415,285],vaus[:,415,285])
plt.xlim(-0.7,0.7)
plt.ylim(-0.7,0.7)
plt.plot(CX[j,i]*np.cos(M2freq*x)+SX[j,i]*np.sin(M2freq*x),
CY[j,i]*np.cos(M2freq*x)+SY[j,i]*np.sin(M2freq*x))
print theta[415,285]
91.9701293886
j=400;i=244
plt.figure(figsize=(5,5))
plt.plot(uaus[:,j,i],vaus[:,j,i])
plt.xlim(-0.7,0.7)
plt.ylim(-0.7,0.7)
plt.plot(CX[j,i]*np.cos(M2freq*x)+SX[j,i]*np.sin(M2freq*x),
CY[j,i]*np.cos(M2freq*x)+SY[j,i]*np.sin(M2freq*x));
ival = 424
jdeep = 267
jsha = 282;
ts = 240 # last half of series
length = v.shape[0]-ts
dr = v.shape[1]
uaus = np.zeros((length,dr,2))
vaus = np.zeros((length,dr,2))
print vaus.shape
uaus = 0.5*(u[ts:,:,ival,(jdeep+1,jsha+1)]+u[ts:,:,ival,(jdeep,jsha)])
vaus = 0.5*(v[ts:,:,ival+1,(jdeep,jsha)]+v[ts:,:,ival,(jdeep,jsha)])
(240, 40, 2)
vM2amp, vM2pha, vK1amp, vK1pha, uM2amp, uM2pha, uK1amp, uK1pha = fittit (uaus, vaus,
time, 0, dr,
0, 2)
CX, SX, CY, SY, ap, am, ep, em, major, minor, theta = ellipse_params (uM2amp, uM2pha, vM2amp, vM2pha)
vaus done
# correct theta (harder)
# in j,i space theta = angle up from j toward i
# in i,j space phi = angle up from i toward j
# but theta is angle from j toward i
fig = plt.figure(figsize=(5.5,6.4))
ax = fig.add_subplot(121)
scale = 7
interval = 2
for j in (0,1):
for i in np.arange(0,dr,interval):
if ap[i,j] > am[i,j]:
thec = 'b'
else:
thec = 'r'
ells = Ellipse(xy=(6*j,i), width=scale*major[i,j],
height=scale*minor[i,j], angle=-theta[i,j]-29,
color=thec)
ax.add_artist(ells)
ells.set_facecolor(thec)
ax.set_ylim(-2,40)
ax.invert_yaxis()
ax.set_xlim(-2,8)
ell = Ellipse(xy=(3,38),width = scale*0.5, height = scale*0.1, angle =-45,
color='g')
ax.add_artist(ell)
ax.set_title("M2 tidal ellipses at two ADCP stations")
ax.set_ylabel("Vertical Grid Point")
print "Green ellipse = 0.5 m/s, red is clockwise"
Green ellipse = 0.5 m/s, red is clockwise
print " Central East"
print "Depth Major-Axis Minor-Axis Angle Major-Axis Minor-Axis Angle "
print " (m) (m/s) (m/s) (deg. (m/s) (m/s) (deg."
print " ccw E) ccw E)"
for i in np.arange(0,dr,interval):
if major[i,1] > 0 :
print '{0:5.1f} {1:.2f} {2:.2f} {3:.0f} {4:.2f} {5:.2f} {6:.0f}'.format(depths[i],
major[i,0], minor[i,0], theta[i,0]+29, major[i,1], minor[i,1], theta[i,1]+29)
elif major[i,0] > 0 :
print '{0:5.1f} {1:.2f} {2:.2f} {3:.0f}'.format(depths[i],
major[i,0], minor[i,0], theta[i,0]+29)
Central East Depth Major-Axis Minor-Axis Angle Major-Axis Minor-Axis Angle (m) (m/s) (m/s) (deg. (m/s) (m/s) (deg. ccw E) ccw E) 0.5 0.24 -0.12 152 0.27 -0.13 129 2.5 0.19 -0.08 150 0.22 -0.08 129 4.5 0.17 -0.07 142 0.18 -0.02 133 6.5 0.19 -0.08 138 0.13 0.03 134 8.5 0.20 -0.09 138 0.11 0.06 99 10.5 0.20 -0.09 139 0.14 0.02 86 12.5 0.20 -0.08 139 0.18 -0.01 93 14.6 0.19 -0.07 139 0.20 -0.04 101 16.8 0.18 -0.06 139 0.22 -0.05 107 19.5 0.17 -0.04 138 0.24 -0.05 111 24.1 0.16 -0.03 135 0.24 -0.05 113 34.7 0.16 -0.03 130 0.22 -0.03 114 58.5 0.17 -0.02 128 0.20 -0.00 122 98.1 0.14 0.00 124 0.19 -0.01 129 147.1 0.12 0.03 128 0.19 0.05 86 199.6 0.13 0.04 134 253.1 0.15 0.03 121 306.8 0.12 0.00 119
j=7*4;i=1
plt.figure(figsize=(5,5))
plt.plot(uaus[:,j,i],vaus[:,j,i])
plt.xlim(-0.7,0.7)
plt.ylim(-0.7,0.7)
plt.plot(CX[j,i]*np.cos(M2freq*x)+SX[j,i]*np.sin(M2freq*x),
CY[j,i]*np.cos(M2freq*x)+SY[j,i]*np.sin(M2freq*x));