Density animation at y = 350 and y = 180 through various wind events.
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from salishsea_tools import (nc_tools, viz_tools, gsw_calls)
import numpy.ma as ma
from matplotlib import animation, rc
%matplotlib inline
mesh_mask = nc.Dataset('/home/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask_downbyone2.nc')
new_domain = mesh_mask.variables['tmask'][0,:,334:898,114:398]
def calc_rho(Sal, TempC, P):
""" Calculate rho: Based on SOG code
"""
# Calculate the square root of the salinities
sqrSal = np.sqrt(Sal)
# Calculate the density profile at the grid point depths
# Pure water density at atmospheric pressure
# (Bigg P.H., (1967) Br. J. Applied Physics 8 pp 521-537)
R1 = ((((6.536332e-9 * TempC - 1.120083e-6) * TempC + 1.001685e-4)
* TempC - 9.095290e-3) * TempC + 6.793952e-2) * TempC - 28.263737
R2 = (((5.3875e-9 * TempC - 8.2467e-7) * TempC + 7.6438e-5)
* TempC - 4.0899e-3) * TempC + 8.24493e-1
R3 = (-1.6546e-6 * TempC + 1.0227e-4) * TempC - 5.72466e-3
# International one-atmosphere equation of state of seawater
SIG = (4.8314e-4 * Sal + R3 * sqrSal + R2) * Sal + R1
# Specific volume at atmospheric pressure
V350P = 1.0 / 1028.1063
SVA = -SIG * V350P / (1028.1063 + SIG)
# Density anomoly at atmospheric pressure
rho = 28.106331 - SVA / (V350P * (V350P + SVA)) + 1000
return rho
Jan04 = nc.Dataset('/ocean/vdo/MEOPAR/completed-runs/SalishSeaLake/LakeJan0.4/SalishSea_1h_20170101_20170107_grid_T.nc')
fig, ax = plt.subplots(1,1, figsize=(12,10))
ax.pcolormesh(ma.masked_array(Jan04.variables['sossheig'][0,:,:], mask = 1-new_domain[0,:,:]))
viz_tools.set_aspect(ax)
ax.plot([0,283], [180, 180], 'r-')
ax.plot([0,283], [350,350], 'r-')
[<matplotlib.lines.Line2D at 0x7fb3b307bcc0>]
mask = 1 - new_domain[:,350,:]
mask2 = 1 - new_domain[:,180,:]
rc('animation', html='html5')
def animatedensity(file, yslice, xlims, ylims, mesh_mask):
pressure1 = gsw_calls.generic_gsw_caller('gsw_p_from_z.m',
[-np.expand_dims(file.variables['deptht'][:], 1) * np.ones(284),
(np.ones([284,40]) *np.expand_dims(file.variables['nav_lat'][yslice,:][:],1)).T])
fig,ax = plt.subplots(figsize=(8,8))
deptht = file.variables['deptht'][:]
yslicemask = 1 - mesh_mask[:,yslice,:]
def animate30(i):
ax.clear()
density1 = calc_rho(file.variables['vosaline'][i,:,yslice,:],
file.variables['votemper'][i,:,yslice,:],
pressure1)
masked_density1 = ma.masked_array(density1, mask = yslicemask)
den = ax.contourf(np.arange(0,284), deptht, masked_density1,
levels = np.linspace(1021.5, 1024.5, 11) #, colors='black'
)
ax.set_ylim(ylims)
ax.set_xlim(xlims)
ax.set_title('Y = ' + str(yslice) + ', hour = %03d'%(i))
return ax
interval = 0.25#in seconds
ani40 = animation.FuncAnimation(fig,animate30,frames=144,interval=interval*1e+3, repeat=False)
return ani40
animatedensity(Jan04, 180, (75,258), (130,0), new_domain)
animatedensity(Jan04, 350, (10,85), (130,0), new_domain)
JanW04 = nc.Dataset('/ocean/vdo/MEOPAR/completed-runs/SalishSeaLake/LakeJan-0.4/SalishSea_1h_20170101_20170107_grid_T.nc')
animatedensity(JanW04, 180, (75,258), (130,0), new_domain)
animatedensity(JanW04, 350, (10,85), (130,0), new_domain)
Jun04 = nc.Dataset('/ocean/vdo/MEOPAR/completed-runs/SalishSeaLake/LakeJun0.4/SalishSea_1h_20160601_20160607_grid_T.nc')
animatedensity(Jun04, 180, (75,258), (130,0), new_domain)
animatedensity(Jun04, 350, (10,85), (130,0), new_domain)
JunW04 = nc.Dataset('/ocean/vdo/MEOPAR/completed-runs/SalishSeaLake/LakeJun-0.4/SalishSea_1h_20160601_20160607_grid_T.nc')
animatedensity(JunW04, 180, (75,258), (130,0), new_domain)
animatedensity(JunW04, 350, (10,85), (130,0), new_domain)