import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from matplotlib import animation
from cmocean import cm
import os
import gsw
from datetime import datetime, timedelta
from dateutil.parser import parse
from salishsea_tools import viz_tools, geo_tools
from itertools import repeat
from warnings import simplefilter
from tqdm import tqdm_notebook as tqdm
from dynmodes import dynmodes
from IPython.display import HTML
%matplotlib inline
plt.rcParams['font.size'] = 12
plt.rcParams['animation.html'] = 'html5'
simplefilter('ignore')
Build HRDPS mask
def build_HRDPS_mask(lon_HRDPS, lat_HRDPS, lon_NEMO, lat_NEMO, mask_NEMO):
"""
"""
# Preallocate
m, n = lon_HRDPS.shape
mask_HRDPS = np.zeros(m*n, dtype=int)
# Evaluate each point on GEM grid
for index, lon, lat in zip(tqdm(range(m*n)), lon_HRDPS.flatten() - 360, lat_HRDPS.flatten()):
j, i = geo_tools.find_closest_model_point(lon, lat, lon_NEMO, lat_NEMO)
if j is np.nan or i is np.nan:
mask_HRDPS[index] = 0
else:
mask_HRDPS[index] = mask_NEMO[j, i]
# Reshape
mask_HRDPS = mask_HRDPS.reshape(m, n)
return mask_HRDPS
Scaling diagram
def scaling_diagram():
"""Plot scaling diagram
"""
x = np.arange(0, 1.01, 0.01)
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.xaxis.set_ticks([0, 0.1, 0.3, 1])
ax.yaxis.set_ticks([0, 0.7, 0.9, 1])
ax.xaxis.set_ticklabels([0, '$x_u$', '$R$', '$L$'])
ax.yaxis.set_ticklabels(['$H$', '$z_h$', '$z_u$', 0])
ax.plot([0.0, 1.0], [0.7, 0.7], 'k-')
#ax.plot([0.0, 0.6], [0.9, 0.7], 'm--')
#ax.plot([0.1, 0.6], [1.0, 0.7], 'm--')
ax.plot(x, 0.2 * np.exp(-x / 0.3) + 0.7, 'k--')
ax.plot(x, 0.3 * np.exp((0.1 - x) / 0.3) + 0.7, 'k--')
ax.text(0.63, 0.89, '$\\frac{\\tau}{\\rho_0f}$', size=18)
ax.text(0.05, 0.75, '$\\frac{z_u}{z_h} < 1$', size=18)
ax.text(0.27, 0.91, '$\\frac{z_u}{z_h} > 1$', size=18)
ax.arrow(0.6, 0.95, 0.1, 0, head_width=0.02, color='k')
return fig
Plot domain
def plot_domain(x, y, tmask, window, gs=None, figsize=(10, 10)):
"""Plot domain in plan view (either xy or lon/lat)
"""
if gs is not None:
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(gs[:, 0])
else:
fig, ax = plt.subplots(1, 1, figsize=figsize)
ax.contourf(x, y, tmask, levels=[-0.01, 0.01], colors='Burlywood')
ax.contour(x, y, tmask, levels=[-0.01, 0.01], colors='k')
ax.set_xlim(window[:2])
ax.set_ylim(window[2:])
viz_tools.set_aspect(ax)
return fig, ax
Make scaling figure
def make_scaling_figure():
"""Make scaling figure layout
"""
# Make figure layout
fig, axs = plt.subplots(2, 3, figsize=(17, 10))
xlabel = '$\\frac{\\tau T}{\\rho z_hfR}$'
ylabels = ['$z_u/z_h$', '$x_u/R$']
for i, row, ylabel in zip(range(2), axs, ylabels):
for ax, hour in zip(row, hours):
ax.set_xlim([0, 8])
ax.set_ylim([0, 5])
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if i == 0:
ax.plot([0, 1], [0, 1], 'k--')
ax.plot([1, 8], [1, 1], 'k--')
ax.set_title(f'{hour} hours')
else:
ax.plot([1, 6], [0, 5], 'k--')
ax.plot([1, 3.5], [0, 5], 'm--')
return fig, axs
Make hindcast prefix
def make_prefix(date, paths, res='h'):
"""Make prefix
"""
path = paths['NEMO']
if 'date_cutoff' in paths and date >= parse(paths['date_cutoff']):
path = paths['NEMO_cutoff']
fn = '_'.join([f'SalishSea_1{res}', *repeat(date.strftime('%Y%m%d'), 2)])
prefix = os.path.join(path, date.strftime('%d%b%y').lower(), fn)
return prefix
Rossby radius
def calc_deformation_radius(rho, xy, e1t=440, z_h=None, angle=np.pi/2, tol=10):
"""Calculate the baroclinic Rossby deformation radius
"""
# Scale factor
scalefac = np.sin(angle) * e1t
# 2-layer stratification
if z_h is not None:
g_prime = const['g'] * (const['rho_0'] - rho) / const['rho_0']
R = np.sqrt(g_prime * z_h * (const['H'] - z_h) / const['H']) / const['f']
H, H_old = const['H'], 0
while abs(H - H_old) > tol:
H_old = H
H = bathy.Bathymetry[xy[1], xy[0]:]
H = float(H[~np.isnan(H)][:int(R / scalefac)].mean())
R = np.sqrt(g_prime * z_h * (H - z_h) / H) / const['f']
# Continuous stratification
else:
N_int = np.sqrt(const['g'] / const['rho_0'] * np.diff(rho, axis=0) * np.diff(depth2d, axis=0)).sum(axis=0)
R_full = (N_int / (np.pi * const['f'])).compressed()
R, R_old = np.median(R_full), 0
while abs(R - R_old) > tol:
R_old = R
R = np.median(R_full[:max(1, int(R / scalefac))])
return R
Calculate density ρ
def calc_rho(data, depth, loc, tmask, time=False):
"""Calculate the density, rho
"""
if time: mskslc = slice(1, None)
else: mskslc = slice(None)
condition = np.broadcast_to(tmask[loc[mskslc]], data[tracers[0]][loc].shape) == 0
rho = gsw.rho(*[np.ma.masked_where(condition, data[k][loc]) for k in tracers], depth)
return rho
Definitions
# Definitions
HRDPS = xr.open_dataset('https://salishsea.eos.ubc.ca/erddap/griddap/ubcSSaSurfaceAtmosphereFieldsV1')
HRDPS_grid = xr.open_dataset('https://salishsea.eos.ubc.ca/erddap/griddap/ubcSSaAtmosphereGridV1')
NEMO_grid = xr.open_dataset('https://salishsea.eos.ubc.ca/erddap/griddap/ubcSSnBathymetryV17-02')
bathy = xr.open_dataset('/data/bmoorema/MEOPAR/grid/bathymetry_201702.nc')
mask = xr.open_dataset('/data/bmoorema/MEOPAR/grid/mesh_mask201702.nc')
const = {'g': 9.81, 'f': 1.11e-4, 'rho_0': 1024, 'H': 200, 'R': 5e3, 'z_h': 10}
depth2d = np.expand_dims(mask.gdept_1d[0, :], axis=1)
tracers = ['vosaline', 'votemper']
subdomain = [114, 398, 334, 898]
hours = [12, 24, 48]
sections = [
{'y': 450, 'x': (235, 315), 'yw': (137, 145), 'xw': (141, 149), 'a': np.arctan(10), 'c': 'r'},
{'y': 500, 'x': (205, 295), 'yw': (145, 153), 'xw': (134, 141), 'a': np.arctan( 1), 'c': 'darkorange'},
{'y': 540, 'x': (175, 255), 'yw': (151, 158), 'xw': (127, 133), 'a': np.arctan( 1), 'c': 'gold'},
{'y': 590, 'x': (130, 205), 'yw': (159, 166), 'xw': (114, 121), 'a': np.arctan( 2), 'c': 'darkslategray'},
{'y': 640, 'x': (130, 180), 'yw': (169, 174), 'xw': (109, 116), 'a': np.arctan(10), 'c': 'c'},
{'y': 680, 'x': (130, 200), 'yw': (177, 181), 'xw': (108, 115), 'a': np.arctan( 5), 'c': 'cyan'},
]
mask_HRDPS = build_HRDPS_mask(
HRDPS_grid.longitude.values - 360,
HRDPS_grid.longitude.values,
NEMO_grid.longitude.values,
NEMO_grid.latitude.values,
mask.tmask[0, 0, ...].values
)
HBox(children=(IntProgress(value=0, max=68096), HTML(value='')))
# Plot scaling diagram
fig = scaling_diagram()
fig.savefig('/home/bmoorema/Desktop/scaling_diagram.pdf', bbox_inches='tight')
By conservation of mass
When zh/zh<1
∫∞0zuexp{−xR}dx=2zuR2=zuR⏟upwelling area=1ρ0f∫tFt0τdt⏟cumulative Ekman transportwhere
R=1f(g′zh(H−zh)H)12Thus
zuzh=1ρ0zhfR∫tFt0τdtWhen zu/zh=1
∫∞xuzhexp{xu−xR}dx+zhxu=zh(R+xu)⏟upwelling area=1ρ0f∫tFt0τdt⏟cumulative Ekman transportor
zh(R+xu2)=1ρ0f∫tFt0τdtConsidering again R
xuR=1ρ0zhfR∫tFt0τdt−1or
xuR=2ρ0zhfR∫tFt0τdt−2# Plot section locations
gs = plt.GridSpec(6, 2, hspace=0.3)
fig, ax = plot_domain(mask.x, mask.y, mask.tmask[0, 0, ...], [110, 350, 350, 800], gs=gs, figsize=(12, 10))
# Plot section bathymetry
for row, sec in zip(range(6), sections):
line = []
for coord, slope in zip([sec['x'][0], sec['y']], [np.tan(sec['a']), -1]):
line.append([coord + w / slope for w in [-10, 10]])
ax_sec = fig.add_subplot(gs[5-row, 1])
ax.plot(sec['x'], [sec['y'], sec['y']], '-', color=sec['c'])
ax.plot(*line, '--', linewidth=2, color='lime')
ax_sec.contourf(mask.x, mask.gdept_1d[0, :], mask.tmask[0, :, sec['y'], :], levels=[-0.01, 0.01], colors='Burlywood')
ax_sec.contour(mask.x, mask.gdept_1d[0, :], mask.tmask[0, :, sec['y'], :], levels=[-0.01, 0.01], colors='k')
ax_sec.set_xlim(sec['x'])
ax_sec.set_ylim([420, 0])
Plot wind sections
xpoints = [
(141, 149),
(134, 141),
(127, 133),
(114, 121),
(109, 116),
(108, 115),
]
ypoints = [
(137, 145),
(145, 153),
(151, 158),
(159, 166),
(169, 174),
(177, 181),
]
xpoints = [np.round(np.linspace(*x, 10)).astype(int) for x in xpoints]
ypoints = [np.round(np.linspace(*y, 10)).astype(int) for y in ypoints]
# Plot wind secions
tmask = np.ma.masked_where(mask.nav_lon == 0, mask.tmask[0, 0, ...])
fig, ax = plot_domain(mask.nav_lon, mask.nav_lat, tmask, [-125.2, -123, 48.9, 50.1], figsize=(10, 6))
for sec, xpoint, ypoint in zip(sections, xpoints, ypoints):
ax.plot(mask.nav_lon[sec['y'], slice(*sec['x'])], mask.nav_lat[sec['y'], slice(*sec['x'])], '-', color=sec['c'])
ax.plot(
HRDPS_grid.longitude.values[ypoint, xpoint] - 360, HRDPS_grid.latitude.values[ypoint, xpoint],
'o', color=sec['c'], markeredgecolor='k'
)
#for n, x in enumerate(sec['xw']):
# lonlat = (HRDPS_grid.longitude[sec['yw'], x].values - 360, HRDPS_grid.latitude[sec['yw'], x].values)
# ax.plot(*lonlat, 'ko', markerfacecolor=sec['c'])
# Define path and filename string
path = '/data/bmoorema/results/Lake/S4d'
fn = 'SalishSeaIdeal_1h_20170701_20170706_grid_T.nc'
tmask = mask.tmask[0, :, slice(*subdomain[2:]), slice(*subdomain[:2])]
# Define scaling parameters
param = {
'u_wind': [5, 10, 15],
'tau': [0.034, 0.123, 0.280],
'z_h': [10, 15, 20],
'rho_surf': [1016, 1018, 1020],
}
# Initialize storage lists
R = [0, 0, 0, 0, 0, 0]
Idealized scaling calculations
def calc_scaling_ideal():
"""Calculate scaling from idealized runs
"""
# Make figure layout
fig, axs = make_scaling_figure()
# Loop though z_h and rho_surf
for z_h in tqdm(param['z_h']):
for rho_s in param['rho_surf']:
# Calculate deformation radius
for i, sec in zip(range(6), sections):
R[i] = calc_deformation_radius(rho_s, (sec['x'][0], sec['y']), z_h=z_h, angle=sec['a'])
# Loop through wind speed
for U, tau in zip(param['u_wind'], param['tau']):
# Open results record
runID = f'SalishSeaPond_S4d{U:02d}ms_halocline{z_h:2d}m_rhosurf{rho_s:4d}'
with xr.open_dataset(os.path.join(path, runID, fn)) as data:
# Loop through sections
for sec, r in zip(sections, R):
scale = np.cos(np.arctan(2) - sec['a']) * tau * 3600 / (const['rho_0'] * z_h * const['f'] * r)
xy = (sec['y'] - subdomain[2], slice(*[x - subdomain[0] for x in sec['x']]))
rho_t0 = calc_rho(data, depth2d, (0,slice(None))+xy, tmask, time=True)
rho_halo = np.median(rho_t0[abs(depth2d[:, 0] - z_h).argmin(), :])
# Loop through hours
for col, hour in zip(axs.T, hours):
scale_T = (hour - 3) * scale
rho = calc_rho(data, 0, (hour,0)+xy, tmask, time=True)
z_u = depth2d[int(np.median(abs(rho_t0 - rho.max()).argmin(axis=0))), 0] / z_h
x_u = np.sin(sec['a']) * (rho >= rho_halo).sum() * 440 / r
col[0].plot(scale_T, z_u, 'o', color=sec['c'], markeredgecolor='k')
col[1].plot(scale_T, x_u, 'o', color=sec['c'], markeredgecolor='k')
return fig
fig = calc_scaling_ideal()
HBox(children=(IntProgress(value=0, max=3), HTML(value='')))
# Define paths and variables
paths = {
'NEMO': '/results/SalishSea/hindcast.201812',
'NEMO_cutoff': '/results2/SalishSea/hindcast.201812_annex',
'date_cutoff': '2016 Nov 21',
}
tmask = mask.tmask
depthw = (depth2d[1:, 0] + depth2d[:-1, 0]) / 2
# Initialize storage lists
rho_t0 = [0, 0, 0, 0, 0, 0]
rho_halo = [0, 0, 0, 0, 0, 0]
R = [0, 0, 0, 0, 0, 0]
z_h = [0, 0, 0, 0, 0, 0]
Define upwelling event start dates
# Upwelling event date windows
events = [parse(date) for date in [
'2015 Mar 11 04:00', #1
'2015 Mar 18 18:00', #2
'2015 Mar 24 22:00', #3
'2015 Mar 28 21:00', #4
'2015 Apr 25 23:00', #5
'2016 Apr 11 13:00', #6
'2017 Apr 03 23:00', #7
'2017 May 30 15:00', #8
'2017 Sep 08 23:00', #9
'2017 Sep 17 02:00', #10
'2017 Oct 21 05:00', #11
'2018 Mar 22 07:00', #12
'2018 Apr 03 10:00', #13
'2018 Apr 12 23:00', #14
'2018 Sep 07 13:00', #15
'2018 Sep 14 18:00', #16
'2018 Oct 23 06:00', #17
]]
Hindcast scaling calculations
def calc_scaling_hindcast():
"""Calculate scaling from hindcast record
"""
# Make figure layout
fig, axs = make_scaling_figure()
# Loop though events and hours
for event in tqdm(events):
for col, hour in enumerate([0] + hours):
date = event + timedelta(hours=hour)
# Open hindcast record
with xr.open_dataset(make_prefix(date, paths) + '_grid_T.nc') as data:
# Loop through sections
for n, sec in enumerate(sections):
yx = (sec['y'], slice(*sec['x']))
e1t = float(mask.e1t[0, sec['y'], sec['x'][0]])
# Calculations at t=0
if hour == 0:
loc = (slice(None),) + yx
rho_t0[n] = calc_rho(data.sel(time_counter=date, method='nearest'), depth2d, loc, tmask[0, ...])
N2 = const['g'] / const['rho_0'] * np.diff(rho_t0[n], axis=0) / np.diff(depth2d, axis=0)
modes = dynmodes(np.ma.median(N2[:, :20], axis=1).compressed()[:24], depthw[:24], 1)
z_h[n] = depthw[int(np.where(np.diff(np.signbit(modes[1][0, :])))[0])]
rho_halo[n] = np.median(rho_t0[n][abs(depth2d[:, 0] - z_h[n]).argmin(axis=0), :])
rho_surf = rho_t0[n][:int(z_h[n]), :].mean()
R[n] = calc_deformation_radius(rho_surf, (sec['x'][0], sec['y']), e1t=e1t, z_h=z_h[n], angle=sec['a'])
# Calculations at hours
else:
# Calculate upwelling metrics
rho = calc_rho(data.sel(time_counter=date, method='nearest'), 0, (0,) + yx, tmask[0, ...])
rho_max = rho[:max(1, int(R[n] / e1t))].max()
z_u = depth2d[int(np.median(abs(rho_t0[n] - rho_max).argmin(axis=0))), 0]
x_u = np.sin(sec['a']) * (rho >= rho_halo[n]).sum() * e1t
# Calculate scaling
#loc = (slice(None), sec['yw'], slice(*sec['xw']))
xlocs = np.round(np.linspace(*sec['xw'], 10)).astype(int)
ylocs = np.round(np.linspace(*sec['yw'], 10)).astype(int)
u, v = [HRDPS.sel(time=slice(event, date))[k].values[:, ylocs, xlocs] for k in ['u_wind', 'v_wind']]
for angle, ec, o1 in zip([True, False], ['k', None], [1, 0]):
if angle: coeff = np.cos(np.pi * (22 / 180 + 1) - np.arctan2(v, u) - sec['a'])
else: coeff = 1
tau = sum(1.225e-3 * (coeff * (u**2 + v**2)).mean(axis=1)) * 3600
for z, r, c, o2 in zip([z_h[n], const['z_h']], [R[n], const['R']], [sec['c'], 'lightgray'], [2, 0]):
scaling = tau / (const['rho_0'] * z * const['f'] * r)
axs[0, col-1].plot(scaling, z_u / z, 'o', color=c, markeredgecolor=ec, zorder=o1+o2)
axs[1, col-1].plot(scaling, x_u / r, 'o', color=c, markeredgecolor=ec, zorder=o1+o2)
return fig
fig = calc_scaling_hindcast()
#fig.savefig('/home/bmoorema/Desktop/scaling_breakdown.pdf', bbox_inches='tight')
HBox(children=(IntProgress(value=0, max=17), HTML(value='')))
# Make figure layout
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.set_xlabel('Maximum $N^2$ profile depth [m]')
ax.set_ylabel('Halocline depth $z_h$ [m]')
# Loop though events and hours
for event in tqdm(events):
# Open hindcast record
with xr.open_dataset(make_prefix(event, paths) + '_grid_T.nc').sel(time_counter=event, method='nearest') as data:
# Loop through sections
for sec in [sections[2]]:
rho = calc_rho(data, depth2d, (slice(None), sec['y'], slice(*sec['x'])), tmask[0, ...])
N2 = const['g'] / const['rho_0'] * np.diff(rho, axis=0) / np.diff(depth2d, axis=0)
for index in range(10, 30):
modes = dynmodes(np.ma.median(N2[:, :20], axis=1).compressed()[:index], depthw[:index], 1)
for mode, color in zip(modes[1], ['k', 'r']):
z_h = depthw[int(np.where(np.diff(np.signbit(mode)))[0])]
ax.plot(depthw[index], z_h, 'o', color=color)
HBox(children=(IntProgress(value=0, max=17), HTML(value='')))
# Make figure layout
fig, axs = plt.subplots(4, 2, figsize=(17, 10))
# Loop though events and hours
for row, event in zip(tqdm(axs), events[:4]):
date = (event - timedelta(hours=6), event + timedelta(hours=24))
time = HRDPS.sel(time=slice(*date)).time
for ax, ylim in zip(row, [(0, 15), (-60, 300)]):
ax.set_xlim(date)
ax.set_ylim(ylim)
ax.xaxis.set_major_locator(mdates.HourLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%H'))
ax.plot([event, event], ylim, 'k--')
for sec in sections:
loc = (slice(None), sec['yw'], slice(*sec['xw']))
u, v = [HRDPS.sel(time=slice(*date))[k][loc].values for k in ['u_wind', 'v_wind']]
mag = np.sqrt(u**2 + v**2).mean(axis=1)
angle = 180 * (1 - (np.arctan2(v, u).mean(axis=1) + sec['a']) / np.pi) + 22
row[0].plot(time, mag, '-', color=sec['c'])
row[1].plot(time, angle, '-', color=sec['c'])
row[1].plot(date, [0, 0], 'k--')
HBox(children=(IntProgress(value=0, max=4), HTML(value='')))
u, z_h = 10, 10
path = '/data/bmoorema/results/Lake/N4d'
fn = 'SalishSeaIdeal_1h_20170701_20170706_grid_T.nc'
fig, axs = plt.subplots(1, 2, figsize=(17, 10))
for ax, rho_0 in zip(axs, [1016, 1020]):
runID = f'SalishSeaPond_N4d{u:02d}ms_halocline{z_h}m_rhosurf{rho_0}'
ds = xr.open_dataset(os.path.join(path, runID, fn))
ax.contourf(bathy.x[slice(*subdomain[:2])], bathy.y[slice(*subdomain[2:])], ds.tracer20m[48, 0, ...], cmap='Purples')
ax.contourf(bathy.x, bathy.y, mask.tmask[0, 0, ...], levels=[-0.01, 0.01], colors='Burlywood')
ax.contour(bathy.x, bathy.y, mask.tmask[0, 0, ...], levels=[-0.01, 0.01], colors='k')
ax.set_xlim([110, 350])
ax.set_ylim([400, 750])
viz_tools.set_aspect(ax)