Notebook

This notebook will load data, perform a tidal analyis, compare with observations, plot the results, and save the analysis in a spreadsheet.

In [30]:

# imports
%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import netCDF4 as NC
from scipy.optimize import curve_fit
from salishsea_tools import tidetools
from salishsea_tools import viz_tools
from salishsea_tools import bathy_tools
import collections
import pandas as pd
import csv
import math

from __future__ import division

Run Details¶

First, let's define the run that we will be analyzing. We can analyze a different run by changing runname in the cell below. A spreadsheet called tide_runs.ods contains a list of runs that we can look at.

In [31]:

# pathname for data - all of the tide runs are stored in this directory
path = '/data/nsoontie/MEOPAR/SalishSea/results/tides/'
path = '../../myresults/'

#the run we want to analyze
runname = 'CBase_FTK1b'

#joining the two string together
name = path +runname +'/'

print name

../../myresults/CBase_FTK1b/

We'll also load the bathymetry data in case we want to look at that. The package tidetools has a function get_SS_bathy_data() that returns bathymetry and grid data.

In [32]:

# grid
grid = NC.Dataset('../../nemo-forcing/grid/bathy_meter_SalishSea2.nc')
bathy, X, Y = tidetools.get_bathy_data(grid)

Observations¶

Next, we can load some observations from a text file: /data/nsoontie/MEOPAR/analysis/compare_tides/obs_tidal_wlev_const_all.csv Note: This file contains a mix of M2/K1 measurements from Foreman et al (1995), US tidal harmonics, Foreman et al (2004) and Foreman et al (2012) (for Northern tides).

In [33]:

filename = '/data/nsoontie/MEOPAR/analysis/compare_tides/obs_tidal_wlev_const_all.csv'
filename = '../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'})
print harm_obs

                  site       lat       lon  M2_amp  M2_pha  K1_amp  K1_pha
0                Sooke  48.36700  123.7330    43.8   282.7    56.9   266.4
1         Port Angeles  48.12500  123.4400    51.8   307.4    66.9   261.4
2           Pedder Bay  48.33100  123.5490    34.2   308.0    62.7   269.0
3            Esquimalt  48.43300  123.4330    36.7   317.1    64.3   268.1
4         Clover Point  48.40500  123.3470    40.3   320.3    64.2   269.8
5             Victoria  48.41700  123.3670    37.3   316.1    62.7   269.2
6        Finnerty Cove  48.47300  123.2950    44.7   357.7    70.8   277.5
7        Port Townsend  48.14500  122.7550    65.2   350.0    75.0   270.8
8               Sidney  48.65000  123.4000    55.4     5.9    76.7   277.6
9         Patricia Bay  48.65000  123.4500    60.3    14.4    76.0   281.3
10           Maple Bay  48.81700  123.6170    68.5    17.0    79.3   281.2
11     Fulford Harbour  48.76700  123.4500    58.2    12.7    75.3   280.0
12           Ladysmith  48.98300  123.8000    70.8    16.3    79.8   281.8
13        Patos Island  48.78300  122.9670    68.0    25.0    79.0   285.6
14       Tumbo Channel  48.79200  123.1080    72.6    31.0    81.1   286.9
15          Whaler Bay  48.88500  123.3250    83.4    32.9    84.7   287.5
16           Silva Bay  49.15300  123.7000    92.2    32.0    86.5   286.7
17            Ferndale  48.83300  122.7170    72.3    23.8    80.1   283.6
18              Blaine  48.99000  122.7600    77.4    25.1    82.3   284.3
19          Tsawwassen  48.99000  123.1330    81.1    27.8    83.4   284.8
20           Sandheads  49.10000  123.3000    86.9    30.9    83.7   286.5
21          Point Grey  49.25000  123.2670    94.5    33.9    90.6   287.0
22      Point Atkinson  49.33300  123.2500    91.8    31.2    86.2   286.1
23            Squamish  49.70000  123.1500    94.2    31.2    87.4   286.8
24     Gibsons Landing  49.40000  123.5000    94.7    30.1    87.2   285.2
25        Halfmoon Bay  49.51700  123.9170    96.4    31.5    88.0   285.8
26     Irvines Landing  49.63300  124.0500    98.8    31.9    88.0   286.7
27          Winchelsea  49.30000  124.0830    95.2    32.6    87.5   286.7
28       Northwest Bay  49.30000  124.2000    95.6    32.7    87.2   286.7
29        Cherry Point  48.86300  122.7570    73.2    21.8    81.5   281.9
..                 ...       ...       ...     ...     ...     ...     ...
47      Sneeoosh Point  48.40000  122.5467   102.6    18.3    78.4   282.0
48          Turner Bay  48.44500  122.5550    94.4    16.7    75.4   281.4
49     Armitage Island  48.53500  122.7967    57.3     0.5    75.6   276.4
50      Friday Harbour  48.54670  123.0100    56.5     9.7    75.8   278.8
51          Richardson  48.44670  122.9000    52.2   340.1    71.3   270.9
52        Cherry Point  48.86330  122.7567    73.4    22.8    81.7   282.8
53              Blaine  48.99167  122.7650    76.3    24.8    78.4   286.3
54        Port Renfrew  48.55000  124.4300    70.8   241.1    45.3   254.1
55        Little River  49.74000  124.9200    99.4    32.9    90.2   287.0
56         Twin Islets  50.03000  124.9300   101.3    35.4    90.4   287.5
57      Campbell River  50.04000  125.2400    82.5    18.4    84.6   284.0
58     Seymour Narrows  50.13000  125.3400    94.6   320.1    69.2   272.1
59            Owen Bay  50.31000  125.2200    85.0   319.9    67.8   272.7
60             Big Bay  50.36000  125.1300    75.5    14.9    83.3   283.5
61       Chatham Point  50.33000  125.4400    90.3   305.1    65.4   270.5
62        Yorke Island  50.44000  125.9700   117.1   271.8    55.8   260.0
63        Powell River  49.86000  124.5500   100.7    34.3    90.4   286.6
64                Lund  49.98000  124.7600   102.2    35.4    88.9   287.9
65         Nymphe Cove  50.13000  125.3600    61.5   350.4    77.0   279.9
66           Brown Bay  50.16000  125.3700    93.5   315.9    67.9   270.1
67      Maude Island E  50.13000  125.3300    55.6     7.4    81.1   283.9
68     Welsford Island  50.22000  125.1300    99.4    35.1    91.1   286.9
69         Redonda Bay  50.26000  124.9900    97.5    36.7    87.1   287.4
70     Channel Islands  50.31000  124.7500   102.6    35.9    89.9   288.0
71      Turnback Point  50.42000  125.1200   102.0    37.0    91.7   287.6
72          Orford Bay  50.59000  124.8600   101.5    37.2    90.3   288.1
73  Waddington Harbour  50.87000  124.8700   103.4    38.0    89.2   288.2
74           Shoal Bay  50.46000  125.3600    89.9   307.5    66.6   269.6
75          Kelsey Bay  50.39000  125.9600   117.0   276.3    57.7   261.4
76              Tacoma  47.26670  122.4133   113.9    11.8    83.8   277.9

[77 rows x 7 columns]

This is a list of observations that we can compare with our model output. Now we have a struc object called harm_obs that contains the data printed above.

Model¶

We don't have model output at all of the above locations. The model outputs are listed below. There is a location.nc file in the run directory for each of the stations listed below.

In [34]:

  stations =  ['PortRenfrew','SheringhamPoint','PedderBay', 'Esquimalt',
             'Victoria','CloverPoint','FinnertyCove', 'FulfordHarbour',
            'TumboChannel','PatosIsland','WhalerBay', 'Tsawwassen',
              'Sandheads', 'PointGrey','PointAtkinson','GibsonsLanding',
             'HalfmoonBay','IrvinesLanding','PowellRiver', 'Lund',
              'TwinIslets','CampbellRiver','MaudeIslandE', 'NympheCove',
              'SeymourNarrows','BrownBay','ChathamPoint','KelseyBay','YorkeIsland']
numsta=len(stations)
#again with spaces because the text file likes that
stations_obs =  ['Port Renfrew','Sheringham Point','Pedder Bay', 'Esquimalt',
                   'Victoria','Clover Point','Finnerty Cove', 'Fulford Harbour',
                    'Tumbo Channel','Patos Island','Whaler Bay', 'Tsawwassen',
                   'Sandheads', 'Point Grey','Point Atkinson','Gibsons Landing',
                    'Halfmoon Bay','Irvines Landing','Powell River', 'Lund',
                    'Twin Islets','Campbell River','Maude Island E', 'Nymphe Cove',
                    'Seymour Narrows','Brown Bay','Chatham Point','Kelsey Bay','Yorke Island']

Next, we can plot these locations on a map of our domain.

In [35]:

fig,ax=plt.subplots(1, 1, figsize=(8, 10))
ax.pcolormesh(X,Y,bathy,cmap='winter_r')

for stn in range(numsta):
    location = stations_obs[stn]
    lon=-harm_obs.lon[harm_obs.site==location]
    lat=harm_obs.lat[harm_obs.site==location]
    ax.plot(lon,lat,'.k',label=location)
    ax.annotate(stn, xy = (lon,lat), xytext = (5,5),ha = 'right', va = 'bottom',
        textcoords = 'offset points')
    print stn, location
    
ax.axis([-126.1,-122,47,51])

0 Port Renfrew
1 Sheringham Point
2 Pedder Bay
3 Esquimalt
4 Victoria
5 Clover Point
6 Finnerty Cove
7 Fulford Harbour
8 Tumbo Channel
9 Patos Island
10 Whaler Bay
11 Tsawwassen
12 Sandheads
13 Point Grey
14 Point Atkinson
15 Gibsons Landing
16 Halfmoon Bay
17 Irvines Landing
18 Powell River
19 Lund
20 Twin Islets
21 Campbell River
22 Maude Island E
23 Nymphe Cove
24 Seymour Narrows
25 Brown Bay
26 Chatham Point
27 Kelsey Bay
28 Yorke Island

Out[35]:

[-126.1, -122, 47, 51]

Note: Some day it would be worthwhile to have the numbers increase from Juan de Fuca to Johnstone Strait. And to place them more carefully so that they don't overlap.

Tidal Harmonics¶

We need a way of determing the amplitude and phase of M2/K1 from our model output. We will do this by fitting our model water levels to cosine curves with the known frequency of M2/K1.

In [36]:

#constants and fitting
# M2
M2freq = 28.984106 # degrees per hour
M2freq = M2freq*np.pi/180. # radians per hour
#K1
K1freq = 15.041069*np.pi/180.

# initialize fit
M2amp = 1.; M2pha = 180.
K1amp = 1.; K1pha = 180.

# initial phase calculation
# our start is currently Oct 26, 2002
# data for phase output from bdytides.F90; found in ocean.output
K1ft = 1.050578
K1uvt = 296.314842
M2ft = 0.987843
M2uvt = 245.888564

In [37]:

# function for fit
def double(x, M2amp, M2pha, K1amp, K1pha):
    return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
            K1amp*np.cos(K1freq*x-K1pha*np.pi/180.))

Now we can apply this fit to our model output.

In [38]:

#allocate space for our arrays
M2_amp=[]
M2_pha=[]
K1_amp=[]
K1_pha=[]

M2_amp_obs= np.zeros(numsta)
M2_pha_obs=np.zeros(numsta)
K1_amp_obs=np.zeros(numsta)
K1_pha_obs=np.zeros(numsta)

ts = 240
  
for stn in range(numsta):
    fT1 = NC.Dataset(name+stations[stn]+'.nc','r')
    time = fT1.variables["time_counter"][:]/3600.  # want hours not seconds
    ssh = fT1.variables["sossheig"][:,0,0]

    fitted, cov = curve_fit(double,time[ts:],ssh[ts:]) 
    if fitted[0] < 0:
        fitted[0] = -fitted[0]
        fitted[1] = fitted[1]+180

    M2_amp.append(fitted[0]*M2ft)
    pha = fitted[1]+M2uvt
    if  pha > 360:
        pha=pha-360
    M2_pha.append(pha)

    K1_amp.append(fitted[2]*K1ft)
    pha= fitted[3]+K1uvt
    if  pha > 360:
        pha=pha-360
    K1_pha.append(pha)   

    #now the observations
    location=stations_obs[stn]
    M2_amp_obs[stn]=harm_obs.M2_amp[harm_obs.site==location]/100
    M2_pha_obs[stn]=harm_obs.M2_pha[harm_obs.site==location]
    K1_amp_obs[stn]=harm_obs.K1_amp[harm_obs.site==location]/100
    K1_pha_obs[stn]=harm_obs.K1_pha[harm_obs.site==location]

The model data is saved in lists M2_amp, M2_pha, K1_amp, K1_pha. We have also saved the observations in M2_amp_obs, etc.

We can compare model and observations by plotting.

In [39]:

#Plotting M2
labels=['JdF/Islands','SoG','North']
split1=8; split2=20
fig=tidetools.plot_scatter_pha_amp(M2_amp,M2_amp_obs,M2_pha,M2_pha_obs,'M2',figsize=(14,6),
                                   split1=split1,split2=split2, labels=labels)

ax_amp,ax_pha = fig.axes
min_value, max_value = ax_amp.set_xlim(0, 1.2)
ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2)

min_value, max_value = ax_pha.set_xlim(0, 360)
ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2)

Out[39]:

[<matplotlib.lines.Line2D at 0x10bddfd90>]

In [40]:

#Plotting - K1

fig=tidetools.plot_scatter_pha_amp(K1_amp,K1_amp_obs,K1_pha,K1_pha_obs,'K1',figsize=(14,6),
                                   split1=split1, split2=split2, labels=labels)

ax_amp,ax_pha = fig.axes
min_value, max_value = ax_amp.set_xlim(0, 1.2)
ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2)

min_value, max_value = ax_pha.set_xlim(0, 360)
ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2)

Out[40]:

[<matplotlib.lines.Line2D at 0x107f8a290>]

The model performs well when the dots are close to the red line.

Statistics¶

We would like to save some statistics so that we can determine which runs give us the best match with observations. So, we will define some functions that will help us calculate statistics.

Mean Error (absolute value)¶

In [41]:

def mean(diff):
    return np.mean(abs(diff))

RMS Error¶

In [42]:

def rms(diff):
    return np.sqrt(np.mean(diff**2))

Complex differences¶

This is a way of measuring distances in the complex plane. We can think of our tidal amplitude and phase as a point on the complex plane. So we would like to measure the distance between a point given by the model and a point given by the observations. The function below does this.

In [43]:

def complex_diff(Ao,go,Am,gm):
    #calculates complex differences between observations and model
    #Ao, go - amplitude and phase from observations
    #Am, gm - amplitude and phase from model
    D = np.sqrt((Ao*np.cos(np.pi*go/180)-Am*np.cos(np.pi*gm/180))**2 + 
                (Ao*np.sin(np.pi*go/180)-Am*np.sin(np.pi*gm/180))**2)
    
    return D

Some other things we will look at are

$R = \frac{A_m}{A_o}$, the ratio of modelled to observed amplitude and

$\Delta \phi = \phi_m - \phi_o$, the difference betwen modelled and observed phase.

In [44]:

#R
R_M2 = M2_amp/M2_amp_obs
R_K1 = K1_amp/K1_amp_obs
#delta phi (adjust so between -180, 180)
Dphi_M2=M2_pha-M2_pha_obs; 
Dphi_M2 = Dphi_M2 -360*(Dphi_M2>180) + 360*(Dphi_M2<-180)
Dphi_K1=K1_pha-K1_pha_obs
Dphi_K1 = Dphi_K1 -360*(Dphi_K1>180) + 360*(Dphi_K1<-180)
#Complex differences
D_M2= complex_diff(np.array(M2_amp_obs),np.array(M2_pha_obs), np.array(M2_amp),np.array(M2_pha))
D_K1= complex_diff(np.array(K1_amp_obs),np.array(K1_pha_obs), np.array(K1_amp),np.array(K1_pha))

Saving the results¶

We will now save these statistics in a spreadsheet

In [45]:

outfile = runname+'.csv'

with open(outfile, 'wb') as csvfile:
    writer = csv.writer(csvfile, delimiter=',')
    writer.writerow([
            'Station Name', 
            'R (M2)', 'Delta phi (M2)', 'D (M2)',
            'R (K1)', 'Delta phi (K1)', 'D (K1)'
        ])
    for stn in range(numsta):
        location = stations_obs[stn]
        writer.writerow([stations_obs[stn],
                        R_M2[stn], Dphi_M2[stn], D_M2[stn],
                        R_K1[stn], Dphi_K1[stn], D_K1[stn]])

    #write averages and rms
    writer.writerow(['Mean Difference',
                    mean(M2_amp-M2_amp_obs),mean(Dphi_M2),mean(D_M2), 
                    mean(K1_amp-K1_amp_obs),mean(Dphi_K1),mean(D_K1)])
    writer.writerow(['RMS Difference',
                    rms(M2_amp-M2_amp_obs),rms(Dphi_M2),rms(D_M2), 
                    rms(K1_amp-K1_amp_obs),rms(Dphi_K1),rms(D_K1)])
    #without the north
    writer.writerow(['Mean Difference no North',
                    mean(M2_amp[:split2]-M2_amp_obs[:split2]),mean(Dphi_M2[:split2]),mean(D_M2[:split2]), 
                    mean(K1_amp[:split2]-K1_amp_obs[:split2]),mean(Dphi_K1[:split2]),mean(D_K1[:split2])])
    writer.writerow(['RMS Difference no North',
                    rms(M2_amp[:split2]-M2_amp_obs[:split2]),rms(Dphi_M2[:split2]),rms(D_M2[:split2]), 
                    rms(K1_amp[:split2]-K1_amp_obs[:split2]),rms(Dphi_K1[:split2]),rms(D_K1[:split2])])

