Notebook

This notebook will plot the spatial extent of the surge on Feb. 4 2006

In [1]:

from __future__ import division
import netCDF4 as NC
import numpy as np


%matplotlib inline
import matplotlib.pyplot as plt
from salishsea_tools import viz_tools
from salishsea_tools import stormtools
from salishsea_tools import nc_tools
import arrow
import datetime
import pytz
import matplotlib.dates as mdates
from datetime import timedelta

In [2]:

fB = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
lats = fB.variables['nav_lat']
lons = fB.variables['nav_lon']
D = fB.variables['Bathymetry']

In [3]:

#Load the data
path='/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/feb2006/'

runs = {'all_forcing','tidesonly', 'surgeonly','riversonly'}

fUs={}; fVs={}; fTs={};

for key in runs:
    fUs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_U.nc','r');
    fVs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_V.nc','r');
    fTs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_T.nc','r');

Wite out the times closest to Feb 4 2006 surge days

In [4]:

plot_ts= nc_tools.timestamp(fTs['all_forcing'],[18,19,20,21,22,23, 24,25,26])
plot_ts[-1]

Out[4]:

<Arrow [2006-02-05T10:00:00+00:00]>

In [5]:

fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20,21]
vmin=0; vmax=0.75
cmap='Reds'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Surge (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    
fig, axs = plt.subplots(1, 4, figsize=(20,6))
depthlevel=0
timesteps = [22,23,24,25]
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    cbar.set_label('Surge (m)')

Struggled with using axs and zip when there was more than one row of subplots.

These are 4 hour averages.

Surges highest in the inlets perhaps due to atmosperic pressure? Surge is also quite high in SoG, espeically the northern SoG. Islands and Victoria have smaller surges.

Compare with surgeonly case¶

In [6]:

fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20,21]
vmin=0; vmax=0.75
cmap='Reds'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['surgeonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    
fig, axs = plt.subplots(1, 4, figsize=(20,6))
depthlevel=0
timesteps = [22,23,24,25]
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['surgeonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])

Difference between anomaly and surgeonly run¶

Take the difference between modelled anomaly (all_forcing-tidesonly) and a run with surgeonly.

In [7]:

fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20,21]
vmin=-0.1; vmax=0.1
cmap='bwr'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly']-Es['surgeonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Surge (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    
fig, axs = plt.subplots(1, 4, figsize=(20,6))
depthlevel=0
timesteps = [22,23,24,25]
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly']-Es['surgeonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    cbar.set_label('Surge (m)')

Behaviour changes after the islands. Also in the North, after Seymour Narrows.

Generally, modlled anomaly higher than surgeonly run in JdF, but opposite in SoG. Differences are typically less than 10cm, except maybe in Johnstone Strait.

Are differences due to tide-surge interactions?

Timing of surge at Point Atkinson¶

In [8]:

start='01-Feb-2006'
stations = {'PointAtkinson': 7795, 'Victoria': 7120, 'PatriciaBay': 7277, 'CampbellRiver': 8074, 
            'CrescentBeach': 1, 'WhiteRock': 1, 'BoundaryBay': 1}
datums = {'PointAtkinson': 3.09, 'Victoria': 1.881, 'PatriciaBay': 2.256, 'CampbellRiver': 2.916, 
          'CrescentBeach': 2.44, 'WhiteRock': 2.85, 'BoundaryBay': 1 }
thalwegs = ['Thalweg1','Thalweg2', 'Thalweg3', 'Thalweg4', 'Thalweg5', 'Thalweg6']
#datums from t_xtide
run_stations={}
us={}; vs={}; lats={}; lons={}; tmps={}; sals={}; sshs={}; ts={};

for run in runs:
    for key in ['PointAtkinson']:
        string = path + run + '/1h_' + key + '.nc'
        run_stations[key] = NC.Dataset(string,'r');
        tim = run_stations[key].variables['time_counter']
        t_count=np.arange(0, tim.shape[0])
        t=nc_tools.timestamp(run_stations[key],t_count)
        tlist=[]
        for a in t:
           tlist.append(a.datetime) 
        ts[run]=tlist
    [us[run], vs[run], lats[run], lons[run], tmps[run], sals[run], sshs[run]] = stormtools.combine_data(run_stations)
    run_stations={}

In [9]:

s = '04-Feb-2006'; e='06-Feb-2006';
unaware=datetime.datetime.strptime(s,"%d-%b-%Y")
sax = unaware.replace(tzinfo=pytz.timezone('utc'))
unaware=datetime.datetime.strptime(e,"%d-%b-%Y")
eax = unaware.replace(tzinfo=pytz.timezone('utc'))

fig, ax = plt.subplots(1, 1, figsize=(15, 3))

hfmt = mdates.DateFormatter('%H')
ax.plot(ts['surgeonly'], sshs['surgeonly']['PointAtkinson'][:,0,0], 'k-',lw=1.5)
ax.set_xlim([sax,eax])
ax.set_ylim([-1.5,1.5])
ax.grid(True,lw=1)
ax.set_ylabel('Residual (m)')
ax.xaxis.set_major_locator(mdates.HourLocator())
ax.xaxis.set_major_formatter(hfmt)

Peak is around 20:00 on Feb 4.

Differences over time¶

How does the surge change over time? Plots of the difference in surge at two consecutive times (surgeonly run)

In [10]:

fig, axs = plt.subplots(1, 3, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20]
vmin=-.5; vmax=0.5
cmap='bwr'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};
    Us2={}; Vs2={}; Es2 ={}; Ss2={}; Ts2={};
    t2=t+1

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    for key in runs:
        [Us2[key], Vs2[key], Es2[key], Ss2[key], Ts2[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t2,depthlevel)
        
    mesh=ax.pcolormesh(Es2['surgeonly']-Es['surgeonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Difference (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title('Times: ' +str(plot_ts[t2-18].hour) +' -' +str(plot_ts[t-18].hour) )

timesteps = [21,22,23]
fig, axs = plt.subplots(1, 3, figsize=(20,6),sharey=True)
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};
    Us2={}; Vs2={}; Es2 ={}; Ss2={}; Ts2={};
    t2=t+1

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    for key in runs:
        [Us2[key], Vs2[key], Es2[key], Ss2[key], Ts2[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t2,depthlevel)
        
    mesh=ax.pcolormesh(Es2['surgeonly']-Es['surgeonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Difference (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title('Times: ' +str(plot_ts[t2-18].hour) +' -' +str(plot_ts[t-18].hour) )
    

Rivers only¶

Rivers only run : no CGRF atmopsheric, no tidal forcing but rivers still on (ln_rnf=true).

OBCs

2D velocities: flather with initial conditions
3D velocities: zero gradient
T+S: flow relexation

Check on ssh field.

In [11]:

fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20,21]
vmin=-0.1; vmax=0.1
cmap='bwr'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['riversonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('SSH (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    
fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [22,23,24,25]
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['riversonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    cbar.set_label('SSH (m)')

Doens't look to be changing over time. Probably because no tides.

Check that salinity field is changing over time.

In [12]:

fig, axs = plt.subplots(1, 3, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20]
vmin=-1; vmax=1
cmap='bwr'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};
    Us2={}; Vs2={}; Es2 ={}; Ss2={}; Ts2={};
    t2=t+1

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    for key in runs:
        [Us2[key], Vs2[key], Es2[key], Ss2[key], Ts2[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t2,depthlevel)
        
    mesh=ax.pcolormesh(Ss2['riversonly']-Ss['riversonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Difference in Salinity')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title('Times: ' +str(plot_ts[t2-18].hour) +' -' +str(plot_ts[t-18].hour) )

timesteps = [21,22,23]
fig, axs = plt.subplots(1, 3, figsize=(20,6),sharey=True)
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};
    Us2={}; Vs2={}; Es2 ={}; Ss2={}; Ts2={};
    t2=t+1

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    for key in runs:
        [Us2[key], Vs2[key], Es2[key], Ss2[key], Ts2[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t2,depthlevel)
        
    mesh=ax.pcolormesh(Ss2['riversonly']-Ss['riversonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Difference in Salinity')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title('Times: ' +str(plot_ts[t2-18].hour) +' -' +str(plot_ts[t-18].hour) )
    

Now adjust the surge field based by removing the elevation due to rivers.

In [13]:

fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20,21]
vmin=0; vmax=0.75
cmap='Reds'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly']-Es['riversonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Surge (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    
fig, axs = plt.subplots(1, 4, figsize=(20,6))
depthlevel=0
timesteps = [22,23,24,25]
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly']-Es['riversonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    cbar.set_label('Surge (m)')

** Comments**

wondering if more frequent output would be useful. We have 4 hour output now.
Removing the rivers didn't really make a difference, but at least now we have an idea of how much the rivers change the ssh.
Maybe it is best to remove the average anomaly over a certain period?

Why do the inlets have such a high surge?¶

Could the sea level pressure in these inlets be underestimated by CGRF? Simulations with all_forcing include the atmospheric pressure effect. If this is too low in the inlets, then the ssh will be elevated. The tideonly simulation does not include the atmospheric pressure effect on ssh. Sea level pressure could be underestimated if it is not actually sea level pressure (ie if it is surface pressure).

Let's check the CGRF data

In [14]:

CGRF = '/ocean/dlatorne/MEOPAR/CGRF/NEMO-atmos/slp_y2006m02d04.nc'

C = NC.Dataset(CGRF,'r')
nc_tools.show_variables(C)
nc_tools.show_variable_attrs(C)

[u'atmpres', u'nav_lat', u'nav_lon', u'time_counter']
<type 'netCDF4.Variable'>
float32 atmpres(time_counter, y, x)
    units: Pascals
    missing_value: 1e+20
    valid_min: 0.0
    valid_max: 1e+06
    long_name: surface atm. pressure
    short_name: atmpres
    online_operation: inst(only(x))
    axis: TYX
    scale_factor: 1.0
    add_offset: 0.0
    savelog10: 0.0
unlimited dimensions: time_counter
current shape = (24, 600, 801)
filling on, default _FillValue of 9.96920996839e+36 used

<type 'netCDF4.Variable'>
float32 nav_lat(y, x)
    units: degrees_north
    valid_min: -90.0
    valid_max: 90.0
    long_name: Latitude
    nav_model: Default grid
unlimited dimensions: 
current shape = (600, 801)
filling on, default _FillValue of 9.96920996839e+36 used

<type 'netCDF4.Variable'>
float32 nav_lon(y, x)
    units: degrees_east
    valid_min: 0.0
    valid_max: 360.0
    long_name: Longitude
    nav_model: Default grid
unlimited dimensions: 
current shape = (600, 801)
filling on, default _FillValue of 9.96920996839e+36 used

<type 'netCDF4.Variable'>
float64 time_counter(time_counter)
    calendar: gregorian
    title: Time
    long_name: Time axis
    units: hours since 2006-02-04 00:00:00
    time_origin: 2006-Feb-04 00:00:00
    valid_range: [ 0 23]
unlimited dimensions: time_counter
current shape = (24,)
filling on, default _FillValue of 9.96920996839e+36 used

In [15]:

press=C.variables['atmpres']
lat=C.variables['nav_lat']
lon=C.variables['nav_lon']
nc_tools.show_dimensions(C)

<type 'netCDF4.Dimension'> (unlimited): name = 'time_counter', size = 24

<type 'netCDF4.Dimension'>: name = 'y', size = 600

<type 'netCDF4.Dimension'>: name = 'x', size = 801

In [16]:

timestamps=[18,19,20]
fig, axs = plt.subplots(1, 3, figsize=(20,6),sharey=True)
vmin=80000; vmax=102000

for t,ax in zip(timestamps,axs):
    mesh=ax.pcolormesh(lon[:]-360,lat[:],press[t,:,:],vmin=vmin,vmax=vmax)
    viz_tools.plot_coastline(ax,fB,coords='map')
    cbar=fig.colorbar(mesh,ax=ax)
    ax.set_xlim([-126.5,-120])
    ax.set_ylim([47,52])
    ax.set_title('Feb 4, 2006 at ' + str(t) + ' oclock')
    cbar.set_label('CGRF Pressure Field [Pa]')
    
timestamps=[21,22,23]
fig, axs = plt.subplots(1, 3, figsize=(20,6),sharey=True)
vmin=80000; vmax=102000

for t,ax in zip(timestamps,axs):
    mesh=ax.pcolormesh(lon[:]-360,lat[:],press[t,:,:],vmin=vmin,vmax=vmax)
    viz_tools.plot_coastline(ax,fB,coords='map')
    cbar=fig.colorbar(mesh,ax=ax)
    ax.set_xlim([-126.5,-120])
    ax.set_ylim([47,52])
    ax.set_title('Feb 4, 2006 at ' + str(t) + ' oclock')
    cbar.set_label('CGRF Pressure Field [Pa]')

It looks to me like this is surface pressure. This would expain the higher anomalies in the inlets and also the higher anomaly in Washington just west of Puget Sound since the pressure drops there.

I guess I could still find an average anomaly state as a background.

Finding a background state¶

Ideas:

Could try an output time from my surge simulation, long after the surge has passed through.
State from a similar time of year in the spin ups

Lets try 1) first since it is most accessible at the moment. Look at both surgeonly and all_forcing minus tidesonly

In [17]:

plot_ts= nc_tools.timestamp(fTs['all_forcing'],[38,39,40,41])
fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [38,39,40,41]
vmin=-1; vmax=1
cmap='bwr'
print '1. surgeonly'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['surgeonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-38])

print '2. all_forcing minus tides only' 
fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-38])

1. surgeonly
2. all_forcing minus tides only

Ok, this looks like the typical background surge. So, I will use the anomaly from t=41 as a background.

In [18]:

t=41
Us={}; Vs={}; Es ={}; Ss={}; Ts={};

for key in {'all_forcing', 'tidesonly'}:
    [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
background = Es['all_forcing']-Es['tidesonly'];

plot_ts= nc_tools.timestamp(fTs['all_forcing'],[18,19,20,21,22,23, 24,25,26])
fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [18,19,20,21]
vmin=0; vmax=0.75
cmap='Reds'
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly']-background,cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Surge (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    
fig, axs = plt.subplots(1, 4, figsize=(20,6))
depthlevel=0
timesteps = [22,23,24,25]
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in runs:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing']-Es['tidesonly']-background,cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-18])
    cbar.set_label('Surge (m)')

Other Backgrounds - what is background for the run with no atmospheric forcing?¶

In [19]:

path='/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/feb2006/'    
runs = {'all_forcing','tidesonly', 'weather_only','ssh_only'}

fUs={}; fVs={}; fTs={};

for key in runs:
    fUs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_U.nc','r');
    fVs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_V.nc','r');
    fTs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_T.nc','r');


plot_ts= nc_tools.timestamp(fTs['all_forcing'],[38,39,40,41])
fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [38,39,40,41]
vmin=-1; vmax=1
cmap='bwr'

for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in {'ssh_only','tidesonly'}:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],
                                                                                 fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['ssh_only']-Es['tidesonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-38])
print 'ssh_only'

ssh_only

In [20]:

path='/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/feb2006/'    
runs = {'all_forcing','tidesonly', 'weather_only','ssh_only'}

fUs={}; fVs={}; fTs={};

for key in runs:
    fUs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_U.nc','r');
    fVs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_V.nc','r');
    fTs[key] = NC.Dataset(path + key +'/SalishSea_4h_20060201_20060207_grid_T.nc','r');


plot_ts= nc_tools.timestamp(fTs['all_forcing'],[38,39,40,41])
fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [38,39,40,41]
vmin=-1; vmax=1
cmap='bwr'

for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in {'weather_only','tidesonly'}:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],
                                                                                 fVs[key],fTs[key],t,depthlevel)
        
    mesh=ax.pcolormesh(Es['weather_only']-Es['tidesonly'],cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[t-38])
print 'weather_only'

weather_only

More Freqeuent Output - Dec 2006¶

In [21]:

#Load the data
path='/data/nsoontie/MEOPAR/SalishSea/results/storm-surges/tide_fix/dec2006/'

runs = {'all_forcing/1hour', 'tidesonly/1hour'}

fUs={}; fVs={}; fTs={};

for key in runs:
    fUs[key] = NC.Dataset(path + key +'/SalishSea_1h_20061214_20061215_grid_U.nc','r');
    fVs[key] = NC.Dataset(path + key +'/SalishSea_1h_20061214_20061215_grid_V.nc','r');
    fTs[key] = NC.Dataset(path + key +'/SalishSea_1h_20061214_20061215_grid_T.nc','r');
    
runs = {'all_forcing','tidesonly'}

for key in runs:
    fUs[key] = NC.Dataset(path + key +'/SalishSea_4h_20061211_20061217_grid_U.nc','r');
    fVs[key] = NC.Dataset(path + key +'/SalishSea_4h_20061211_20061217_grid_V.nc','r');
    fTs[key] = NC.Dataset(path + key +'/SalishSea_4h_20061211_20061217_grid_T.nc','r');

Surge starts half way through Dec 14.

In [22]:

ls=np.arange(18,42,1)
plot_ts= nc_tools.timestamp(fTs['all_forcing/1hour'],ls)
plot_ts[-1]

Out[22]:

<Arrow [2006-12-15T17:30:00+00:00]>

In [23]:

t=41
Us={}; Vs={}; Es ={}; Ss={}; Ts={};

for key in {'all_forcing', 'tidesonly'}:
    [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
background = Es['all_forcing']-Es['tidesonly'];

fig, axs = plt.subplots(4, 6, figsize=(30,25),sharey=True)
depthlevel=0
timesteps = ls
vmin=0; vmax=.75
cmap='jet'
count=0
for ax,t in zip (axs.flat,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in {'all_forcing/1hour','tidesonly/1hour'}:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],
                                                                                 t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing/1hour']-Es['tidesonly/1hour']-background,cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Anomaly  (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[count])
    count=count+1

Time increases from top left, moving right through each row.

It might be neat to make a 3D image of the SSH. Maybe even a movie
Highest anamolies occur in SoG.

Peak Surge - Dec 5, 7:30 -10:30¶

Note colorbar has changed

In [24]:

t=41
Us={}; Vs={}; Es ={}; Ss={}; Ts={};

for key in {'all_forcing', 'tidesonly'}:
    [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],t,depthlevel)
background = Es['all_forcing']-Es['tidesonly'];

fig, axs = plt.subplots(1, 4, figsize=(20,6),sharey=True)
depthlevel=0
timesteps = [31,32,33,34]
vmin=.35; vmax=.7
cmap='jet'
count=13
for ax,t in zip (axs,timesteps):
    Us={}; Vs={}; Es ={}; Ss={}; Ts={};

    for key in {'all_forcing/1hour','tidesonly/1hour'}:
        [Us[key], Vs[key], Es[key], Ss[key], Ts[key]] = stormtools.get_variables(fUs[key],fVs[key],fTs[key],
                                                                                 t,depthlevel)
        
    mesh=ax.pcolormesh(Es['all_forcing/1hour']-Es['tidesonly/1hour']-background,cmap=cmap,vmin=vmin,vmax=vmax)
    cbar = fig.colorbar(mesh, ax=ax)
    cbar.set_label('Anomaly  (m)')
    viz_tools.plot_coastline(ax,fB)
    ax.set_title(plot_ts[count])
    count=count+1

CGRF Winds¶

Look at the direction of the CGRF winds. Trying to determine why the mainland side has higher surges.

In [25]:

CGRF = '/ocean/dlatorne/MEOPAR/CGRF/NEMO-atmos/'

U= NC.Dataset(CGRF +'u10_y2006m12d15.nc','r')
V= NC.Dataset(CGRF +'v10_y2006m12d15.nc','r')

u=U.variables['u_wind'][:]
v=V.variables['v_wind'][:]

tims=U.variables['time_counter'][:]

In [26]:

#need to truncrate CGRF domain
xs=454; xe=472
ys=518; ye=531

latU=U.variables['nav_lat'][:]
lonU=U.variables['nav_lon'][:]
lonU=lonU-360

latV=V.variables['nav_lat'][:]
lonV=V.variables['nav_lon'][:]
lonV=lonV-360

print lonU[xs,ys], lonU[xe,ye]
print latU[xs,ys], latU[xe,ye]

print lonV[xs,ys], lonV[xe,ye]
print latV[xs,ys], latV[xe,ye]

-126.9 -121.05
46.35 51.75
-126.9 -121.05
46.35 51.75

In [27]:

lat=latU[xs:xe,ys:ye]
lon=lonU[xs:xe,ys:ye]
u=u[:,xs:xe,ys:ye]
v=v[:,xs:xe,ys:ye]

speeds = np.sqrt(np.square(u) + np.square(v))
print tims

[  0.   1.   2.   3.   4.   5.   6.   7.   8.   9.  10.  11.  12.  13.  14.
  15.  16.  17.  18.  19.  20.  21.  22.  23.]

In [31]:

ts=np.arange(3,15,1)
plot_ts= nc_tools.timestamp(U,ts)
fig, axs = plt.subplots(2, 6, figsize=(30,20 ))
count=0

for t,ax in zip(ts,axs.flat):
    viz_tools.plot_coastline(ax,fB,coords='map')
    quiver = ax.quiver(lon, lat, u[t,:,:], v[t,:,:],color='blue',scale=40,pivot='mid', width=0.01)
    ax.set_title(plot_ts[count])
    ax.set_xlim([-125,-122.5])
    ax.set_ylim([48,50])
    ax.quiverkey(quiver, -124.5, 48.9, 5, '5 m/s', coordinates='data', color='black', labelcolor='black')
    count=count+1

The time of first plot here corresponds to the 2nd row, 4th column of storm surge plot. Sorry for the mismatch in rows.

The max surge occurs around 9:00, at which time winds around Point Atkinson are mostly southerly. Then they change to westerly but after the max surge has occured.

In [28]: