This notebook creates monthly forcing files for the sea surface height (hourly frequency) at Tofino. We want sea surface height forcing separated from the tides. We will take the observed sea surface height and remove the tidal predictions.

I will begin by looking at Jan 1, 2006 through to to Dec. 31, 2006. There was another great storm surge in Dec. 2006. Since data NEMO in UTC but the observations come in PST I will download data rom Dec 31, 2005 to Jan 2, 2007.

Storm surge dates: Nov 15, 2006, Dec 14, 2006

Observations: DFO

Tidal predictions: t_tide

In [1]:
%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import netCDF4 as NC
import numpy as np
import datetime
import requests

from salishsea_tools import nc_tools
from salishsea_tools import tidetools
import os
import csv

DFO observed water level

Begin by collecting water level data from the DFO.

In order to generate forcing files for a different set of dates, modify "start" and "end". The time period needs to be at least one month long. Due to local time vs UTC, "start" needs to be one day before the month of interest and "end" needs to be two days after the desired end date.

In [2]:
#Tofino mmsl over 1982-2001 (in metres)
msl = 2.103
In [3]:
yr=2009

statID = '8615'#Tofino station ID number
start = '31-Dec-' + str(yr-1) #start date 
end = '02-Jan-' + str(yr+1) #end date

#tidetools.get_dfo_wlev(statID,start,end)
#observations and tidal prediction stored in storm-surge/Revisions/tides/forcing
filename = '/data/nsoontie/MEOPAR/storm-surge/Revisions/tides/forcing/wlev_'+str(statID) + '_' + start +'_' +end +'.csv'

import datetime
import pandas as pd
import pytz

def dateParserMeasured(s):
    #convert the string to a datetime object
    unaware = datetime.datetime.strptime(s, "%Y/%m/%d %H:%M")
    #add in the local time zone (Canada/Pacific)
    aware = unaware.replace(tzinfo=pytz.timezone('Canada/Pacific'))
    #convert to UTC
    return aware.astimezone(pytz.timezone('UTC'))
    
wlev_meas = pd.read_csv(filename,skiprows=7,parse_dates=[0],date_parser=dateParserMeasured)
wlev_meas = wlev_meas.rename(columns={'Obs_date': 'time', 'SLEV(metres)': 'slev'})

print(wlev_meas.time[0].tzinfo)

#the average sea surface height over this time. 
sshavg = wlev_meas.slev.mean()
print sshavg

plt.figure(figsize=(14,6))
pylab.subplot(1,2,1)
plt.plot(wlev_meas.time,wlev_meas.slev)
plt.title('Tofino Observations' + start+' to ' + end)
plt.plot(wlev_meas.time,sshavg*np.ones(wlev_meas.time.size),'--k')
plt.xlabel('Time (UTC)')
plt.ylabel('Water Level Height (m)')

pylab.subplot(1,2,2)
plt.plot(wlev_meas.time,wlev_meas.slev-sshavg)
plt.title('Tofino Observations, mean removed: ' + start+' to ' + end)
plt.xlabel('Time (UTC)')
plt.ylabel('Water Level Anomoly (m)')
UTC
2.07770095368
Out[3]:
<matplotlib.text.Text at 0x7f8620d4d350>

Additionally, download measurements from the year prior. These measurements will be used to calculate harmonics and tidal predictions for the year of interest.

Predictions from t_tide

Run t_tide on the time series to get the tidal predictions.

This is a little bit clunky, but I have modified Kate's get_tidal_anomaly.m matlab script to use t_tide. This script saves the predcitions in location_t_tide_compare8_dates_snr2_filter.csv file. It also saves the harmoncis in location_harmonics_dates_filter.csv. Before this part of the notebook is executed, we must go to matlab to generate these files:

eg.: get_ttide_8_filter('wlev_8615_31-Dec-2006_02-Jan-2007.csv', 'Tofino', '31-Dec-2005', '02-Jan-2007');

The first string is a filename that stores the water level measurements. A harmonic analysis is performed on this time series. The second string indicates the location and the last two string indicate the beginning and end dates of the tidal predications.

The matlab script get_ttide_8_filter.m is stored in the storm-surge/Revisions/tides/ repository. It is recommended to add this directory to your path and run the script in the directory where you would like the predictions to be stored. You also need t_tide added you your path.

This has already been done for 2006, 2009 and 2012. The results are stored in storm-surge/Revisions/tides/forcing/

For now, load the data produced by the matlab script.

In [4]:
filename = '/data/nsoontie/MEOPAR/storm-surge/Revisions/tides/forcing/Tofino_t_tide_compare8_' + start+'_'+end+'_snr2_filter.csv'


import datetime
import pandas as pd
import pytz

def dateParserMeasured2(s):
    #convert the string to a datetime object
    unaware = datetime.datetime.strptime(s, "%d-%b-%Y %H:%M:%S ")
    #add in the local time zone (Canada/Pacific)
    aware = unaware.replace(tzinfo=pytz.timezone('Canada/Pacific'))
    #convert to UTC
    return aware.astimezone(pytz.timezone('UTC'))

#read msl
line_number = 1     
with open(filename, 'rb') as f:
    mycsv = csv.reader(f); mycsv = list(mycsv)
    
ttide = pd.read_csv(filename,skiprows=3,parse_dates=[0],date_parser=dateParserMeasured2)
ttide = ttide.rename(columns={'Time_Local ': 'time', ' pred_8 ': 'pred_8', ' pred_all ': 'pred_all'})

ssanomaly = np.zeros(len(wlev_meas.time))

for i in np.arange(0,len(wlev_meas.time)):
    #check that there is a corresponding time 
    #if any(wlev_pred.time == wlev_meas.time[i]):
    ssanomaly[i] =(wlev_meas.slev[i] - (ttide.pred_all[ttide.time==wlev_meas.time[i]]+msl))
    if not(ssanomaly[i]):
        ssanomaly[i]=float('Nan')
        
In [5]:
    print 'Comparison between ttide predictions and measurements'
    
    plt.figure(figsize=(13,6))
    plt.subplot(2,1,1)
    plt.plot(wlev_meas.time,wlev_meas.slev,'b',label='measured')
    plt.plot(ttide.time,ttide.pred_all+msl,'m',label='predicted')
    #plt.xlim(datetime.datetime(req_year,1,1,0,0,0),datetime.datetime(req_year,12,31,23,59,59))
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
    plt.xlabel('Time (UTC)')
    plt.ylabel('Elevation (m)')
    
    plt.subplot(2,1,2)
    plt.plot(wlev_meas.time,ssanomaly,'g',label='anomaly (meas-pred)')
    plt.ylim(-1,1)
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
    plt.xlabel('Time (UTC)')
    plt.ylabel('Elevation (m)')
    plt.plot(wlev_meas.time,0*np.ones(wlev_meas.time.size),'--r')
    
    ssarray = np.array(ssanomaly)
    mean_anom= ssarray.mean()
    print 'Mean anomaly: ' + str(mean_anom)
Comparison between ttide predictions and measurements
Mean anomaly: -0.0252720551771

Save as a NetCDF

Follow Susan's SSH notebook for this. Make a new netcdf file for each month. I will save the ssanonmaly at each point along the western boundary.

Get some preliminary data first: bathymetry, indices of the edges.

In [6]:
startj = 384
endj = 471
lengthj = endj-startj
r = 1

fB = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
lat = fB.variables['nav_lat'][:]
lon = fB.variables['nav_lon'][:]
fB.close()
print lat.shape
(898, 398)

A function for saving the netcdf file.

In [7]:
def prepare_netcdf(ssh_file, count, ssh,month,year):
    # dataset attributes
    nc_tools.init_dataset_attrs(
    ssh_file, 
    title='Juan de Fuca SSH hourly values', 
    notebook_name='SSH_Tofino', 
    nc_filepath='/data/nsoontie/MEOPAR/NEMO-forcing/open_boundaries/west/ssh/' + filename,
    comment='SSH from Tofino: measured - tides (from t_tide)')

    
    #dimensions
    ssh_file.createDimension('xbT', lengthj*r)
    ssh_file.createDimension('yb', 1)
    ssh_file.createDimension('time_counter', None)
    # variables
    # time_counter
    time_counter = ssh_file.createVariable('time_counter', 'float32', ('time_counter'))
    time_counter.long_name = 'Time axis'
    time_counter.axis = 'T'
    time_counter.units = 'hour since 00:00:00 on 01/' + str(month) +'/'+ str(year)
    # nav_lat and nav_lon
    nav_lat = ssh_file.createVariable('nav_lat','float32',('yb','xbT'))
    nav_lat.long_name = 'Latitude'
    nav_lat.units = 'degrees_north'
    nav_lon = ssh_file.createVariable('nav_lon','float32',('yb','xbT'))
    nav_lon.long_name = 'Longitude'
    nav_lon.units = 'degrees_east'
    # ssh
    sossheig = ssh_file.createVariable('sossheig', 'float32', 
                               ('time_counter','yb','xbT'), zlib=True)
    sossheig.units = 'm'
    sossheig.long_name = 'Sea surface height'   
    sossheig.coordinates = 'nav_lon nav_lat time_counter'
    sossheig.grid = 'SalishSea2'
    # vobtcrtx, vobtcrty
    vobtcrtx = ssh_file.createVariable('vobtcrtx', 'float32',
                                   ('time_counter','yb','xbT'), zlib=True)
    vobtcrtx.units = 'm/s'
    vobtcrtx.long_name = 'Barotropic U Velocity- ZEROD'   
    vobtcrtx.grid = 'SalishSea2'
    vobtcrty = ssh_file.createVariable('vobtcrty', 'float32',
                                   ('time_counter','yb','xbT'), zlib=True)
    vobtcrty.units = 'm/s'
    vobtcrty.long_name = 'Barotropic V Velocity- ZEROD'   
    vobtcrty.grid = 'SalishSea2'
    # nbidta, ndjdta, ndrdta
    nbidta = ssh_file.createVariable('nbidta', 'int32' , ('yb','xbT'), zlib=True)
    nbidta.long_name = 'i grid position'
    nbidta.units = 1
    nbjdta = ssh_file.createVariable('nbjdta', 'int32' , ('yb','xbT'), zlib=True)
    nbjdta.long_name = 'j grid position'
    nbjdta.units = 1
    nbrdta = ssh_file.createVariable('nbrdta', 'int32' , ('yb','xbT'), zlib=True)
    nbrdta.long_name = 'position from boundary'
    nbrdta.units = 1
    
    for ir in range(0,r):
        nav_lat[0,ir*lengthj:(ir+1)*lengthj] = lat[startj:endj,ir]
        nav_lon[0,ir*lengthj:(ir+1)*lengthj] = lon[startj:endj,ir]
        nbidta[0,ir*lengthj:(ir+1)*lengthj] = ir
        nbjdta[0,ir*lengthj:(ir+1)*lengthj] = range(startj,endj)
        nbrdta[0,ir*lengthj:(ir+1)*lengthj] = ir
        
    for ib in range(0,lengthj*r):
        sossheig[0:count,0,ib] = ssh[0:count]
        time_counter[0:count] = range(1,count+1)
        vobtcrtx[0:count,0,ib] = 0*np.ones(count)
        vobtcrty[0:count,0,ib] = 0*np.ones(count)

Loop through the ssanomaly and time arrays to create the netcdf files. Save them month by month.

This is really clunky for now. I need the netcdf files in UTC, but all of the data files are in local time. So I will just ignore/delete the first and last month's files since these are incomplete. There is definitely room for improvement with this part of the code.

In [8]:
#setting up the loop through data

#length of data:
l=len(ssanomaly)
hour_start = wlev_meas.time[0].hour;  yr_start= wlev_meas.time[0].year; month_start =wlev_meas.time[0].month; 
yr_end= wlev_meas.time[l-1].year;  month_end =wlev_meas.time[l-1].month; 
print 'Date range: ' +str(wlev_meas.time[0]) + ' to ' + str(wlev_meas.time[l-1]) + ' UTC'

#files that will be deleted if hour does not begin at zero. 
if month_start<10:
    filename = 'ssh_y' +str(yr_start)+'m0'+str(month_start)+'.nc'
else:
    filename = 'ssh_y' +str(yr_start)+'m'+str(month_start)+'.nc'
file1=filename
if month_end<10:
    file2 = 'ssh_y' +str(yr_end)+'m0'+str(month_end)+'.nc'
else:
    file2 = 'ssh_y' +str(yr_end)+'m'+str(month_end)+'.nc'

m=month_start; yr=yr_start; count =0; ssh=[];

for ii in range(0, l):
    yr=wlev_meas.time[ii].year
    month = wlev_meas.time[ii].month
    #If we entered a new month, save all the data.
    if month != m:
        #save the month's data:
        ssh_file = NC.Dataset(filename, 'w', zlib=True)
        print filename
        print count       
        #save data in netcdf    
        prepare_netcdf(ssh_file,count,ssh, m, yr)
        #Set up the next month's file
        m=month;
        if m<10:
            filename = 'ssh_y' +str(yr)+'m0'+str(m)+'.nc'
        else:
            filename = 'ssh_y' +str(yr)+'m'+str(m)+'.nc'
        #clear the ssh and count
        ssh = []; count=0
        ssh_file.close()
 
    ssh.append(ssanomaly[ii])
    count=count+1
    
    #If we are at the last item, save all the data.
    if ii == l-1:
        ssh_file = NC.Dataset(filename, 'w', zlib=True)
        print filename
        print count       
        prepare_netcdf(ssh_file,count,ssh,month,yr)
        #clear the ssh and count
        ssh = []; count=0
        ssh_file.close()
    
#If the hour is not zero then we have incomplete monthly data. Delete the incomplete files. 
#Note: this is not a proper check for incomplete data. For example, if the day is not the first day of the month...
if hour_start != 0:
    print '**** Removing ' + file1 + ' and ' + file2
    os.remove(file1); os.remove(file2)
    
Date range: 2008-12-31 08:12:00+00:00 to 2010-01-02 08:12:00+00:00 UTC
ssh_y2008m12.nc
16
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:52] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m01.nc
743
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:53] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m02.nc
672
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:53] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m03.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:54] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m04.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:55] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m05.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:56] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m06.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:56] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m07.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:57] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m08.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:58] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m09.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:59] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m10.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:45:59] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m11.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:46:00] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2009m12.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:46:01] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
ssh_y2010m01.nc
33
file format: NETCDF4
Conventions: CF-1.6
title: Juan de Fuca SSH hourly values
institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
source: https://github.com/SalishSeaCast/tools/blob/master/I_ForcingFiles/OBC/SSH_Tofino.ipynb
references: REQUIRED
history: [2015-04-23 14:46:02] Created netCDF4 zlib=True dataset.
comment: SSH from Tofino: measured - tides (from t_tide)
**** Removing ssh_y2008m12.nc and ssh_y2010m01.nc

Check the NetCDF files

A quick check that the data is saved in the netcdf file as expected.

In [10]:
filename = 'ssh_y2009m12.nc'

f=  NC.Dataset(filename,'r');

for dim in f.dimensions.values():
    print dim
    
ssh=f.variables['sossheig'];
us=f.variables['vobtcrtx'];
vs=f.variables['vobtcrty'];
print us.shape
print us[1,0,:]
print vs[:,0,:]
    
mn = ssh[:].min(); print mn
mx = ssh[:].max(); print mx

fig, ((ax_net,ax_data)) = plt.subplots(1, 2, figsize=(14,4))
ax_net.plot(ssh[:,0,0])
ax_net.set_xlabel(filename)
ax_net.axis([10*24, 17*24, mn,mx])

ax_data.plot(wlev_meas.time,ssanomaly)
ax_data.set_xlim(['Dec 11 2009', 'Dec 18 2009'])
ax_data.set_ylim([mn,mx])
<type 'netCDF4.Dimension'>: name = 'xbT', size = 87

<type 'netCDF4.Dimension'>: name = 'yb', size = 1

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

(744, 1, 87)
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[[ 0.  0.  0. ...,  0.  0.  0.]
 [ 0.  0.  0. ...,  0.  0.  0.]
 [ 0.  0.  0. ...,  0.  0.  0.]
 ..., 
 [ 0.  0.  0. ...,  0.  0.  0.]
 [ 0.  0.  0. ...,  0.  0.  0.]
 [ 0.  0.  0. ...,  0.  0.  0.]]
-0.156304
0.479705
Out[10]:
(-0.156304, 0.47970501)

This looks good so let's move the files out of my working directory and into the forcing directory.

In [11]:
srcdir='.'
dstdir='/data/nsoontie/MEOPAR/NEMO-forcing/open_boundaries/west/ssh/'

import shutil

for basename in os.listdir(srcdir):
    if basename.endswith('.nc'):
        pathname = os.path.join(srcdir, basename)
        if os.path.isfile(pathname):
            shutil.copy2(pathname, dstdir)

Improvements to this code

There are a few ways to improve the functionality of this code.

  1. Include the matlab call to t_xtide within the notebook somwehere. At the moment, the call to matlab is performed outside of this script.

  2. Improve the method of handling local time and UTC. Right now, I have to download extra data and then delete the incomplete month files.

In [ ]: