Look at full-depth CTD casts from Line W transects, estimate vertical scales of eddy over slope. The following cruises are potentially useful.
* 2009 aug-sept
* 2007 apr maybe
* 2006 apr
* 2004 sept
* 2002 oct
* 2001 oct
* 1997 dec - cold streamer perhaps?
* 1997 feb - cold streamer again?
Year 2007 has CTD data as .nc files for each station. Let's look at that first
2007 apr has crossing of streamer + cyclonic eddy.
In [3]:
%matplotlib inline
%load_ext secnum
%secnum
import numpy as np
import scipy as sp
import pandas as pd
import h5py
import netCDF4
import sys
sys.path.append('/media/data/Work/tools/python/dc/')
from dcutils import *
import dcocean
import numexpr as ne
from scipy.io import loadmat
import scipy.interpolate
import scipy.integrate
from datetime import datetime
import time
from IPython import display
import pylab
import seawater as sw
# matplotlib section
import matplotlib as mpl
import matplotlib.pyplot as plt
# big figures
mpl.rcParams['savefig.dpi'] = 2 * mpl.rcParams['savefig.dpi']
# better color cycling
import brewer2mpl
colors = brewer2mpl.get_map('Set2', 'qualitative', 8).mpl_colors
mpl.rcParams['axes.color_cycle'] = colors
import glob
1997 feb, dec and 2001 don't have many cross-sections. Neither does 2002oct
In [4]:
dirname = 'data/linew/'
years = ['2009/', '2007/', '2006/', '2004/', '2002/']
prefix = ['2009sep', '2007apr', '2006apr', '2004sep', '2002oct']
# ring center profile location
centers = [9010, 9008, 9013, 9009, 9009]
# climatological data
woa = loadmat('../datasets/woa05.mat')
# profile locations - access as linew[9001] for e.g.
linew1 = {
# ID lat lon depth
9001: [40.290, -70.200, 90 ],
9002: [40.141, -70.116, 120 ],
9003: [40.000, -70.000, 160 ],
9004: [39.906, -69.930, 700 ],
9005: [39.859, -69.901, 1000],
9006: [39.790, -69.853, 1500],
9007: [39.700, -69.800, 2000],
9008: [39.477, -69.639, 2400],
9009: [39.350, -69.537, 2500],
9010: [39.087, -69.354, 3000],
9011: [38.820, -69.200, 3300],
9012: [38.558, -69.010, 3400],
9013: [38.322, -68.869, 3800],
9014: [38.090, -68.700, 4300],
9015: [37.860, -68.515, 4400],
9016: [37.621, -68.381, 4600],
9017: [37.380, -68.214, 4700],
9018: [37.139, -66.060, 4900],
9019: [36.898, -67.899, 4950],
9020: [36.658, -67.740, 4950],
9021: [36.206, -67.452, 5000],
9022: [35.709, -67.158, 5100],
9023: [35.228, -66.871, 5100],
9024: [34.740, -66.580, 5200],
9025: [34.260, -66.292, 5215],
9026: [33.785, -65.995, 5100]}
# screw the above, use pandas for nice interface to this array
linew = np.array([
[9001, 40.290, -70.200, 90 ],
[9002, 40.141, -70.116, 120 ],
[9003, 40.000, -70.000, 160 ],
[9004, 39.906, -69.930, 700 ],
[9005, 39.859, -69.901, 1000],
[9006, 39.790, -69.853, 1500],
[9007, 39.700, -69.800, 2000],
[9008, 39.477, -69.639, 2400],
[9009, 39.350, -69.537, 2500],
[9010, 39.087, -69.354, 3000],
[9011, 38.820, -69.200, 3300],
[9012, 38.558, -69.010, 3400],
[9013, 38.322, -68.869, 3800],
[9014, 38.090, -68.700, 4300],
[9015, 37.860, -68.515, 4400],
[9016, 37.621, -68.381, 4600],
[9017, 37.380, -68.214, 4700],
[9018, 37.139, -68.060, 4900],
[9019, 36.898, -67.899, 4950],
[9020, 36.658, -67.740, 4950],
[9021, 36.206, -67.452, 5000],
[9022, 35.709, -67.158, 5100],
[9023, 35.228, -66.871, 5100],
[9024, 34.740, -66.580, 5200],
[9025, 34.260, -66.292, 5215],
[9026, 33.785, -65.995, 5100]])
lwdata = pd.DataFrame(linew, index=linew[:,0], columns=('ID', 'lat', 'lon', 'depth'))
In [5]:
trackdist, ang = sw.dist(lwdata.lat, lwdata.lon, units='km')
trackdist = np.insert(np.cumsum(trackdist), 0, 0)
# store profiles, z vector for eddy center locations
centerprof = [] # profiles at eddy center
centerzz = []; # z-vector at center
Hmodes = []; # horizontal velocity mode
Vmodes = []; # vertical velocity mode
cmodes = []; # gravity wave speed
zmodes = []; # z-vector for calculated modes
zstr = []; # stretched z co-ordinate
for ii in range(len(years)):
# list to store anomaly locations
anom = []
rhoctd = []
azz = []
adist = []
#plt.figure(1)
folder = dirname + years[ii]
# iterate over standard locations
for counter, ff in enumerate(lwdata.ID):
# get CTD cast file name
globstr = folder + prefix[ii] + '*' + str(int(ff)) + '*.ctd'
try:
fname = glob.glob(globstr)[0]
except IndexError:
# some stations are missing in some years
# in that case, the [0] indexing does not work
# since fname=[]
continue
# read in CTD data: pressure, temperature, salinity, oxygen
data = np.loadtxt(fname, skiprows=6)
# calculate density for line-W profile
rho = sw.pden(data[:,2], data[:,1], data[:,0])
#plt.figure(1)
#plt.plot(rho, -1*data[:,0])
# get background density from WOA atlas
xi = find_nearest(woa['X'] - 360, lwdata.loc[ff].lon)
yi = find_nearest(woa['Y'], lwdata.loc[ff].lat)
twoa = woa['temp'][:, yi, xi]
swoa = woa['sal'][:, yi, xi]
rhowoa = sw.pden(swoa, twoa, np.squeeze(sw.pres(woa['Z'], woa['Y'][yi])))
#plt.hold(True)
#plt.plot(rhowoa, -1*woa['Z'])
# interpolate woa data to cruise grid
ifunc = sp.interpolate.interp1d(np.squeeze(woa['Z']), rhowoa)
rhowi = ifunc(data[:,0])
# subtract to get anomaly
drho = rho - rhowi
# if center profile, save data
if counter+1 == centers[ii]-9000:
centerprof.append(drho)
centerzz.append(-1*data[:,0])
# calculate vertical mode at center of eddy
# no negative sign since rhowoa[0] = surface
N2 = 9.81/1025 * np.diff(rhowoa)/np.diff(np.squeeze(woa['Z']))
#N2 = 9.81/1025 * np.diff(rhowi)/np.diff(np.squeeze(data[:,0]))
# remove NaNs from N2
inan = np.isnan(N2).argmax() - 1
if inan == -1:
inan = len(N2)
zmode = np.float64(np.squeeze(woa['Z'][0:inan+1]))
N = np.sqrt(N2[0:inan]);
# save stretched co-ordinate
ifuncN = sp.interpolate.interp1d(dcocean.avg1(zmode), N,
bounds_error=False, fill_value=N[-1])
zstretch = 1/np.sqrt(N.min()) * sp.integrate.cumtrapz(ifuncN(data[:,0]), data[:,0])
zstr.append(zstretch)
#zmode = data[0:inan+1,0]
vmode, hmode, cmode = dcocean.vertmode(N2[0:inan], zmode, 3, 0)
Vmodes.append(vmode)
Hmodes.append(hmode)
cmodes.append(cmode)
zmodes.append(dcocean.avg1(zmode))
# save data for gridding later
#anom.append(drho)
rhoctd.append(rho)
azz.append(-1*data[:,0])
adist.append(trackdist[counter] * np.ones_like(data[:,0]))
#axx.append(lwdata.loc[ff].lon * np.ones_like(data[:,0]))
#ayy.append(lwdata.loc[ff].lat * np.ones_like(data[:,0]))
# plot anomaly
#plt.figure(2)
#plt.hold(True)
#plt.plot(drho, -1*data[:,0])
# grid anomaly onto regular grid for each cruise
adistvec = np.hstack(adist)
azzvec = np.hstack(azz)
#xgrd = trackdist
#zgrd = np.linspace(adistvec.min(), adistvec.max(), 60)
xgrd, zgrd = np.meshgrid(np.linspace(adistvec.min(), adistvec.max(), 50),
np.linspace(azzvec.min(), azzvec.max(), 2000))
rhogrid = sp.interpolate.griddata((np.hstack(adist), np.hstack(azz)),
np.hstack(rhoctd), (xgrd, zgrd))
#anomgrid = sp.interpolate.griddata((np.hstack(adist), np.hstack(azz)),
# np.hstack(anom), (xgrd, zgrd))
plotrho = np.ma.masked_array(rhogrid-1000, np.isnan(rhogrid))
#plotanom = np.ma.masked_array(anomgrid, np.isnan(anomgrid))
# plot!
fig = plt.figure()
handle = plt.pcolor(xgrd, zgrd, plotrho, cmap='RdYlBu_r')
ax = plt.gca()
fig.colorbar(handle)
plt.hold(True)
plt.contour(xgrd, zgrd, plotrho,
levels=[24, 25, 26, 26.5, 27, 27.5, 28], colors='k')
plt.ylim((-1000,0))
plt.xlim((0, 800))
ax.set_xlabel('Cross-track dist (km)')
ax.set_ylabel('Pressure (dbar)')
ax2 = ax.twiny()
ax2.set_xticks(trackdist);
ax2.set_xticklabels([int(x) for x in lwdata.ID-9000],
rotation='vertical', size='x-small');
ax2.set_title(prefix[ii])
plt.tight_layout()
In [16]:
# plot center profiles
fig = plt.figure()
fig.hold(True)
ax1 = plt.subplot(1,2,1)
ax2 = plt.subplot(1,2,2)
for ii in range(len(centerprof)):
ax1.plot(centerprof[ii], centerzz[ii])
rhomode = np.diff(Hmodes[ii][:,0])/np.diff(zmodes[ii])
ax2.plot(rhomode, -1*avg1(zmodes[ii]))
ax1.set_title('Profiles of density anomaly at ring center')
ax1.set_ylabel('Pressure (dbar)')
ax1.set_xlabel('Density')
ax1.axvline(0)
ax1.set_ylim([-2000, 0])
ax1.axhline(-700)
ax2.axhline(-700)
ax2.axvline(0)
ax2.set_title('density mode shape')
ax2.set_ylim([-2000, 0])
plt.tight_layout()
ax1.legend(years, loc=3)
# plot center profiles
fig = plt.figure()
fig.hold(True)
ax1 = plt.subplot(1,2,1)
ax2 = plt.subplot(1,2,2)
#ax = fig.gca()
for ii in range(len(centerprof)):
ax1.plot(dcocean.avg1(centerprof[ii]), -1*zstr[ii])
rhomode = np.diff(Hmodes[ii][:,0])/np.diff(zstr[ii])
ax2.plot(rhomode, -1*avg1(zstr[ii]))
ax1.set_title('Density profile - stretched co-ordinate')
ax1.set_ylabel('Pressure (dbar)')
ax1.set_xlabel('Density')
ax1.axvline(0)
#ax1.set_ylim([-2000, 0])
#ax1.axhline(-700)
ax2.axhline(-700)
ax2.axvline(0)
ax2.set_title('density mode shape')
ax2.set_ylim([-2000, 0])
plt.tight_layout()
The ring "centers" in a water depth of
In [20]:
[lwdata.depth[x] for x in centers]
Out[20]:
In [ ]:
In [ ]:
trackdist, ang = sw.dist(lwdata.lat, lwdata.lon, units='km')
trackdist = np.insert(np.cumsum(trackdist), 0, 0)
for ii in range(len(years)):
# list to store anomaly locations
uladcp = []
vladcp = []
azz = []
adist = []
#plt.figure(1)
folder = dirname + years[ii]
# iterate over standard locations
for counter, ff in enumerate(lwdata.ID):
# get CTD cast file name
globstr = folder + prefix[ii] + '*' + str(int(ff)) + '*.ladcp.txt'
try:
fname = glob.glob(globstr)[0]
except IndexError:
# some stations are missing in some years
# in that case, the [0] indexing does not work
# since fname=[]
continue
# read in LADCP data: z, u, v, ev(????)
data = np.loadtxt(fname, skiprows=7)
depth = data[:,0]
u = data[:,1]
v = data[:,2]
# save data for gridding later
#anom.append(drho)
uladcp.append(u)
vladcp.append(v)
azz.append(-1*data[:,0])
adist.append(trackdist[counter] * np.ones_like(data[:,0]))
#axx.append(lwdata.loc[ff].lon * np.ones_like(data[:,0]))
#ayy.append(lwdata.loc[ff].lat * np.ones_like(data[:,0]))
# plot anomaly
#plt.figure(2)
#plt.hold(True)
#plt.plot(drho, -1*data[:,0])
# grid anomaly onto regular grid for each cruise
adistvec = np.hstack(adist)
azzvec = np.hstack(azz)
xgrd, zgrd = np.meshgrid(np.linspace(adistvec.min(), adistvec.max(), 50),
np.linspace(azzvec.min(), azzvec.max(), 2000))
ugrid = sp.interpolate.griddata((np.hstack(adist), np.hstack(azz)),
np.hstack(uladcp), (xgrd, zgrd))
vgrid = sp.interpolate.griddata((np.hstack(adist), np.hstack(azz)),
np.hstack(vladcp), (xgrd, zgrd))
#anomgrid = sp.interpolate.griddata((np.hstack(adist), np.hstack(azz)),
# np.hstack(anom), (xgrd, zgrd))
plotu = np.ma.masked_array(ugrid, np.isnan(ugrid))
plotv = np.ma.masked_array(vgrid, np.isnan(vgrid))
plotvel = plotu**2 + plotv**2
# plot!
fig = plt.figure()
handle = plt.pcolor(xgrd, zgrd, plotvel, cmap='RdYlBu_r')
ax = plt.gca()
fig.colorbar(handle)
#plt.hold(True)
#plt.contour(xgrd, zgrd, plotrho,
# levels=[24, 25, 26, 26.5, 27, 27.5, 28], colors='k')
#plt.ylim((-1000,0))
plt.xlim((0, 800))
ax.set_xlabel('Cross-track dist (km)')
ax.set_ylabel('Pressure (dbar)')
ax2 = ax.twiny()
ax2.set_xticks(trackdist);
ax2.set_xticklabels([int(x) for x in lwdata.ID-9000],
rotation='vertical', size='x-small');
ax2.set_title(prefix[ii])
plt.tight_layout()
In [ ]:
plotvel = plotu**2 + plotv**2
plt.figure()
plt.plot(xgrd, plotvel)
In [ ]: