Notebook

In [1]:

import numpy as np
import netCDF4 as nc
import matplotlib.mlab as mlab
from matplotlib import pyplot as plt
from matplotlib import animation
import matplotlib.cm as cm
import glob
import os
import matplotlib.gridspec as gridspec
from salishsea_tools import (nc_tools,viz_tools,stormtools)

%matplotlib inline

In [2]:

bathy = nc.Dataset('/ocean/imachuca/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')

In [3]:

path = '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/'
files = glob.glob(path+'*/SalishSea_1d_*grid_T.nc')
files.sort(key=os.path.basename)
#38 files, 3 are size 5 and 35 are 10. estimate 350+15 = 365 days
files

Out[3]:

['/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/1jan5jan/SalishSea_1d_20030101_20030105_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/6jan10jan/SalishSea_1d_20030106_20030110_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/11jan20jan/SalishSea_1d_20030111_20030120_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/21jan30jan/SalishSea_1d_20030121_20030130_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/31jan9feb/SalishSea_1d_20030131_20030209_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/10feb19feb/SalishSea_1d_20030210_20030219_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/20feb1mar/SalishSea_1d_20030220_20030301_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/2mar11mar/SalishSea_1d_20030302_20030311_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/12mar21mar/SalishSea_1d_20030312_20030321_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/22mar31mar/SalishSea_1d_20030322_20030331_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/1apr10apr/SalishSea_1d_20030401_20030410_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/11apr20apr/SalishSea_1d_20030411_20030420_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/21apr30apr/SalishSea_1d_20030421_20030430_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/1may10may/SalishSea_1d_20030501_20030510_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/11may20may/SalishSea_1d_20030511_20030520_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/21may30may/SalishSea_1d_20030521_20030530_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/31may9jun/SalishSea_1d_20030531_20030609_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/10jun19jun/SalishSea_1d_20030610_20030619_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/20jun29jun/SalishSea_1d_20030620_20030629_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/30jun9jul/SalishSea_1d_20030630_20030709_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/10jul19jul/SalishSea_1d_20030710_20030719_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/20jul29jul/SalishSea_1d_20030720_20030729_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/30jul8aug/SalishSea_1d_20030730_20030808_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/9aug18aug/SalishSea_1d_20030809_20030818_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/19aug28aug/SalishSea_1d_20030819_20030828_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/29aug7sep/SalishSea_1d_20030829_20030907_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/8sep17sep/SalishSea_1d_20030908_20030917_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/18sep27sep/SalishSea_1d_20030918_20030927_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/28sep7oct/SalishSea_1d_20030928_20031007_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/8oct17oct/SalishSea_1d_20031008_20031017_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/18oct27oct/SalishSea_1d_20031018_20031027_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/28oct6nov/SalishSea_1d_20031028_20031106_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/7nov16nov/SalishSea_1d_20031107_20031116_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/17nov26nov/SalishSea_1d_20031117_20031126_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/27nov6dec/SalishSea_1d_20031127_20031206_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/7dec16dec/SalishSea_1d_20031207_20031216_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/17dec26dec/SalishSea_1d_20031217_20031226_grid_T.nc',
 '/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/27dec31dec/SalishSea_1d_20031227_20031231_grid_T.nc']

In [4]:

tracers = nc.Dataset(files[0])
nc_tools.show_variables(tracers)

[u'nav_lon', u'nav_lat', u'deptht', u'time_counter', u'time_counter_bnds', u'sossheig', u'votemper', u'vosaline', u'rain_rate', u'snow_rate']

In [5]:

tracers = nc.MFDataset(files)

In [6]:

timesteps = tracers.variables['time_counter']
sal = tracers.variables['vosaline']
ssh = tracers.variables['sossheig']
tem = tracers.variables['votemper']
zlv = tracers.variables['deptht']
zlv = -zlv[:]

In [7]:

tem.shape

Out[7]:

(365, 40, 898, 398)

In [8]:

!head thalweg.txt
thalweg = np.loadtxt('thalweg.txt', dtype=int, unpack=True)
thalweg[1].shape

Out[8]:

(718,)

Tracers¶

In [9]:

def sal_fxn(ax, sal, t, part):
    sal0 = sal[:,0,:,:]
    cmin = sal0.min()
    cmax = sal0.max()
    
    cmap = plt.get_cmap('ocean_r')
    cmap.set_bad('burlywood')
    ax.set_axis_bgcolor('burlywood')
    viz_tools.set_aspect(ax)
    
    sal_t = np.ma.masked_values(sal0[t], 0)
    #mesh = ax.pcolormesh(sal_t, cmap=cmap)
    levels = np.arange(0, 33, 4)
    mesh = ax.contourf(sal_t, levels, cmap=cmap, extend='both')
    ax.set_xlabel('x Index'); ax.set_ylabel('y Index')
    ax.set_title('{label}'.format(label=sal.long_name.title()))
    if part == 'init':
        cbar = fig.colorbar(mesh, ax=ax)
        cbar.set_label('[{units}]'.format(units=sal.units))
    else:
        pass

In [10]:

def tem_fxn(ax, tem, t, part):
    tem0 = tem[:,0,:,:]
    cmin = tem0.min()
    cmax = tem0.max()
    
    cmap = plt.get_cmap('gist_rainbow_r')
    cmap.set_bad('burlywood')
    ax.set_axis_bgcolor('burlywood')
    viz_tools.set_aspect(ax)

    tem_t = np.ma.masked_values(tem0[t], 0)
    #mesh = ax.pcolormesh(tem_t, cmap=cmap)
    levels = np.arange(round(cmin), round(cmax), 2)
    mesh = ax.contourf(tem_t, levels, cmap=cmap, extend='both')
    ax.set_xlabel('x Index'); ax.set_ylabel('y Index')
    ax.set_title('{label}'.format(label=tem.long_name.title()))
    if part == 'init':
        cbar = fig.colorbar(mesh, ax=ax)
        cbar.set_label('[{units}]'.format(units=tem.units))
    else:
        pass

In [11]:

def ssh_fxn(ax, ssh, t, part):
    cmin = ssh[:].min()
    cmax = ssh[:].max()
    
    cmap = plt.get_cmap('terrain_r')
    cmap.set_bad('burlywood')
    ax.set_axis_bgcolor('burlywood')
    viz_tools.set_aspect(ax)
    
    ssh_t = np.ma.masked_values(ssh[t], 0)
    #mesh = ax.pcolormesh(ssh_t, cmap=cmap)
    levels = np.arange(round(cmin, 1), round(cmax, 1)+0.5, 0.25)
    mesh = ax.contourf(ssh_t, levels, cmap=cmap, extend='both')
    ax.set_xlabel('x Index'); ax.set_ylabel('y Index')
    ax.set_title('{label}'.format(label=ssh.long_name.title()))
    if part == 'init':
        cbar = fig.colorbar(mesh, ax=ax)
        cbar.set_label('[{units}]'.format(units=ssh.units))
        cbar.set_ticks(np.arange(-1.4, 1.9, 0.4))
    else:
        pass

Thalweg¶

In [12]:

def twg_fxn(ax, grid):
    
    cmap = plt.get_cmap('winter_r')
    cmap.set_bad('burlywood')
    viz_tools.set_aspect(ax)
    ax.set_axis_bgcolor('burlywood')
    
    bathy = grid.variables['Bathymetry']
    x_slice = np.arange(bathy.shape[1])
    y_slice = np.arange(0, 898)
    
    ax.pcolormesh(x_slice, y_slice, bathy[y_slice, x_slice], cmap=cmap)
    ax.plot(thalweg[1], thalweg[0], linestyle='-', marker='.', color='red', label='Thalweg Points')
    ax.set_xlabel('x Index'); ax.set_ylabel('y Index')
    ax.set_title('Thalweg Points')

In [13]:

def Twg_fxn(ax, sal, zlv, t, part):
    
    cmap = plt.get_cmap('hsv')
    cmap.set_bad('burlywood')
    ax.set_axis_bgcolor('burlywood')
    
    x, z = np.meshgrid(np.arange(thalweg.shape[1]), zlv)
    thw = np.ma.masked_values(sal[t, :, thalweg[0], thalweg[1]], 0)
    cs = [26,27,28,29,30,30.2,30.4,30.6,30.8,31,32,33,34]
    mesh=ax.contourf(x, z, thw, cs, cmap=cmap, extend='both')
    ax.set_xlim(0, thalweg[0][-1])
    ax.set_xlabel('Thalweg Points'); ax.set_ylabel('Depth')
    ax.set_title('{label} Profile along Thalweg'.format(label=sal.long_name.title()))
    
    if part == 'init':
        cbar = fig.colorbar(mesh, ax=ax)
        cbar.set_ticks([26,27,28,29,30,30.2,30.4,30.6,30.8,31,32,33,34])
        cbar.set_label('[{units}]'.format(units=sal.units))
    else:
        pass

In [17]:

fig = plt.figure(figsize=(18,12))
gs = gridspec.GridSpec(2, 4, height_ratios=[1.5,1])
ax0 = fig.add_subplot(gs[0, 0])
ax1 = fig.add_subplot(gs[0, 1])
ax2 = fig.add_subplot(gs[0, 2])
ax3 = fig.add_subplot(gs[0, 3])
ax4 = fig.add_subplot(gs[1, :])
twg_fxn(ax3, bathy)

def init():
    sal_fxn(ax0, sal, 0, 'init')
    tem_fxn(ax1, tem, 0, 'init')
    ssh_fxn(ax2, ssh, 0, 'init')
    Twg_fxn(ax4, sal, zlv, 0, 'init')
    timestamp = nc_tools.timestamp(tracers,0)
    plt.suptitle(timestamp.format('MMM DD'), fontsize=20, verticalalignment='bottom')
    
def animate(t):
    sal_fxn(ax0, sal, t, 'anim')
    tem_fxn(ax1, tem, t, 'anim')
    ssh_fxn(ax2, ssh, t, 'anim')
    Twg_fxn(ax4, sal, zlv, t, 'anim')
    timestamp = nc_tools.timestamp(tracers,t)
    plt.suptitle(timestamp.format('MMM DD'), fontsize=20, verticalalignment='bottom')

anim = animation.FuncAnimation(fig, animate, init_func=init, frames=365)
mywriter = animation.FFMpegWriter()
anim.save('Seminar_Profiles.mp4',writer=mywriter,fps=1, dpi=200, bitrate=1000000, codec="libx264")

jan 31
feb 28
mar 31
apr 30
may 31
jun 30
jul 31
aug 31
sep 30
oct 31
nov 30
dec 31

In [14]: