Notebook to compare DFO tides with my tides.
DFO tides are from Anne's 07795const.wlev constituent file. The predictions are calculated with ttide using these constituents.
See /data/nsoontie/annes_tides/
In [1]:
import matplotlib.pyplot as plt
import pandas as pd
from salishsea_tools import stormtools
import datetime
from dateutil import tz
import numpy as np
%matplotlib inline
In [2]:
def date_parser(s):
unaware =datetime.datetime.strptime(s, '%Y-%m-%d %H:%M')
aware = unaware.replace(tzinfo=tz.tzutc())
return aware
In [3]:
DFO_file='/data/nsoontie/MEOPAR/analysis/Nancy/tides/Observed-DFO.csv'
DFO=pd.read_csv(DFO_file,parse_dates=[1],header=0,names=['station','time','tides','e1','e2'],date_parser=date_parser)
In [4]:
my_file1='/data/nsoontie/MEOPAR/analysis/Nancy/tides/Point Atkinson_t_tide_compare8_31-Dec-2013_02-Jan-2015.csv'
ttide1,msl1=stormtools.load_tidal_predictions(my_file1)
my_file1='/data/nsoontie/anne_tides/Point Atkinson_atide_compare8_31-Dec-2013_02-Dec-2015.csv'
ttide2,msl2=stormtools.load_tidal_predictions(my_file1)
In [5]:
fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.plot(DFO.time,DFO.tides,'-o',label='DFO')
ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all')
ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8')
ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes')
ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes')
ax.set_xlim(datetime.datetime(2014,12,7),datetime.datetime(2014,12,12))
ax.legend(loc=0)
Out[5]:
<matplotlib.legend.Legend at 0x7fdc7db91d50>
Zoom on Dec 8 lows
In [6]:
fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.plot(DFO.time,DFO.tides,'-o',label='DFO')
ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all')
ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8')
ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes')
ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes')
ax.set_xlim(datetime.datetime(2014,12,8),datetime.datetime(2014,12,9))
ax.set_ylim([0,1])
ax.legend(loc=0)
ax.grid()
Predictions with Anne's constituents line up with DFO exactly.
Zoom on Dec 8 highs
In [7]:
fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.plot(DFO.time,DFO.tides,'-o',label='DFO')
ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all')
ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8')
ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes')
ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes')
ax.set_xlim(datetime.datetime(2014,12,8),datetime.datetime(2014,12,9))
ax.set_ylim([4,5])
ax.legend(loc=0)
ax.grid()
Zoom on Dec 8 middles
In [8]:
fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.plot(DFO.time,DFO.tides,'-o',label='DFO')
ax.plot(ttide1.time,ttide1.pred_all +msl1,'-o',label='pred_all')
ax.plot(ttide1.time,ttide1.pred_8+msl1,label='pred_8')
ax.plot(ttide2.time,ttide2.pred_all +msl2,'-o',label='pred_all annes')
ax.plot(ttide2.time,ttide2.pred_8+msl2,label='pred_8 annes')
ax.set_xlim(datetime.datetime(2014,12,8),datetime.datetime(2014,12,9))
ax.set_ylim([3,4])
ax.legend(loc=0)
ax.grid()
Looks like we are shifted in our tidal predictions. Eyeballing, we are 10-15cm low (pred_all)
How do means, maxes and mins compare for Dec 8-12? Look at Annes' and DFO
In [9]:
t1=datetime.datetime(2014,12,8); t2=datetime.datetime(2014,12,12)
t1 = t1.replace(tzinfo=tz.tzutc())
t2 = t2.replace(tzinfo=tz.tzutc())
subDFO=np.array(DFO.tides)
timDFO=np.array(DFO.time)
indices = np.where(np.logical_and(timDFO >= t1, timDFO <= t2))
subDFO=subDFO[indices];
timDFO=timDFO[indices];
meanDFO=np.nanmean(subDFO); maxDFO=np.nanmax(subDFO); minDFO=np.nanmin(subDFO)
submine=np.array(ttide2.pred_all);
timmine=np.array(ttide2.time);
indices = np.where(np.logical_and(timmine >= t1, timmine <= t2))
submine=submine[indices] +msl2;
timmine=timmine[indices];
meanmine=np.nanmean(submine); maxmine=np.nanmax(submine); minmine=np.nanmin(submine)
In [10]:
print 'DFO: Mean {}, Max {}, Min {}'.format(meanDFO,maxDFO,minDFO)
print 'Mine: Mean {}, Max {}, Min {}'.format(meanmine,maxmine,minmine)
DFO: Mean 3.13881439722, Max 4.77, Min 0.681 Mine: Mean 3.14256752577, Max 4.762454, Min 0.683941
In [11]:
print 'Differences: Mean {}, Max {}, Min {}'.format(meanDFO-meanmine,maxDFO-maxmine,minDFO-minmine)
Differences: Mean -0.00375312854856, Max 0.007546, Min -0.002941
These differences are so small and are likely due to differences in timing.
Effect on tidal corrections¶
Look to see the effect on the model corrections in the storm surge preditions
In [12]:
from salishsea_tools.nowcast import figures
import datetime
from glob import glob
import os
import netCDF4 as nc
import scipy.io as sio
%matplotlib inline
In [13]:
def results_dataset(period, grid, results_dir):
"""Return the results dataset for period (e.g. 1h or 1d)
and grid (e.g. grid_T, grid_U) from results_dir.
"""
filename_pattern = 'SalishSea_{period}_*_{grid}.nc'
filepaths = glob(os.path.join(results_dir, filename_pattern.format(period=period, grid=grid)))
return nc.Dataset(filepaths[0])
In [14]:
run_date = datetime.datetime(2015,2,8)
# Results dataset location
results_home = '/data/dlatorne/MEOPAR/SalishSea/nowcast/'
results_dir = os.path.join(results_home, run_date.strftime('%d%b%y').lower())
# model winds
model_path = '/ocean/sallen/allen/research/MEOPAR/Operational/'
coastline = sio.loadmat('/ocean/rich/more/mmapbase/bcgeo/PNW.mat')
In [15]:
grid_T_hr = results_dataset('1h', 'grid_T', results_dir)
bathy = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
In [16]:
from salishsea_tools.nowcast import analyze
from salishsea_tools import tidetools, stormtools
import matplotlib.pyplot as plt
from dateutil import tz
In [17]:
t_orig=datetime.datetime(2014,12,1); t_final=datetime.datetime(2015,1,31)
path='/data/dlatorne/MEOPAR/SalishSea/nowcast/'
files_all=analyze.get_filenames(t_orig,t_final,'1h','grid_T','nowcast')
In [18]:
def get_new_tides(name):
if name == 'Point Atkinson':
path = ('/data/nsoontie/anne_tides/')
filename ='{}_atide_compare8_31-Dec-2013_02-Dec-2015.csv'.format(name)
else:
# Tide file covers 2014 and 2015. Harmonics from a 2013 time series.
path = (
'/data/nsoontie/MEOPAR/tools/SalishSeaTools/salishsea_tools/'
'nowcast/tidal_predictions/')
filename = '{}_t_tide_compare8_31-Dec-2013_02-Dec-2015.csv'.format(name)
tfile = os.path.join(path, filename)
ttide, msl = stormtools.load_tidal_predictions(tfile)
return ttide
In [19]:
lats,lons=figures.station_coords()
b,X,Y= tidetools.get_bathy_data(bathy)
name='Point Atkinson'
j,i=tidetools.find_closest_model_point(lons[name],lats[name],X,Y,b,)
ssh, time=analyze.combine_files(files_all,'sossheig','None',j,i)
ttide_new=get_new_tides(name)
ttide_old = figures.get_tides(name)
In [20]:
msl=3.09
sdt = t_orig.replace(tzinfo=tz.tzutc())
edt = t_final +datetime.timedelta(days=1)
edt=edt.replace(tzinfo=tz.tzutc())
ssh_corr_new = stormtools.correct_model(ssh,ttide_new,sdt,edt)
ssh_corr_old = stormtools.correct_model(ssh,ttide_old,sdt,edt)
obs=stormtools.load_observations((t_orig-datetime.timedelta(days=1)).strftime('%d-%b-%Y'),
(t_final+datetime.timedelta(days=1)).strftime('%d-%b-%Y'),'PointAtkinson')
#obs=figures.load_PA_observations()
In [21]:
fig,ax=plt.subplots(1,1, figsize=(15,5))
ax.plot(obs.time,obs.slev,'g',lw=1.5,label='observations')
ax.plot(time,ssh+msl,'k--',lw=1.5,label='model')
ax.plot(time,ssh_corr_new+msl,'b',lw=1.5,label='corrected model - anns tides')
ax.plot(time,ssh_corr_old+msl,'r',lw=1.5,label='corrected model - my tides')
plt.legend(loc=0)
td=datetime.timedelta(days=17)
win=datetime.timedelta(days=5)
ax.set_xlim([t_orig+td,t_orig+td+win])
ax.grid()
In [22]:
def interp_to_model_time(time_model,varp,tp):
""" interpolate to model time
"""
# Strategy: convert times to seconds past a reference value.
# Use this as the independent variable in interpolation.
#Set epoc (reference) time.
epoc=time_model[0];
# Determine tp times wrt epc
tp_wrt_epoc=[]
for t in tp:
tp_wrt_epoc.append((t-epoc).total_seconds())
# Interpolate observations to model times
varp_interp=[]
for t in time_model:
mod_wrt_epoc=(t-epoc).total_seconds();
varp_interp.append(np.interp(mod_wrt_epoc,tp_wrt_epoc,varp))
return varp_interp
In [23]:
obs_interp=interp_to_model_time(time,obs.slev,obs.time)
In [24]:
fig,ax=plt.subplots(1,1, figsize=(15,5))
ax.plot(time,obs_interp,'g',lw=1.5,label='observations')
ax.plot(time,ssh+msl,'k--',lw=1.5,label='model')
ax.plot(time,ssh_corr_new+msl,'b',lw=1.5,label='corrected model - anns tides')
ax.plot(time,ssh_corr_old+msl,'r',lw=1.5,label='corrected model - my tides')
plt.legend(loc=0)
td=datetime.timedelta(days=17)
win=datetime.timedelta(days=5)
#ax.set_xlim([t_orig+td,t_orig+td+win])
ax.grid()
Look at errors
In [25]:
errors={}
errors['new']=ssh_corr_new+msl-obs_interp
errors['old'] = ssh_corr_old +msl - obs_interp
errors['none'] = ssh +msl -obs_interp
for n in ['new','old','none']:
print n, 'mean error:', np.mean(errors[n]), 'max error:', np.max(errors[n]), 'min error:', np.min(errors[n])
new mean error: -0.00621937923476 max error: 0.288683449532 min error: -0.357990750374 old mean error: -0.109088827823 max error: 0.18190589049 min error: -0.431898750374 none mean error: -0.0636528634417 max error: 0.226186265945 min error: -0.364891257286
Compare with Foreman's constituents¶
In [26]:
filename = '/data/nsoontie/MEOPAR/analysis/compare_tides/obs_tidal_wlev_const_all.csv'
harm_obs = pd.read_csv(filename,sep=';',header=0)
harm_obs = harm_obs.rename(columns={'Site': 'site', 'Lat': 'lat', 'Lon': 'lon',
'M2 amp': 'M2_amp', 'M2 phase (deg UT)': 'M2_pha',
'K1 amp': 'K1_amp', 'K1 phase (deg UT)': 'K1_pha'})
filename = '/data/nsoontie/MEOPAR/analysis/Idalia/other_constituents.csv'
harm_other = pd.read_csv(filename,sep=',',header=0)
harm_other = harm_other.rename(columns={'Site': 'site', 'Lat': 'lat', 'Lon': 'lon',
'O1 amp': 'O1_amp', 'O1 phase (deg UT)': 'O1_pha',
'P1 amp': 'P1_amp', 'P1 phase (deg UT)': 'P1_pha',
'Q1 amp': 'Q1_amp', 'Q1 phase (deg UT)': 'Q1_pha',
'S2 amp': 'S2_amp', 'S2 phase (deg UT)': 'S2_pha',
'N2 amp': 'N2_amp', 'N2 phase (deg UT)': 'N2_pha',
'K2 amp': 'K2_amp', 'K2 phase (deg UT)': 'K2_pha'})
In [27]:
s='Point Atkinson'
harm_other[harm_other.site==s]
Out[27]:
| site | lat | lon | O1_amp | O1_pha | P1_amp | P1_pha | Q1_amp | Q1_pha | S2_amp | S2_pha | N2_amp | N2_pha | K2_amp | K2_pha | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 15 | Point Atkinson | 49.334 | -123.25 | 48.3 | 263.2 | 26.8 | 283.1 | 7.7 | 258.8 | 22.9 | 59.9 | 18.4 | 2.9 | 6.2 | 59.9 |
In [28]:
ampsF={}; phasF={}
for t in ['K1','M2']:
astr=t + '_amp'
pstr=t + '_pha'
ampsF[t]= harm_obs[harm_obs.site==s][astr].values
phasF[t]= harm_obs[harm_obs.site==s][pstr].values
tides=['O1','P1','Q1','S2','N2','K2']
for t in tides:
astr=t + '_amp'
pstr=t + '_pha'
ampsF[t]= harm_other[harm_other.site==s][astr].values
phasF[t]= harm_other[harm_other.site==s][pstr].values
In [29]:
annes='/data/nsoontie/anne_tides/07795const.wlev'
atides=pd.read_table(annes,skiprows=24,names=['per','amp','pha','a','b','c','d','e','f','g'],delim_whitespace=True)
In [30]:
ampsA={};phasA={}
tides=['K1','O1','P1','Q1','M2','S2','N2','K2']
for t in tides:
pha=atides.ix[t].pha +8.0*360/atides.ix[t].per #convert to UTC
ampsA[t]=atides.ix[t].amp*100
phasA[t]= pha -360*(pha>360)
In [31]:
fig,axs=plt.subplots(2,1)
for t in tides:
axs[0].plot(ampsF[t],ampsA[t],'o',label=t)
axs[1].plot(phasF[t],phasA[t],'o',label=t)
axs[0].plot([0,100],[0,100],'r')
axs[1].plot([0,360],[0,360],'r')
axs[0].set_xlabel('Foreman Amp')
axs[0].set_ylabel('Annes Amp')
axs[1].set_xlabel('Foreman Phase')
axs[1].set_ylabel('Annes Phase')
axs[0].legend(loc=0)
Out[31]:
<matplotlib.legend.Legend at 0x7fdc74068450>
In [32]:
print 'Constitents, ampA, ampF, phaA, phaF'
for t in tides:
print t, ampsA[t], ampsF[t], phasA[t], phasF[t]
Constitents, ampA, ampF, phaA, phaF K1 86.22 [ 86.2] 286.240910003 [ 286.1] O1 47.58 [ 48.3] 263.365760874 [ 263.2] P1 27.18 [ 26.8] 284.740905011 [ 283.1] Q1 7.77 [ 7.7] 257.020710138 [ 258.8] M2 91.91 [ 91.8] 31.0053892601 [ 31.2] S2 23.23 [ 22.9] 59.48 [ 59.9] N2 19.42 [ 18.4] 3.15409543372 [ 2.9] K2 6.22 [ 6.2] 59.511820005 [ 59.9]
There is a reasonable match between Anne's and Foreman's for the 8 constituents.
In [32]: