This notebook creates monthly forcing files for the sea surface height (hourly frequency) at Port Hardy, in the northern part of the domain. 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 NEMO is in UTC but the observations come in PST I will download data rom Dec 31, 2005 to Jan 2, 2007.

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

Observations: DFO

Tidal predictions: t_ttide

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
In [2]:
#Port Hardy mmsl over 1982-2001 (in metres)
msl = 2.896

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 [3]:
yr=2009

statID = '8408'#Port Hardystation 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)
# now there is a csv file in this directory with the sear surface height information. Let's take a look at the plots.
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('Port Hardy 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('Port Hardy Observations, mean removed: ' + start+' to ' + end)
plt.xlabel('Time (UTC)')
plt.ylabel('Water Level Anomoly (m)')
UTC
2.86029745686
Out[3]:
<matplotlib.text.Text at 0x7ff18ee55390>

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_8408_31-Dec-2006_02-Jan-2007.csv', 'Port Hardy', '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/PortHardy_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.0356977391008

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]:
starti = 32
endi = 62
lengthi = endi-starti
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='Port Hardy SSH hourly values', 
    notebook_name='SSH_PortHardy', 
    nc_filepath='/data/nsoontie/MEOPAR/NEMO-forcing/open_boundaries/north/ssh/' + filename,
    comment='SSH from Port Hardy: measured - tides (from t_tide)')

    
    #dimensions
    ssh_file.createDimension('xbT', lengthi*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):
        j = 897-ir
        nav_lat[0,ir*lengthi:(ir+1)*lengthi] = lat[j,starti:endi]
        nav_lon[0,ir*lengthi:(ir+1)*lengthi] = lon[j,starti:endi]
        nbidta[0,ir*lengthi:(ir+1)*lengthi] =  range(starti,endi)
        nbjdta[0,ir*lengthi:(ir+1)*lengthi] = j
        nbrdta[0,ir*lengthi:(ir+1)*lengthi] = ir
        
    for ib in range(0,lengthi*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 = 'sshNorth_y' +str(yr_start)+'m0'+str(month_start)+'.nc'
else:
    filename = 'sshNorth_y' +str(yr_start)+'m'+str(month_start)+'.nc'
file1=filename
if month_end<10:
    file2 = 'sshNorth_y' +str(yr_end)+'m0'+str(month_end)+'.nc'
else:
    file2 = 'sshNorth_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 = 'sshNorth_y' +str(yr)+'m0'+str(m)+'.nc'
        else:
            filename = 'sshNorth_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
sshNorth_y2008m12.nc
16
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:52] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m01.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:52] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m02.nc
671
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:53] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m03.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:53] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m04.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:54] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m05.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:54] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m06.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:55] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m07.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:55] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m08.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:56] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m09.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:56] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m10.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:56] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m11.nc
720
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:57] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2009m12.nc
744
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:57] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
sshNorth_y2010m01.nc
33
file format: NETCDF4
Conventions: CF-1.6
title: Port Hardy 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_PortHardy.ipynb
references: REQUIRED
history: [2015-04-23 16:58:58] Created netCDF4 zlib=True dataset.
comment: SSH from Port Hardy: measured - tides (from t_tide)
**** Removing sshNorth_y2008m12.nc and sshNorth_y2010m01.nc

Check the NetCDF files

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

In [9]:
filename = 'sshNorth_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[:,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,1])
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 = 30

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

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

(744, 1, 30)
[[ 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.270444
0.407423
Out[9]:
(-0.27044401, 0.40742299)

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

In [10]:
srcdir='.'
dstdir='/data/nsoontie/MEOPAR/NEMO-forcing/open_boundaries/north/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)
In [ ]: