Plots comparing sea surface heights for an all forcing run and a tides only run over a period of 2 days and 1 hour intervals.¶

Imports¶

In [1]:

from __future__ import division, print_function

%matplotlib inline

import matplotlib.pyplot as plt

import numpy as np

import netCDF4 as nc

from salishsea_tools import (
    nc_tools,
    viz_tools,)

Result File¶

1. All Forcing¶

Hourly values of T-grid variables (i.e. sea surface height) for Dec 14 - Dec 15, 2006.

In [2]:

tracers1 = nc.Dataset('/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/dec2006/all_forcing/1hour/SalishSea_1h_20061214_20061215_grid_T.nc')

2. Tides Only¶

Hourly values of T-grid variables (i.e. sea surface height) for Dec 14 - Dec 15, 2006.

In [3]:

tracers2 = nc.Dataset('/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/dec2006/tidesonly/1hour/SalishSea_1h_20061214_20061215_grid_T.nc')

Dimensions of the dataset: 2 lateral dimensions, depth size = 40, time counter = 48 (hourly results for 2 days), time counter bounds size = 2.

In [4]:

nc_tools.show_variables(tracers1)
nc_tools.show_variables(tracers2)

[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']
[u'nav_lon', u'nav_lat', u'deptht', u'time_counter', u'time_counter_bnds', u'sossheig', u'votemper', u'vosaline']

Variables: nav_lon=longitude, nav_lat=latitude, deptht=depths, time_counter=time values, time counter_bnds=start and end times of each interval in time_counter, and sossheig=sea surface height (contains 2 dimensions). Variables have the same dimensions in both files.

In [5]:

timesteps = tracers1.variables['time_counter']
#timesteps will be used for plot

Defining Sea Surface Height¶

Sossheigh(time_counter=48 hrs, y=898, x=398). Interval is output interval=3600 s. Interval operation is time step for calculation = 10 s. FillValue and filling off means ssh data does not automatically mask out land areas (land = 0).

In [6]:

#Defining sea surface height as ssh1 and ssh2 for "all forcing" and "tides only" respectively
ssh1 = tracers1.variables['sossheig']
ssh2 = tracers2.variables['sossheig']

Function: Sea Surface Height Map¶

Defining a function that creates a sea surface height map for any time index. Here we only look at ssh1.

In [7]:

def ssh_fxn(timeindex):
    
    ssh1_t = np.ma.masked_values(ssh1[timeindex], 0)

    fig, ax = plt.subplots(1, 1, figsize=(10, 8))
    viz_tools.set_aspect(ax)
    cmap = plt.get_cmap('CMRmap')
    cmap.set_bad('burlywood')
    mesh = ax.pcolormesh(ssh1_t, cmap=cmap)
    cbar = fig.colorbar(mesh)

    ax.set_xlabel('x Index')
    ax.set_ylabel('y Index')
    ax.set_title('t = {:.1f}h'.format((timesteps[timeindex] / 3600)-72))
    cbar.set_label('{label} [{units}]'.format(label=ssh1.long_name.title(), units=ssh1.units))

Without a land mask, the land is shaded with the colour representing 0 for that time. The numpy masked array fills in the missing land data. the index for ssh belongs to the hour we're plotting. first hour = 0. last hour = 47.

In [8]:

timeindex=1
ssh_fxn(timeindex)

Row 1: All forcing. Row 2: tides only. Row 3: difference.

Sea Surface Height: Comparing runs¶

In [11]:

timesteps_comp = (0, 1, 2)
ssh_comp = (ssh1,ssh2)
fig, axs = plt.subplots(1, 3, figsize=(16, 8), sharey=True)


for ax, ssh in zip(axs, ssh_comp): 
    fig, axs = plt.subplots(1, 3, figsize=(16, 8), sharey=True)
    cmap = plt.get_cmap('jet')
    cmap.set_bad('burlywood')
    for ax, t in zip(axs, timesteps_comp):
        ssh_t = np.ma.masked_values(ssh[t], 0)
        viz_tools.set_aspect(ax)
        mesh = ax.pcolormesh(ssh_t, cmap=cmap)
        mesh.set_clim(vmin=-0.75,vmax=0.3)
        cbar = fig.colorbar(mesh, ax=ax)
        ax.set_xlabel('x Index')
    axs[0].set_ylabel('y Index')
        

        
fig, axs = plt.subplots(1, 3, figsize=(16, 8), sharey=True)
cmap = plt.get_cmap('jet')
cmap.set_bad('burlywood')
for ax, t in zip(axs, timesteps_comp):
    ssh3_t = np.ma.masked_values((ssh1[t]-ssh2[t]), 0)
    viz_tools.set_aspect(ax)
    mesh = ax.pcolormesh(ssh3_t, cmap=cmap)
    mesh.set_clim(vmin=-0.75,vmax=0.3)
    cbar = fig.colorbar(mesh, ax=ax)
    ax.set_xlabel('x Index')

axs[0].set_ylabel('y Index')

Out[11]:

<matplotlib.text.Text at 0x7fe1e541b810>

In [9]: