In [3]:
%matplotlib inline
%load_ext secnum
%secnum
import numpy as np
import h5py
import netCDF4
import numexpr as ne
from scipy.io import loadmat
from datetime import datetime
import time
from IPython import display
import pylab
import pyroms
# 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
#mpl.rcParams['lines.linewidth'] = 1.5
# need wider colorbar
# change default colormap
# change default fonts
# zoomable inline plots
#mport mpld3
#mpld3.disable_notebook()
# test for estimating baroclinicity
def baroclinicity(u, z):
u = np.ma.masked_array(u, np.isnan(u))
umean = u.mean()
#urms = np.sqrt((u**2).mean())
#uarms = np.sqrt( ((u-umean)**2).mean() )
# integrate kinetic energy
ketot = np.abs( 0.5 * np.trapz(u**2, x=z.transpose()))
if np.ma.is_masked(u):
length = len(u.compressed())
else:
length = len(u)
# fill in with ones to make depth averaged velocity
keda = np.abs(0.5 * np.trapz((umean**2) * u/u, x=z.transpose()))
bc = (ketot-keda)/ketot
#bc = uarms/umean
if np.ma.is_masked(u):
#print( 'Total KE = {0:.2f} | DA KE = {1:0.2f}'.format(ketot[0], keda[0]) )
return bc[0]
else:
#print( 'Total KE = {0:.2f} | DA KE = {1:0.2f}'.format(ketot, keda) )
return bc
In [6]:
lbdir = 'data/leebrink2010/'
# For context, let's look at some SST, SSH plots
SSTname = lbdir + 'SST.nc4'
SSTmname = lbdir + 'SST-mean.nc4'
SSHname = lbdir + 'SSH-daily.nc'
bathyname = '../datasets/ETOPO2v2g_f4.nc'
sstfile = h5py.File(SSTname,'r')
sstmfile = h5py.File(SSTmname,'r');
bathyfile = netCDF4.Dataset(bathyname,'r')
class Bunch(object):
def __init__(self, **kwds):
self.__dict__.update(kwds)
sst = Bunch(lon=sstfile['lon'][:], lat=sstfile['lat'][:],
time=sstfile['time'][:], anom=sstfile['anom'],
mean=sstmfile['sst']);
sst.ntime = sst.anom.shape[0]
# ugh, for some reason it doesn't do the scale factor
sst.anom = sst.anom.attrs['scale_factor'] * sst.anom;
sst.anom = sst.anom + sstfile['anom'].attrs['add_offset']
sst.mean = sst.mean.attrs['scale_factor'] * sst.mean;
sst.mean = sst.mean + sstmfile['sst'].attrs['add_offset']
anomrange = sstfile['anom'].attrs['actual_range']
meanrange = sstmfile['sst'].attrs['actual_range']
# set time right
sstepoch = datetime(1800, 1, 1);
sst.time = sst.time + sstepoch.toordinal();
sst.lon = sst.lon - 360;
# get SSH data
sshfile = netCDF4.Dataset(SSHname, 'r')
sshepoch = datetime(1950,1,1)
ssh = Bunch(lon=sshfile.variables['NbLongitudes'][:], lat=sshfile.variables['NbLatitudes'][:],
time=sshfile.variables['time'][:], H=sshfile.variables['Grid_0001'][:]);
ssh.lon = ssh.lon-360;
ssh.time = sst.time + sshepoch.toordinal()
# get ETOPO data
bathy = Bunch(lon=bathyfile.variables['x'][:], lat=bathyfile.variables['y'][:],
z=bathyfile.variables['z'][:])
#bathy.lon = bathy.lon + 180 - 360;
# define plot ranges
def get_range(lon, lat):
import numpy as np
latmin = 36; latmax = 46;
lonmin = -72; lonmax = -60;
lats = np.nonzero( np.logical_and( lat < latmax, lat > latmin ))[0]
lons = np.nonzero( np.logical_and( lon < lonmax, lon > lonmin ))[0]
return (lons, lats)
# get plot ranges for all variables
sshlons, sshlats = get_range(ssh.lon, ssh.lat)
sstlons, sstlats = get_range(sst.lon, sst.lat)
bathylons, bathylats = get_range(bathy.lon, bathy.lat)
In [ ]:
%matplotlib inline
fig = plt.figure()
# reads SST, SSH & bathymetric data
for ii in range(10): #sst.ntime):
plt.clf()
# SST plot
plt.subplot(1,2,1)
plt.pcolormesh(sst.lon[sstlons], sst.lat[sstlats],
sst.anom[ii,sstlats,sstlons[0]:sstlons[-1]]
+ sst.mean[ii,sstlats,sstlons[0]:sstlons[-1]],
cmap='RdBu_r')
plt.colorbar()
plt.xlim(-72,-60)
plt.clim(anomrange + meanrange)
plt.hold(True)
datestr = 'Full SST ' + \
datetime.strftime( datetime.fromordinal(np.int(sst.time[ii])), '%d-%m-%Y')
plt.title(datestr)
plt.hold(True)
plt.contour(bathy.lon[bathylons], bathy.lat[bathylats],
bathy.z[bathylats,bathylons[0]:bathylons[-1]+1],
[ -200, -1000, -2000],
colors='k', linestyles='solid', alpha=0.3)
# SSH plot
plt.subplot(1,2,2)
plt.pcolor(ssh.lon[sshlons], ssh.lat[sshlats],
ssh.H[ii,sshlons,sshlats[0]:sshlats[-1]].transpose(), cmap='RdBu_r')
plt.clim(np.array([-1,1]) * abs(ssh.H[:]).max())
plt.title('SSH (cm)');
plt.colorbar()
plt.hold(True)
cs = plt.contour(bathy.lon[bathylons], bathy.lat[bathylats],
bathy.z[bathylats,bathylons[0]:bathylons[-1]+1],
[ -200, -1000, -1800],
colors='k', linestyles='solid', alpha=0.3)
#plt.clabel(cs)
#fig.show()
#print(ii)
#time.sleep(0.5);
# update inline images
display.clear_output()
display.display(plt.gcf())
plt.close()
In [ ]:
lbdir = 'data/leebrink2010/'
g3name = lbdir + 'globec3.mat'
g3uvname = lbdir + 'g3_t2_intensive_uv.mat'
# load GLOBEC3 data
g3 = loadmat(g3name)
g3uv = loadmat(g3uvname)
# first, process time
g3epoch = datetime(1997,1,1).toordinal()
g3time = g3['txy'][:,0]
g3uvtime = g3uv['dd'][:]
# figure out indices for each day
# store this in dictionary - g3days
# first for T,S
g3days = {}
for ii in range(178,185):
i1 = np.nonzero(g3time == ii)
i2 = np.nonzero(g3time == ii+1)
key = datetime.strftime( datetime.fromordinal(g3epoch + ii), '%d%b')
g3days[key] = np.arange(i1[0],i2[0]-1)
# now for u,v
g3uvdays = {}
for ii in range(179,185):
i1 = np.nonzero(g3uvtime > ii)[0]
i2 = np.nonzero(g3uvtime > ii+1)[0]
key = datetime.strftime( datetime.fromordinal(g3epoch + ii), '%d%b')
g3uvdays[key] = np.arange(i1[0],i2[0]-1)
# manually add in 27 June since first time record is at the middle of the 27th
g3days['27Jun'] = np.arange(0, g3days['28Jun'][0]-1)
g3uvdays['27Jun'] = 0
g3uvdays['28Jun'] = np.arange(0, g3uvdays['29Jun'][0]-1)
In [ ]:
# let's try a zslice as an experiment
day = '29Jun'
tsindices = g3days[day]
uvindices = g3uvdays[day]
plt.clf()
plt.scatter(g3['txy'][tsindices,1], g3['txy'][tsindices,2], c=g3['theta'][tsindices,8])
plt.colorbar()
plt.hold(True)
plt.quiver(g3uv['longitude'][uvindices], g3uv['latitude'][uvindices],
g3uv['u'][uvindices,1], g3uv['v'][uvindices,1])
For now, what I know is that the vortex on the shelf is barotropic - at least in vorticity. To estimate vorticity, I think I have 2 options:
How do I estimate baroclinicity?
Function baroclinicity defined earlier calculates the ratio of (Total K.E - Depth-averaged K.E)/Total K.E - i.e., it estimates the fraction of total KE contributed by the velocity anomalies (defined w.r.t the depth average).
In [ ]:
# Let's see how barotropic in (u,v) the filament is
# use Jul 04 data - I've manually picked out the points
ind1 = g3uvdays['04Jul'][270:290]
ind2 = g3uvdays['04Jul'][315:330]
indices = np.concatenate( (ind1,ind2), 1)
plt.figure()
# let's try a zslice as an experiment
day = '04Jul'
tsindices = g3days[day]
uvindices = indices
plt.clf()
plt.scatter(g3['txy'][tsindices,1], g3['txy'][tsindices,2], c=g3['theta'][tsindices,8])
plt.colorbar()
plt.hold(True)
plt.quiver(g3uv['longitude'][uvindices], g3uv['latitude'][uvindices],
g3uv['u'][uvindices,1], g3uv['v'][uvindices,1])
# depth profiles of velocities
plt.figure()
plt.subplot(2,2,1)
plt.title('u (m/s) - outside')
plt.hold(True)
plt.subplot(2,2,2)
plt.hold(True)
plt.subplot(2,2,3)
plt.hold(True)
plt.subplot(2,2,4)
plt.hold(True)
BCu1 = BCv1 = np.zeros(len(ind1))
BCu2 = BCv2 = np.zeros(len(ind2))
depth = -1*g3uv['depth']
for ii in ind1:
plt.subplot(2,2,1)
u = g3uv['u'][ii,:]
plt.plot(u, depth);
BCu1[ii-ind1[0]] = baroclinicity(u, depth)
plt.subplot(2,2,2)
v = g3uv['v'][ii,:]
plt.plot(v, depth)
BCv1[ii-ind1[0]] = baroclinicity(v, depth)
plt.subplot(2,2,1)
plt.title('u (m/s) out - BC = {0:.2f} +- {1:.2f}'.format(BCu1.mean(), BCu1.std()))
plt.subplot(2,2,2)
plt.title('v (m/s) out - BC = {0:.2f} +- {1:.2f}'.format(BCv1.mean(), BCv1.std()))
for ii in ind2:
plt.subplot(2,2,3)
u = g3uv['u'][ii,:]
plt.plot(u, depth);
BCu2[ii-ind2[0]] = baroclinicity(u, depth)
plt.subplot(2,2,4)
v = g3uv['v'][ii,:]
plt.plot(v, depth)
BCv2[ii-ind2[0]] = baroclinicity(v, depth)
plt.subplot(2,2,3)
plt.title('u (m/s) cen - BC = {0:.2f} +- {1:.2f}'.format(BCu2.mean(), BCu2.std()))
plt.subplot(2,2,4)
plt.title('v (m/s) cen - BC = {0:.2f} +- {1:.2f}'.format(BCv2.mean(), BCv2.std()))
for ii in range(1,5):
plt.subplot(2,2,ii)
plt.xlim(-1,0.2)
plt.tight_layout()
In [ ]:
# test on actual data
u = g3uv['u'][ind1[0]]
depth = -1*g3uv['depth'];
umean = np.ma.masked_array(u, np.isnan(u)).mean()
uarms = np.sqrt( ( (u-umean)**2).mean() )
plt.subplot(1,2,1)
plt.plot(u-umean, depth)
#plt.axvline(umean)
plt.axvline(0, color='k')
plt.title('BC = {0:.3f}'.format(baroclinicity(u, depth)))
# test for idealized BT + BC1 vel profile
z = -1*np.arange(1,90)
xzU = 1 + 1 * np.cos(np.pi * z/z.min()) #+ 0.1 * np.cos(2*np.pi * z/z.min())
Urms = np.sqrt( (U**2).mean())
Umean = Umean
#Umean = U.mean()
Uarms = np.sqrt( ((U-Umean)**2).mean() )
plt.subplot(1,2,2)
plt.plot(U-Umean,z)
plt.axvline(0, color='k')
plt.title('BC = {0:.3f}'.format(baroclinicity(U, z)))
Doing this by taking net velocity vector and rotating to get along-track and cross-track velocities.
There are some problems currently:
In [ ]:
pylab.rcParams['figure.figsize'] = (7.0, 4.0)
u = g3uv['u'][:]
v = g3uv['v'][:]
lon = g3uv['longitude'][:]
lat = g3uv['latitude'][:]
# calculate cruise track angle and distances using seawater toolbox
import seawater as sw
dist,angle = sw.dist(lat, lon, units='km')
# returns in km, convert that to m
dist = dist*1000
# returns angle in degrees, convert to radians
angle = np.pi/180 * angle;
daystr = '04Jul'
tsindices = g3days[daystr]
uvindices = g3uvdays[daystr]
# along-track velocity
us = u[1:] * np.cos(angle) - v[1:] * np.sin(angle)
# cross-track velocity
un = u[1:] * np.sin(angle) + v[1:] * np.cos(angle)
# convert to masked arrays
us = np.ma.masked_array(us, np.isnan(us))
un = np.ma.masked_array(un, np.isnan(un))
vor = np.diff(un, axis=0)/np.diff(dist, axis=0)
# map
plt.subplot(1,2,1)
plt.hold(True)
plt.scatter(g3['txy'][tsindices,1], g3['txy'][tsindices,2], c=g3['theta'][tsindices,8], alpha=1)
plt.colorbar()
plt.subplot(1,2,2)
plt.scatter(lon[uvindices], lat[uvindices], c=vor[uvindices,10]/1e-4, cmap='RdBu_r')
plt.clim(-50,50)
plt.colorbar()
# depth profiles
plt.figure()
plt.subplot(1,2,1)
for ii in ind1:
plt.plot(vor[ii,:]/1e-4, depth);
plt.xlim(-200,200)
plt.title('Vorticity - outside')
plt.subplot(1,2,2)
for ii in ind2:
plt.plot(vor[ii,:]/1e-4, depth)
plt.xlim(-200,200)
plt.title('Vorticity - center')