Now there is a csv file in this directory with data about this run. It should be called runname.csv (where runname is the string we defined at the beginning of the notebook).

Things to try:

Add the complex differences information printed in the .csv file to tide_runs.ods. Also, if you notice any discrepancies, you can correct them. (Check M2 amplitude at Yorke Island)
Work through the notebook with a different run. There is a list of runs in tide_runs.ods
Commit and push any changes you've made to this notebook and tide_runs.odt.

Try this: * hg status (see what changes have been made) * hg in * hg commit mynotebook.ipynb (write a commit message and then save and exit) * hg commit tide_runs.odt * hg pull --rebase * pg push

Add any new csv files to the repository.

Try this: * hg add filename.csv * hg commit filename.csv * hg pull --rebase * hg push

Repeat with all the runs listed in tide_runs.odt

Plots comparing measured and model amplitudes and phases¶

In [46]:

plt.figure(figsize=(20,12))

plt.subplot(3,2,1)
plt.plot(M2_amp, '-bo', label = 'model')
plt.plot(M2_amp_obs, 'r-o', label = 'observation')
plt.title('M2 Amplitude')
plt.legend( loc='upper left' )

plt.subplot(3,2,2)
plt.plot(K1_amp, '-bo', label = 'model')
plt.plot(K1_amp_obs, 'r-o', label = 'observation')
plt.title('K1 Amplitude')

plt.subplot(3,2,3)
# use the un-wrap function to plot the M2 phase more smoothly
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha)*np.pi/180.)
plt.plot(M2_pha_uw, '-bo', label = 'model')
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
plt.plot(M2_pha_uw, 'r-o', label = 'observation')
plt.title('M2 Phase')

plt.subplot(3,2,4)
plt.plot(K1_pha, '-bo', label = 'model')
plt.plot(K1_pha_obs, 'r-o', label = 'observation')
plt.title('K1 Phase')

plt.subplot(3,2,5)
plt.plot(D_M2, '-bo', label = 'M2')
plt.plot(D_K1, '-go', label = 'K1')
plt.plot((0,30),(0.05,0.05),'k')
plt.plot((0,30),(0.10,0.10),'r')
plt.title('D error')
plt.legend( loc='upper left' )

Out[46]:

<matplotlib.legend.Legend at 0x1083b5490>

In [47]:

plt.figure(figsize=(6,6))
plt.plot(D_M2, '-co', label = 'M2')
plt.plot(D_K1, '-mo', label = 'K1')
plt.plot((0,30),(0.05,0.05),'k')
plt.plot((0,30),(0.10,0.10),'r')
plt.title('D error')
plt.legend( loc='upper left' )
cmap = plt.get_cmap('PuBu')
cmap.set_bad('burlywood')

fig,axs=plt.subplots(1, 2, figsize=(20,12))

error_D = (D_M2, D_K1)
for ax,error_D1 in zip(axs, error_D):
    ax.pcolormesh(X,Y,bathy,cmap='PuBu')

    for stn in range(numsta):
        location = stations_obs[stn]
        lon=-harm_obs.lon[harm_obs.site==location]
        lat=harm_obs.lat[harm_obs.site==location]
        if error_D1 [stn] <= 0.05:
            ax.plot(lon,lat,'og',label=location,markersize=10,markeredgecolor='g')
        if error_D1 [stn] > 0.1:
            ax.plot(lon,lat,'or',label=location,markersize=10,markeredgecolor='r')
        if 0.1 >= error_D1[stn] > 0.05:
            ax.plot(lon,lat,'oy',label=location,markersize=10,markeredgecolor='y')
        
        ax.annotate(stn, xy = (lon,lat), xytext = (5,5),ha = 'right', va = 'bottom',
            textcoords = 'offset points')
    ax.axis([-126.1,-122,47,51])
   

First Plot: M2 D error. Second Plot: K1 D error. Markers increase in size as D error increases.¶

Green: D error <= 0.05, Yellow: 0.05 < D error <= 0.1, Red: D error > 0.1

In [47]: