<jupyter_start><jupyter_text>* This notebook was made to look at tidal heights at Steveston, Deas Island channel, New Westminster and Mission in the deepening river channel case.<jupyter_code>from __future__ import division
from cStringIO import StringIO
from salishsea_tools import (
    nc_tools,
    viz_tools,
    stormtools,
    tidetools,
)
from salishsea_tools.nowcast import figures
from salishsea_tools.nowcast import analyze
from matplotlib import pylab
from dateutil import tz
from matplotlib.backends import backend_agg as backend
from scipy import interpolate as interp
from matplotlib import animation
import scipy.io as sio
import datetime
import glob
import os
import arrow
import matplotlib.cm as cm
import matplotlib.dates as mdates
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import netCDF4 as nc
import numpy as np
import pandas as pd
import requests
import matplotlib.pyplot as plt
import datetime
import arrow
%matplotlib inline
sites = {
    'Sandheads': {
        'lat': 49.10,
        'lon': -123.30,
        'i' : 77, #426 start from 349 fortran start from 350
        'j' : 34, #293 start from 259 fortran start from 260
   
    },
    'Steveston Jetty':{
        'lat': 49.1057,
        'lon': -123.3032,
        'i' :  79,#428
        'j' :  33,#292
    },
    'Deas Island':{
        'lat' : 48.12,
        'lon' : -123.07,
        'i' : 65,#414,
        'j' : 69,#328
    },    
    'New Westminster': {
        'lat': 49.21,
        'lon': -122.89,
        'i': 69,#418,
        'j': 106,#365
        },

    'Mission': {
        'lat': 48.59,
        'lon': -123.2,
        'i':  151,#500,
        'j':108,#367
    },
    'Edge of domain':{
        'lat': 49.580574,
        'lon': -123.01498,
        'i': 151, ##500
        'j': 136, ##395
    },
    'nowcasts edge of domain':{
        'lat': 49.130043,
        'lon': -123.04199,
        'i': 65, ## 414
        'j': 75, ## 334
    },

        }
## location of Steveston from EC is:
#-123.185, 49.1242, so 13km away from
#our location for that in the model 
SITES = {
    '<NAME>': {
        'lat': 49.33,
        'lon': -123.25,
        'msl': 3.09,
        'stn_no': 7795,
        'extreme_ssh': 5.61},
    'Victoria': {
        'lat': 48.41,
        'lon': -123.36,
        'msl': 1.8810,
        'stn_no': 7120,
        'extreme_ssh': 3.76},
    '<NAME>': {
        'lat': 50.04,
        'lon': -125.24,
        'msl': 2.916,
        'stn_no': 8074,
        'extreme_ssh': 5.35},
    'Neah Bay': {
        'lat': 48.4,
        'lon': -124.6,
        'stn_no':  9443090},
    'Friday Harbor': {
        'lat': 48.55,
        'lon': -123.016667,
        'stn_no': 9449880},
    '<NAME>': {
        'lat': 48.866667,
        'lon': -122.766667,
        'stn_no': 9449424},
    'Sandheads': {
        'lat': 49.10,
        'lon': -123.30},
    'Tofino': {
        'stn_no': 8615},
    }
grid_TT = nc.Dataset('/data/jieliu/MEOPAR/river-treatment/14nor_NW/1hplume06160629gridT.nc')
t_orig, t_end, t_nemo = figures.get_model_time_variables(grid_TT)
def zoomin_period_tide(grid_T, stn_name, PST , MSL = 0):
    "Find out tide heights and time in PST or UTC \
     given a time series at stations we got our tidal predictions"
    time_shift = datetime.timedelta(hours=-7)# 7hrs in summer, 8hrs 
    #in winter
    t_orig, t_end, t_nemo = figures.get_model_time_variables(grid_T)
    ttide = figures.get_tides(stn_name)
    time = ttide.time + PST * time_shift
    ssh = ttide.pred_all + SITES[stn_name]['msl'] * MSL
    ssh_zoom = []
    time_zoom = []
    for i in np.arange(len(ttide)):
        if (time[i].year == t_orig.year) & (time[i].month == t_orig.month)\
        & (time[i].day >= t_orig.day) & (time[i].day <= t_end.day):
            ssh_zoom.append(ssh[i])
            time_zoom.append(time[i])
    return ssh_zoom, time_zoom, t_orig, t_end
def find_highlow_tide(grid_T, ssh_zoom, time_zoom, t_end, t_orig, PST ):
    "Determine high and low tide heights in PST or UTC at\
    a given day at stations we got our tidal predictions, unit of \
    return ssh_max, ssh_min is feet rather than meter"
    ssh_max = []
    ssh_min = []
    time_max = []
    time_min = []
    j = 0
    for i in np.arange(t_end.day - t_orig.day +1):
        ssh_max.append(max(ssh_zoom[j:j+24])/0.3048)
        ssh_min.append(min(ssh_zoom[j:j+24])/0.3048)
        time_max.append(time_zoom[ssh_zoom.index(max(ssh_zoom[j:j+24]))])  
        time_min.append(time_zoom[ssh_zoom.index(min(ssh_zoom[j:j+24]))]) 
        j = j +24
    return ssh_max, ssh_min, time_max, time_min
def print_out_diff_PA(ssh_max, ssh_min, meter):
    """This function was made to print out
        tidal range in meters/feet
        at Point Atkinson to determine
        which is spring/neap tide in the strait"""
    diff = []
    if meter ==0:
        for i, j in zip(ssh_max,ssh_min):
            dd =i-j
            diff.append(dd)
    elif meter ==1:
        for i, j in zip(ssh_max,ssh_min):
            dd =(i-j)* 0.3048
            diff.append(dd)
    return diff
ssh_zoom, time_zoom, t_orig, t_end = zoomin_period_tide(grid_TT,'Point Atkinson', 0, MSL = 0)
ssh_max,ssh_min,time_max,time_min=find_highlow_tide(grid_TT, ssh_zoom,time_zoom, t_end, t_orig, 0)
diff = print_out_diff_PA(ssh_max,ssh_min, meter =1)
diff ## tidal range at PA in  meters!!! June 24 has the smallest tidal range that I regard it as neap tide
##combine both max and min tidal anomaly for Point Atkinson
ssh_obs = []
ssh_obs.append(ssh_max)
ssh_obs.append(ssh_min)<jupyter_output><empty_output><jupyter_text>For new river channel model results(without smoothed bathymetry6)<jupyter_code>ssh_mod = grid_TT.variables['sossheig']
time_mod = grid_TT.variables['time_counter'][:]
ssh_mod.shape
t = nc_tools.timestamp(grid_TT, np.arange(time_mod.shape[0]))
for ind in range(len(t)):
    t[ind] = t[ind].datetime
tt = np.array(t) ## convert time_counter into datetime
def readoff_stn_highlow_waterlevel(t_end, t_orig, jnd, ind, ssh_mod, time_mod):
    """This function was used to read through high and low water\
        levels at various stns during a time period, the water levels\
        returned are in feet unit instead of meters"""
    ssh_mod_max = []
    ssh_mod_min = []
    time_mod_max = []
    time_mod_min = []
    j = 0
    for i in np.arange(t_end.day - t_orig.day +1):
        ssh_mod_max.append(max(ssh_mod[j:j+23, jnd, ind])/0.3048)
        ssh_mod_min.append(min(ssh_mod[j:j+23, jnd, ind])/0.3048)
        time_mod_max.append(time_mod[j:j+23][np.argmax(ssh_mod[j:j+23, jnd, ind])])  
        time_mod_min.append(time_mod[j:j+23][np.argmin(ssh_mod[j:j+23, jnd, ind])]) 
        ## need to convert time_mod_min in hour to date!!
        j = j +24
    return ssh_mod_max, ssh_mod_min, time_mod_max, time_mod_min  
def readoff_stn_highlow_waterlevel_new(t_end, t_orig, jnd, ind, ssh_mod, tt):
    """This function was used to read through high and low water\
        levels at various stns during a time period, the water levels\
        returned are in feet unit instead of meters"""
    ssh_mod_max = []
    ssh_mod_min = []
    time_mod_max = []
    time_mod_min = []
    j = 0
    for i in np.arange(t_end.day - t_orig.day +1):
        ssh_mod_max.append(max(ssh_mod[j:j+23, jnd, ind])/0.3048)
        ssh_mod_min.append(min(ssh_mod[j:j+23, jnd, ind])/0.3048)
        time_mod_max.append(tt[j:j+23][np.argmax(ssh_mod[j:j+23, jnd, ind])])  
        time_mod_min.append(tt[j:j+23][np.argmin(ssh_mod[j:j+23, jnd, ind])]) 
        ## need to convert time_mod_min in hour to date!!
        j = j +24
    return ssh_mod_max, ssh_mod_min, time_mod_max, time_mod_min 
def get_time_lag( ttmin, time_min, ttma, time_max ):
    """This function was made to calculate time lag
        of high and low water level between stns and 
        PA"""
    delta_min = []
    delta_max = []
    for t1, t2, t3, t4 in zip( ttmin, time_min,ttma, time_max):
        diff_min = (t1.hour- t2.hour) * 60 + (t1.minute - t2.minute)
        diff_max = (t3.hour- t4.hour) * 60 + (t3.minute - t4.minute)
        delta_min.append(diff_min)
        delta_max.append(diff_max)
    return delta_max, delta_min
def print_out_diff_wl(mima, mimi):
    """To print out difference of high and low
        water levels at each Fraser gauge for
        both models and obs from EC
        Input: mima & mimi are high & low water
        levels for stn.
        Output: array of their difference[unit feet]"""
    df = []
    for i, j in zip(mima, mimi):
        d = i -j
        df.append(d)
    return df
def plot_height_time_stn(ssh_model, ssh_obs, time_lag_mod, stn):
    """This function was made to plot water level height
    and time lag for Fraser gauges with PA tidal height """
    fig, axs = plt.subplots(1,2, figsize = (20,10))
    axs[0].scatter(ssh_model,ssh_obs,color='blue',s=15,edgecolor='none')
    axs[0].text(-2, 4, r'$Q_Hope = 200,000cfs$', fontsize=15)
    axs[0].grid()
    axs[0].set_xlim([-8, 14])
    axs[0].set_ylim([-10, 6])
    axs[0].set_xlabel('Water level at ' + stn + ' [feet]')
    axs[0].set_ylabel('Point Atkinson Height[feet]')
    axs[1].scatter(time_lag_mod,ssh_obs,color='blue',s=15,edgecolor='none')
    axs[1].grid()
    #ax.set_xlim([4.5, 9.5])
    #ax.set_ylim([-10, 6])
    axs[1].set_xlabel('Time lag at ' + stn + ' [Minute]')
    axs[1].set_ylabel('Point Atkinson Height[feet]')
    return fig<jupyter_output><empty_output><jupyter_text>* New Westminster<jupyter_code>## use another time counts[69,106]
tma, tmi, ttma, ttmin = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 69, 106,ssh_mod, tt)
##combine both max and min water level for NEW WESTMINSTER
ssh_model = []
ssh_model.append(tma)
ssh_model.append(tmi)
df = print_out_diff_wl(tma, tmi)
df
## average, max and min of the tidal amplitude at NW
ddf = np.array(df)
mean_ddf_NW = np.mean(ddf)* 0.3048
max_ddf_NW = np.max(ddf)* 0.3048
min_ddf_NW = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_NW))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_NW))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_NW))
delta_max, delta_min = get_time_lag(ttmin, time_min, ttma, time_max)<jupyter_output><empty_output><jupyter_text>* For spring tidal range in the strait about 4.0m, 50min delay at New Westminster.<jupyter_code>delta_max<jupyter_output><empty_output><jupyter_text>* For neap tidal range in the strait about 2m, 100min delay at New Westminster.<jupyter_code>delta_min<jupyter_output><empty_output><jupyter_text>* My question: * 1 different datum based in literture and PA prediction in my case?* 2 ours are hourly result, not accurate enough to calculate minutes.<jupyter_code>##combine both time lag for high and low water level at NEW WESTMINSTER
time_lag_mod = []
time_lag_mod.append(delta_max)
time_lag_mod.append(delta_min)
fig = plot_height_time_stn(ssh_model, ssh_obs, time_lag_mod, 'New Westminster')<jupyter_output><empty_output><jupyter_text>* Mission<jupyter_code>## grid cell for Mission is [151, 108]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 151, 108,ssh_mod, tt)
##combine both max and min water level for MISIION
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Mission
ddf = np.array(df)
mean_ddf_MI = np.mean(ddf)* 0.3048
max_ddf_MI = np.max(ddf)* 0.3048
min_ddf_MI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_MI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_MI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_MI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)<jupyter_output><empty_output><jupyter_text>* 190min delay for high water in the Ages model at Mission<jupyter_code>delta_max<jupyter_output><empty_output><jupyter_text>* 370min delay for low water in the Ages model at Mission<jupyter_code>delta_min
##combine both time lag for high and low water level at Mission
time_lag_mod = []
time_lag_mod.append(delta_max)
time_lag_mod.append(delta_min)
fig = plot_height_time_stn(ssh_model, ssh_obs, time_lag_mod, 'Mission')<jupyter_output><empty_output><jupyter_text>* Tide nearly no dissipation * Deas Island channel<jupyter_code>## grid cell for Deas Island is [65,69]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 65, 69,ssh_mod, tt)
##combine both max and min water level for DEAS ISLAND
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Deas Island channel for model
ddf = np.array(df)
mean_ddf_DI = np.mean(ddf)* 0.3048
max_ddf_DI = np.max(ddf)* 0.3048
min_ddf_DI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_DI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_DI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_DI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)<jupyter_output><empty_output><jupyter_text>* 10min delay at Deas Island for spring tide with tidal range of 5m at PA.<jupyter_code>delta_max<jupyter_output><empty_output><jupyter_text>* 70min delay at Deas Island for neap tides with tidal range of 0.5m at PA.<jupyter_code>delta_min
##combine both time lag for high and low water level at Deas Island
time_lag_mod = []
time_lag_mod.append(delta_max)
time_lag_mod.append(delta_min)
fig = plot_height_time_stn(ssh_model, ssh_obs, time_lag_mod, 'Deas Island')<jupyter_output><empty_output><jupyter_text>* Steveston Jetty<jupyter_code>## grid cell for Steveston Jetty is [79, 33]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 79, 33,ssh_mod, tt)
##combine both max and min water level for Steveston Jetty
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Steveston for model
ddf = np.array(df)
mean_ddf_ST = np.mean(ddf)* 0.3048
max_ddf_ST = np.max(ddf)* 0.3048
min_ddf_ST = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_ST))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_ST))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_ST))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)<jupyter_output><empty_output><jupyter_text>* 5min delay at Steveston for spring tide with tidal range of 5m at PA.<jupyter_code>delta_max<jupyter_output><empty_output><jupyter_text>* 25min delay at Steveston for spring tide with tidal range of 5m at PA.<jupyter_code>delta_min
##combine both time lag for high and low water level at Steveston Jetty
time_lag_mod = []
time_lag_mod.append(delta_max)
time_lag_mod.append(delta_min)
fig = plot_height_time_stn(ssh_model, ssh_obs, time_lag_mod, 'Steveston')<jupyter_output><empty_output><jupyter_text>* Looks good at Steveston* In this case, overdo the tide in the river, which seldom dissipate in the channel. Too high water level at Mission! * Edge of new river channel model domain<jupyter_code>## grid cell for edge of new river treatment is [151, 136]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 151, 136,ssh_mod, tt)
##combine both max and min water level for edge of model domain
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at edge of model domain 
ddf = np.array(df)
mean_ddf_ST = np.mean(ddf)* 0.3048
max_ddf_ST = np.max(ddf)* 0.3048
min_ddf_ST = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_ST))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_ST))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_ST))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min
##combine both time lag for high and low water level at Steveston Jetty
time_lag_mod = []
time_lag_mod.append(delta_max)
time_lag_mod.append(delta_min)
fig = plot_height_time_stn(ssh_model, ssh_obs, time_lag_mod, 'Edge of domain')<jupyter_output><empty_output><jupyter_text>New results in the Fraser River stations(with smoothed bathymetry6)<jupyter_code>grid_T = nc.Dataset('/data/jieliu/MEOPAR/river-treatment/14nor_NWSmoo/1hplume06160629smoogridT.nc')
t_orig, t_end, t_nemo = figures.get_model_time_variables(grid_T)
ssh_zoom, time_zoom, t_orig, t_end = zoomin_period_tide(grid_T,'Point Atkinson', 0, MSL = 0)
ssh_max,ssh_min,time_max,time_min=find_highlow_tide(grid_T, ssh_zoom,time_zoom, t_end, t_orig, 0)
diff = print_out_diff_PA(ssh_max,ssh_min, meter =1)
diff ## tidal range at PA in  meters!!! June 24 has the smallest tidal range that I regard it as neap tide
##combine both max and min tidal anomaly for Point Atkinson
ssh_obs = []
ssh_obs.append(ssh_max)
ssh_obs.append(ssh_min)
ssh_mod = grid_T.variables['sossheig']
time_mod = grid_T.variables['time_counter'][:]
t = nc_tools.timestamp(grid_T, np.arange(time_mod.shape[0]))
for ind in range(len(t)):
    t[ind] = t[ind].datetime
tt = np.array(t) ## convert time_counter into datetime<jupyter_output><empty_output><jupyter_text>* New Westminster<jupyter_code>## use another time counts[69,106]
tma, tmi, ttma, ttmin = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 69, 106,ssh_mod, tt)
##combine both max and min water level for NEW WESTMINSTER
ssh_model = []
ssh_model.append(tma)
ssh_model.append(tmi)
df = print_out_diff_wl(tma, tmi)
df
## average, max and min of the tidal amplitude at NW
ddf = np.array(df)
mean_ddf_NW = np.mean(ddf)* 0.3048
max_ddf_NW = np.max(ddf)* 0.3048
min_ddf_NW = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_NW))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_NW))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_NW))
delta_max, delta_min = get_time_lag(ttmin, time_min, ttma, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Mission<jupyter_code>## grid cell for Mission is [151, 108]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 151, 108,ssh_mod, tt)
##combine both max and min water level for MISIION
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Mission
ddf = np.array(df)
mean_ddf_MI = np.mean(ddf)* 0.3048
max_ddf_MI = np.max(ddf)* 0.3048
min_ddf_MI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_MI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_MI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_MI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Deas Island<jupyter_code>## grid cell for Deas Island is [65,69]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 65, 69,ssh_mod, tt)
##combine both max and min water level for Deas ISALND
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Deas Island channel for model
ddf = np.array(df)
mean_ddf_DI = np.mean(ddf)* 0.3048
max_ddf_DI = np.max(ddf)* 0.3048
min_ddf_DI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_DI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_DI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_DI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Steveston Jetty<jupyter_code>## grid cell for Steveston Jetty is [79, 33]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 79, 33,ssh_mod, tt)
##combine both max and min water level for Steveston Jetty
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Steveston for model
ddf = np.array(df)
mean_ddf_ST = np.mean(ddf)* 0.3048
max_ddf_ST = np.max(ddf)* 0.3048
min_ddf_ST = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_ST))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_ST))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_ST))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Compared tidal heights and time lag results of smoothed bathymetry6 with original bathymetry6, there is very little difference, which could be considered as the same. Therefore, almost no difference was caused by smoothed/unsmoothed bathymetry6.   For nowcast results * Deas Island<jupyter_code>def readoff_stn_highlow_waterlevel_nowcast(tf, to, ssh_daily, tt):
    """This function was used to read through high and low water\
        levels at various stns during a time period, the water levels\
        returned are in feet unit instead of meters in nowcasts results"""
    ssh_mod_max = []
    ssh_mod_min = []
    time_mod_max = []
    time_mod_min = []
    j = 0
    for i in np.arange(tf.day - to.day +1):
        ssh_mod_max.append(max(ssh_daily[j:j+23])/0.3048)
        ssh_mod_min.append(min(ssh_daily[j:j+23])/0.3048)
        time_mod_max.append(tt[j:j+23][np.argmax(ssh_daily[j:j+23])])  
        time_mod_min.append(tt[j:j+23][np.argmin(ssh_daily[j:j+23])]) 
        ## need to convert time_mod_min in hour to date!!
        j = j +24
    return ssh_mod_max, ssh_mod_min, time_mod_max, time_mod_min 
to = datetime.datetime(2015,6,16)
tf = datetime.datetime(2015,6,29)
files = analyze.get_filenames(to,tf,'1h','grid_T','/data/dlatorne/MEOPAR/SalishSea/nowcast/')
##combine files
ssh_daily, time_daily = analyze.combine_files(files, 'sossheig',\
                                    'None',sites['Deas Island']['i'] + 349,sites['Deas Island']['j']+ 259)
now_ssh_max, now_ssh_min, t_now_max, t_now_min = \
readoff_stn_highlow_waterlevel_nowcast(tf, to,ssh_daily, tt)
##combine both max and min water level for Deas Island of nowcast
ssh_model = []
ssh_model.append(now_ssh_max)
ssh_model.append(now_ssh_min)
df = print_out_diff_wl(now_ssh_max, now_ssh_min)
df
## average, max and min of the tidal amplitude at Deas Island of nowcast
ddf = np.array(df)
mean_ddf_DI_now = np.mean(ddf)* 0.3048
max_ddf_DI_now = np.max(ddf)* 0.3048
min_ddf_DI_now = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_DI_now))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_DI_now))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_DI_now))
delta_max, delta_min = get_time_lag(t_now_min, time_min, t_now_max, time_max)<jupyter_output><empty_output><jupyter_text>* Same time lag phase with new model results<jupyter_code>delta_max<jupyter_output><empty_output><jupyter_text>* More time lag than observations and new results, which should be around 40-50min for 3-4m tidal height at PA.<jupyter_code>delta_min<jupyter_output><empty_output><jupyter_text>* In summary, time lag of falling for low water level in Deas Island in nowcasts results is more than in the observations and new model results. Tidal range at Deas Island in the nowcasts is also too small which indicates tides dissipate too much and river not deep enough, so overall, the new river channel looks good for tide! * Edge of river channel for nowcasts<jupyter_code>to = datetime.datetime(2015,6,16)
tf = datetime.datetime(2015,6,29)
files = analyze.get_filenames(to,tf,'1h','grid_T','/data/dlatorne/MEOPAR/SalishSea/nowcast/')
##combine files
ssh_daily, time_daily = analyze.combine_files(files, 'sossheig',\
        'None',sites['nowcasts edge of domain']['i'] + 349,sites['nowcasts edge of domain']['j']+ 259)
now_ssh_max, now_ssh_min, t_now_max, t_now_min = \
readoff_stn_highlow_waterlevel_nowcast(tf, to,ssh_daily, tt)
##combine both max and min water level for edge of river channel of nowcast
ssh_model = []
ssh_model.append(now_ssh_max)
ssh_model.append(now_ssh_min)
df = print_out_diff_wl(now_ssh_max, now_ssh_min)
df
## average, max and min of the tidal amplitude at edge of river channel for nowcast
ddf = np.array(df)
mean_ddf_DI_now = np.mean(ddf)* 0.3048
max_ddf_DI_now = np.max(ddf)* 0.3048
min_ddf_DI_now = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_DI_now))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_DI_now))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_DI_now))
delta_max, delta_min = get_time_lag(t_now_min, time_min, t_now_max, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>Tidal heights inside Fraser River with shallowing of Mission to 3m<jupyter_code>grid_T = nc.Dataset('/data/jieliu/MEOPAR/river-treatment/14nor_Missionto3m/1hplume06160629Mission3mgridT.nc')
t_orig, t_end, t_nemo = figures.get_model_time_variables(grid_T)
ssh_zoom, time_zoom, t_orig, t_end = zoomin_period_tide(grid_T,'Point Atkinson', 0, MSL = 0)
ssh_max,ssh_min,time_max,time_min=find_highlow_tide(grid_T, ssh_zoom,time_zoom, t_end, t_orig, 0)
diff = print_out_diff_PA(ssh_max,ssh_min, meter =1)
diff
##combine both max and min tidal anomaly for Point Atkinson
ssh_obs = []
ssh_obs.append(ssh_max)
ssh_obs.append(ssh_min)
ssh_mod = grid_T.variables['sossheig']
time_mod = grid_T.variables['time_counter'][:]
t = nc_tools.timestamp(grid_T, np.arange(time_mod.shape[0]))
for ind in range(len(t)):
    t[ind] = t[ind].datetime
tt = np.array(t) ## convert time_counter into datetime<jupyter_output><empty_output><jupyter_text>* Deas Island<jupyter_code>## grid cell for Deas Island is [65,69]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 65, 69,ssh_mod, tt)
##combine both max and min water level for DEAS ISLAND
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Deas Island channel for model
ddf = np.array(df)
mean_ddf_DI = np.mean(ddf)* 0.3048
max_ddf_DI = np.max(ddf)* 0.3048
min_ddf_DI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_DI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_DI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_DI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min
##combine both time lag for high and low water level at NEW WESTMINSTER
time_lag_mod = []
time_lag_mod.append(delta_max)
time_lag_mod.append(delta_min)<jupyter_output><empty_output><jupyter_text>* New Westminster<jupyter_code>## use another time counts[69,106]
tma, tmi, ttma, ttmin = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 69, 106,ssh_mod, tt)
##combine both max and min water level for NEW WESTMINSTER
ssh_model = []
ssh_model.append(tma)
ssh_model.append(tmi)
df = print_out_diff_wl(tma, tmi)
df
## average, max and min of the tidal amplitude at NW
ddf = np.array(df)
mean_ddf_NW = np.mean(ddf)* 0.3048
max_ddf_NW = np.max(ddf)* 0.3048
min_ddf_NW = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_NW))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_NW))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_NW))
delta_max, delta_min = get_time_lag(ttmin, time_min, ttma, time_max)
delta_max
delta_min
##combine both time lag for high and low water level at NEW WESTMINSTER
time_lag_mod = []
time_lag_mod.append(delta_max)
time_lag_mod.append(delta_min)
fig = plot_height_time_stn(ssh_model, ssh_obs, time_lag_mod, 'New Westminster')<jupyter_output><empty_output><jupyter_text>* Mission<jupyter_code>## grid cell for Mission is [151, 108]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 151, 108,ssh_mod, tt)
##combine both max and min water level for MISIION
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Mission
ddf = np.array(df)
mean_ddf_MI = np.mean(ddf)* 0.3048
max_ddf_MI = np.max(ddf)* 0.3048
min_ddf_MI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_MI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_MI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_MI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Edge of model domain in the 3m Mission new case<jupyter_code>## grid cell for edge of new river treatment with Mission at 3m is [151, 136]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 151, 136,ssh_mod, tt)
tmima
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at edge of domain
ddf = np.array(df)
mean_ddf_MI = np.mean(ddf)* 0.3048
max_ddf_MI = np.max(ddf)* 0.3048
min_ddf_MI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_MI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_MI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_MI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Relatively higher tidal amplitude in Deas Island, NW than observation, relatively lower tidal amplitude in Mission. Looks okay. But the ssh pcolorplot in the strait is too low which I am a little bit worry that the tides in the Fraser is good but in the strait are not.? * It makes sense that water level at edge of model domain/Mission is high as in Ages paper, it argues that in freshet this should be the hydraulic grade line which is equal to water surface line and it should be around 3m higher than than at Sandheads or in the SoG. But our 3m bathymetry at Mission case makes it too high and if wanted to be accurate, tune back to around 6-8m at Mission, I guess, we will see how 8m at Mission works as follows.  Tidal heights inside Fraser River with shallowing of Mission to 8m<jupyter_code>grid_T = nc.Dataset('/data/jieliu/MEOPAR/river-treatment/14nor_Missionto8m/1hplume06160629Mission8mgridT.nc')
t_orig, t_end, t_nemo = figures.get_model_time_variables(grid_T)
ssh_zoom, time_zoom, t_orig, t_end = zoomin_period_tide(grid_T,'Point Atkinson', 0, MSL = 0)
ssh_max,ssh_min,time_max,time_min=find_highlow_tide(grid_T, ssh_zoom,time_zoom, t_end, t_orig, 0)
diff = print_out_diff_PA(ssh_max,ssh_min, meter =1)
diff
##combine both max and min tidal anomaly for Point Atkinson
ssh_obs = []
ssh_obs.append(ssh_max)
ssh_obs.append(ssh_min)
ssh_mod = grid_T.variables['sossheig']
time_mod = grid_T.variables['time_counter'][:]
t = nc_tools.timestamp(grid_T, np.arange(time_mod.shape[0]))
for ind in range(len(t)):
    t[ind] = t[ind].datetime
tt = np.array(t) ## convert time_counter into datetime<jupyter_output><empty_output><jupyter_text>* Steveston Jetty<jupyter_code>## grid cell for Steveston Jetty is [79, 33]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 79, 33,ssh_mod, tt)
##combine both max and min water level for STEESTON JETTY
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at STEVESTON JETTY for model
ddf = np.array(df)
mean_ddf_DI = np.mean(ddf)* 0.3048
max_ddf_DI = np.max(ddf)* 0.3048
min_ddf_DI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_DI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_DI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_DI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Deas Island<jupyter_code>## grid cell for Deas Island is [65,69]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 65, 69,ssh_mod, tt)
##combine both max and min water level for DEAS ISLAND
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Deas Island channel for model
ddf = np.array(df)
mean_ddf_DI = np.mean(ddf)* 0.3048
max_ddf_DI = np.max(ddf)* 0.3048
min_ddf_DI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_DI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_DI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_DI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* New Westminster<jupyter_code>## use another time counts[69,106]
tma, tmi, ttma, ttmin = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 69, 106,ssh_mod, tt)
##combine both max and min water level for NEW WESTMINSTER
ssh_model = []
ssh_model.append(tma)
ssh_model.append(tmi)
df = print_out_diff_wl(tma, tmi)
df
## average, max and min of the tidal amplitude at NW
ddf = np.array(df)
mean_ddf_NW = np.mean(ddf)* 0.3048
max_ddf_NW = np.max(ddf)* 0.3048
min_ddf_NW = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_NW))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_NW))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_NW))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Mission<jupyter_code>## grid cell for Mission is [151, 108]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 151, 108,ssh_mod, tt)
##combine both max and min water level for MISIION
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at Mission
ddf = np.array(df)
mean_ddf_MI = np.mean(ddf)* 0.3048
max_ddf_MI = np.max(ddf)* 0.3048
min_ddf_MI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_MI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_MI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_MI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Edge of model domain<jupyter_code>## grid cell for edge of new river treatment with Mission at 8m is [151, 136]
mima, mimi, tmima, tmimi = readoff_stn_highlow_waterlevel_new(t_end, t_orig, 151, 136,ssh_mod, tt)
##combine both max and min water level for edge of model domain
ssh_model = []
ssh_model.append(mima)
ssh_model.append(mimi)
df = print_out_diff_wl(mima, mimi)
df
## average, max and min of the tidal amplitude at edge of model domain with Mission of 8m
ddf = np.array(df)
mean_ddf_MI = np.mean(ddf)* 0.3048
max_ddf_MI = np.max(ddf)* 0.3048
min_ddf_MI = np.min(ddf)* 0.3048
print ('Max amplitude of model: {0:.4}'.format(max_ddf_MI))
print ('Min amplitude of model: {0:.4}'.format(min_ddf_MI))
print ('Mean amplitude of model: {0:.4}'.format(mean_ddf_MI))
delta_max, delta_min = get_time_lag(tmimi, time_min, tmima, time_max)
delta_max
delta_min<jupyter_output><empty_output><jupyter_text>* Shallowing Mission to 8 m decreased tidal height at Mission to about 1.2m, still around 0.9 m higher than observations, for Steveston Jetty is good, for Deas Island and New West, still 0.5m and 0.3m higher respectively. In a combination of Mission to 3 m where the tidal height at Mission is around 0.13 m which is 0.2 m lower than observations, the bathymetry after New West should be tunned around to 5m.<jupyter_code>ssh= grid_T.variables['sossheig'][:]
plt.pcolormesh(ssh[330, :,:],vmin=0, vmax=2)
plt.colorbar()
sshh= grid_TT.variables['sossheig'][:]
plt.pcolormesh(sshh[330, :,:],vmin=0, vmax=2)
plt.colorbar()
ssh[330:336, 151, 136]
sshh[330:336, 151, 136]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1211866307761328, 'fraction_numerical': 0.036395244504262, 'mean_word_length': 4.8911463495209775, 'pattern_counts': {'":': 0, '<': 153, '<?xml version=': 0, '>': 153, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '37', 'dir': 'jupyter-structured-clean-dedup', 'id': '252892', 'license': 'Apache-2.0', 'n_tokens_mistral': 14972, 'n_tokens_neox': 13712, 'n_words': 3551, 'path': 'jie/Tidal heights in deepening river channel.ipynb', 'repo_name': 'SalishSeaCast/analysis', 'size': '35686'}