Notebook

This notebook will load data, perform a tidal analyis, compare with observations, plot the results, and save the analysis in a spreadsheet. Eight Tidal Constituents: M2, K1, O1, S2, P1, N2, Q1 and K2 are considered.

In [31]:

# 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 [32]:

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

#the run we want to analyze
#runname = 'corr15'
runname = 'tides_test2'

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

print name

/ocean/imachuca/MEOPAR/SalishSea/results/tides_test2/

pathname for data - all of the tide runs are stored in this directory¶

path = '/data/nsoontie/MEOPAR/SalishSea/results/tides/'¶

path = '/ocean/sallen/allen/research/MEOPAR/myResults/'

the run we want to analyze¶

runname = 'corr15'¶

runname = 'oldtopog'

joining the two string together¶

name = path +runname +'/'

print name

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 [33]:

# grid
grid = NC.Dataset('/ocean/imachuca/MEOPAR/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 [34]:

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.

In [35]:

filename = '../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'})
print harm_other

                  site     lat      lon  O1_amp  O1_pha  P1_amp  P1_pha  \
0             Neah Bay  48.385 -124.616   30.90  231.50   15.50  244.60   
1         Port Renfrew  48.537 -124.476   28.30  234.80   14.07  250.60   
2         Port Angeles  48.129 -123.400   39.10  241.60   20.70  259.40   
3             Victoria  48.413 -123.399   37.00  247.80   19.70  264.60   
4        Port Townsend  48.112 -122.758   45.00  249.90   23.90  268.40   
5               Bangor  47.748 -122.727   46.60  251.90   26.00  273.90   
6              Seattle  47.605 -122.338   45.80  255.40   25.20  274.50   
7               Tacoma  47.267 -122.413   45.90  255.10   25.50  277.20   
8         Cherry Point  48.863 -122.758   45.60  260.00   25.60  281.40   
9        Friday Harbor  48.540 -123.010   42.30  256.40   23.60  274.90   
10       Hanbury Point  48.580 -123.172   43.60  253.60   23.40  271.40   
11              Sidney  48.658 -123.383   44.40  255.80   24.20  275.20   
12     Fulford Harbour  48.765 -123.453   43.00  257.80   23.40  277.80   
13        Patos Island  48.783 -122.967   45.50  262.10   24.50  284.60   
14          Tsawwassen  48.991 -123.137   47.20  261.80   25.90  282.60   
15      Point Atkinson  49.334 -123.250   48.30  263.20   26.80  283.10   
16  Winchelsea Islands  49.300 -124.083   47.70  263.50   27.40  286.20   
17        Little River  49.744 -124.918   49.26  263.94   28.62  285.67   
18         Twin Islets  50.029 -124.936   49.29  264.24   28.62  286.97   
19      Campbell River  50.042 -125.247   48.46  263.74   24.60  280.57   
20     Seymour Narrows  50.135 -125.347   41.27  254.54   21.28  271.47   
21            Owen Bay  50.311 -125.223   38.19  251.34   20.97  267.47   
22             Big Bay  50.394 -125.136   46.63  262.44   25.33  282.07   
23       Chatham Point  50.332 -125.441   37.46  249.04   20.39  265.97   
24        Yorke Island  50.444 -125.975   32.16  241.04   17.10  257.67   
25           Alert Bay  50.588 -126.937   30.60  239.84   16.00  251.77   
26          Port Hardy  50.720 -127.476   29.70  233.50   15.40  245.50   
27       Montagu Point  50.639 -126.213   31.10  237.60   16.60  251.30   
28          Siwash Bay  50.680 -125.763   31.30  239.40   17.10  253.20   
29      Winter Harbour  50.490 -128.044   27.26  231.20   13.39  242.90   
30         Bella Bella  52.177 -128.111   27.80  236.20   14.20  247.20   
31              Tofino  49.144 -125.937   24.50  227.20   12.30  237.90   

    Q1_amp  Q1_pha  S2_amp  S2_pha  N2_amp  N2_pha  K2_amp  K2_pha  
0     5.50  222.10   22.80   272.6   16.60  222.80    6.00  266.40  
1     5.04  225.90   21.04   268.7   15.15  217.30    4.92  263.10  
2     6.60  232.80   14.60   326.4   11.60  280.10    2.70  332.70  
3     6.10  236.00   10.20   332.8    9.10  292.00    2.00  341.90  
4     7.40  243.60   16.80    13.0   14.20  321.80    5.00   18.30  
5     8.00  247.20   25.70    29.5   20.80  333.50    7.30   28.50  
6     7.50  250.60   25.80    37.9   21.20  340.20    7.20   36.50  
7     7.60  250.60   28.20    37.8   22.50  341.20    8.20   39.60  
8     7.60  253.20   17.90    50.3   15.40  354.50    5.00   50.50  
9     6.80  244.00   13.30    34.9   12.20  341.30    3.50   40.60  
10    7.50  247.00   12.70    18.0   11.30  324.90    3.80   37.90  
11    7.50  247.00   13.20    26.8   12.00  334.60    3.80   37.90  
12    7.00  251.60   13.90    37.2   11.90  342.60    3.90   40.00  
13    7.80  253.20   16.70    54.8   14.30  354.20    4.90   58.50  
14    6.90  258.50   20.00    55.0   17.20    0.20    5.60   59.40  
15    7.70  258.80   22.90    59.9   18.40    2.90    6.20   59.90  
16    8.00  257.40   23.60    62.0   20.60    5.60    6.40   64.60  
17    8.38  257.20   25.02    61.6   21.64    5.42    6.80   62.56  
18    7.89  258.59   25.82    64.8   21.82    9.12    6.92   63.66  
19    8.08  252.39   20.27    43.6   19.20    2.82    5.42   49.76  
20    7.25  244.99   30.27   339.6   20.48  290.52    8.29  333.06  
21    6.37  244.89   27.52   339.6   17.89  290.92    6.89  335.26  
22    8.20  224.79   19.29    35.3   15.94  346.02    4.72   35.56  
23    5.82  243.69   29.44   326.8   19.57  276.22    8.05  322.36  
24    5.33  234.89   38.56   301.2   25.73  248.12   10.70  293.76  
25    5.18  231.09   40.63   290.0   26.97  237.72   11.19  279.96  
26    5.00  224.30   42.00   281.4   27.30  227.80   10.90  276.20  
27    5.20  230.40   49.60   292.7   31.60  238.70   12.50  285.50  
28    5.20  232.50   50.60   296.7   32.70  242.50   14.00  290.00  
29    4.89  224.50   29.55   273.1   20.74  219.00    7.87  265.80  
30    4.90  225.60   40.10   280.0   27.10  227.50   10.90  271.10  
31    4.40  219.60   27.90   269.5   20.30  215.60    7.60  261.60

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 [36]:

  stations =  ['PortRenfrew','SheringhamPoint','PedderBay', 'Esquimalt',
             'Victoria','CloverPoint','FinnertyCove', 'FulfordHarbour',
            'TumboChannel','PatosIsland','WhalerBay', 'Tsawwassen',
              'Sandheads', 'PointGrey','PointAtkinson','GibsonsLanding', 'WinchelseaIs',
             'HalfmoonBay','IrvinesLanding','PowellRiver', 'LittleRiver', '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', 'Winchelsea',
                    'Halfmoon Bay','Irvines Landing','Powell River', 'Little 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 [37]:

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 Winchelsea
17 Halfmoon Bay
18 Irvines Landing
19 Powell River
20 Little River
21 Lund
22 Twin Islets
23 Campbell River
24 Maude Island E
25 Nymphe Cove
26 Seymour Narrows
27 Brown Bay
28 Chatham Point
29 Kelsey Bay
30 Yorke Island

Out[37]:

[-126.1, -122, 47, 51]

Note: Some day it would be worthwhile to place the numbers more carefully so that they don't overlap.

Tidal Harmonics¶

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

In [38]:

#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.
#O1
O1freq = 13.943036*np.pi/180.
#S2
S2freq = 30.000002*np.pi/180.
#P1
P1freq = 14.958932*np.pi/180.
#N2
N2freq = 28.439730*np.pi/180.
#Q1
Q1freq = 13.398661*np.pi/180.
#K2
K2freq = 30.082138*np.pi/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
O1ft = 1.081364
O1uvt = 312.950020
S2ft = 1.0
S2uvt = 0.0
P1ft = 1.0
P1uvt = 55.79460
N2ft = 0.98784
N2uvt = 353.570277
Q1ft = 1.081364
Q1uvt = 60.631733
K2ft = 1.114095
K2uvt = 52.129248

# for start of Apr 21, 2003
new = 'true'
if new == 'true':

    K1ft = 1.065505
    K1uvt = 111.481741
    M2ft = 0.982328
    M2uvt = 250.506179
    O1ft = 1.105495
    O1uvt = 142.040782
    S2ft = 1.000000  
    S2uvt = 0.000000
    P1ft = 1.000000
    P1uvt = 241.335269
    N2ft = 0.982328
    N2uvt = 205.684028
    Q1ft = 1.105495 
    Q1uvt = 97.218631
    K2ft = 1.159036 
    K2uvt = 42.361669

In [39]:

# 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.))

In [40]:

# function for fitting 3 frequencies
def triple(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha):
    return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
            K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+
            O1amp*np.cos(O1freq*x-O1pha*np.pi/180.))

In [41]:

# function for fitting 4 frequencies
def quad(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha):
    return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
            K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+
            O1amp*np.cos(O1freq*x-O1pha*np.pi/180.)+
            S2amp*np.cos(S2freq*x-S2pha*np.pi/180.))

In [42]:

# function for fitting 6 frequencies
def sextuple(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha,
                 P1amp, P1pha, N2amp, N2pha):
    return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
            K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+
            O1amp*np.cos(O1freq*x-O1pha*np.pi/180.)+
            S2amp*np.cos(S2freq*x-S2pha*np.pi/180.)+
            P1amp*np.cos(P1freq*x-P1pha*np.pi/180.)+
            N2amp*np.cos(N2freq*x-N2pha*np.pi/180.))

In [43]:

# function for fitting 8 frequencies
def octuple(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha,
                 P1amp, P1pha, N2amp, N2pha, Q1amp, Q1pha, K2amp, K2pha):
    return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
            K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+
            O1amp*np.cos(O1freq*x-O1pha*np.pi/180.)+
            S2amp*np.cos(S2freq*x-S2pha*np.pi/180.)+
            P1amp*np.cos(P1freq*x-P1pha*np.pi/180.)+
            N2amp*np.cos(N2freq*x-N2pha*np.pi/180.)+
            Q1amp*np.cos(Q1freq*x-Q1pha*np.pi/180.)+
            K2amp*np.cos(K2freq*x-K2pha*np.pi/180.))

Now we can apply this fit to our model output.

In [44]:

fig, ax = plt.subplots(1,1,figsize=(12,5))
for stn in (0,4,14,23):
    print stations[stn]
    fT1 = NC.Dataset(name+stations[stn]+'.nc','r')
    time = fT1.variables["time_counter"][:]/3600.  # want hours not seconds
    ssh = fT1.variables["sossheig"][:,0,0]
    ax.plot(time,ssh)

PortRenfrew
Victoria
PointAtkinson
CampbellRiver

In [45]:

print ssh.shape
print 25*48*2, 30*24*4

(3840,)
2400 2880

In [46]:

#allocate space for our arrays
M2_amp=[]; M2_pha=[]; K1_amp=[]; K1_pha=[]
O1_amp=[]; O1_pha=[]; S2_amp=[]; S2_pha=[]
P1_amp=[]; P1_pha=[]; N2_amp=[]; N2_pha=[]
Q1_amp=[]; Q1_pha=[]; K2_amp=[]; K2_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)
O1_amp_obs=np.zeros(numsta); O1_pha_obs=np.zeros(numsta)
S2_amp_obs=np.zeros(numsta); S2_pha_obs=np.zeros(numsta)
P1_amp_obs=np.zeros(numsta); P1_pha_obs=np.zeros(numsta)
N2_amp_obs=np.zeros(numsta); N2_pha_obs=np.zeros(numsta)
Q1_amp_obs=np.zeros(numsta); Q1_pha_obs=np.zeros(numsta)
K2_amp_obs=np.zeros(numsta); K2_pha_obs=np.zeros(numsta)

ts = 240
te = ssh.shape[0]

  
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(octuple,time[ts:te],ssh[ts:te]) 
    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
    elif pha < 0:
        pha = pha+360
    if stn == 6:
        print pha
    M2_pha.append(pha)

    if fitted[2] < 0:
        fitted[2] = - fitted[2]
        fitted[3] = fitted[3] + 180
    K1_amp.append(fitted[2]*K1ft)
    pha = fitted[3] + K1uvt
    if  pha > 360:
        pha = pha-360
    K1_pha.append(pha)  
    
    if fitted[4] < 0:
        fitted[4] = -fitted[4]
        fitted[5] = fitted[5]+180
    O1_amp.append(fitted[4]*O1ft)
    pha= fitted[5]+O1uvt
    if  pha > 360:
        pha=pha-360
    O1_pha.append(pha) 
    
    if fitted[6] < 0:
        fitted[6] = -fitted[6]
        fitted[7] = fitted[7]+180
    S2_amp.append(fitted[6]*S2ft)
    pha= fitted[7]+S2uvt
    if  pha > 360:
        pha=pha-360
    S2_pha.append(pha) 
    
    if fitted[8] < 0:
            fitted[8] = -fitted[8]
            fitted[9] = fitted[9]+180
    P1_amp.append(fitted[8]*P1ft)
    pha= fitted[9]+P1uvt
    if  pha > 360:
        pha=pha-360
    P1_pha.append(pha) 
    
    if fitted[10] < 0:
            fitted[10] = -fitted[10]
            fitted[11] = fitted[11]+180
    N2_amp.append(fitted[10]*N2ft)
    pha= fitted[11]+N2uvt
    if  pha > 360:
        pha=pha-360
    N2_pha.append(pha) 
    
    if fitted[12] < 0:
            fitted[12] = -fitted[12]
            fitted[13] = fitted[13]+180
    Q1_amp.append(fitted[12]*Q1ft)
    pha= fitted[13]+Q1uvt
    if  pha > 360:
        pha=pha-360
    Q1_pha.append(pha) 
    
    if fitted[14] < 0:
            fitted[14] = -fitted[14]
            fitted[15] = fitted[15]+180
    K2_amp.append(fitted[14]*K2ft)
    pha= fitted[15]+K2uvt
    if pha > 360:
        pha = pha-360
    if pha < 0:
        pha = pha + 360
    K2_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]
    
    #O1/S2/P1/N2/Q1/K2 are in the other file
    if (harm_other.site==location).any():
        O1_amp_obs[stn]=harm_other.O1_amp[harm_other.site==location]/100
        O1_pha_obs[stn]=harm_other.O1_pha[harm_other.site==location]
        S2_amp_obs[stn]=harm_other.S2_amp[harm_other.site==location]/100
        S2_pha_obs[stn]=harm_other.S2_pha[harm_other.site==location]
        P1_amp_obs[stn]=harm_other.P1_amp[harm_other.site==location]/100
        P1_pha_obs[stn]=harm_other.P1_pha[harm_other.site==location]
        N2_amp_obs[stn]=harm_other.N2_amp[harm_other.site==location]/100
        N2_pha_obs[stn]=harm_other.N2_pha[harm_other.site==location]
        Q1_amp_obs[stn]=harm_other.Q1_amp[harm_other.site==location]/100
        Q1_pha_obs[stn]=harm_other.Q1_pha[harm_other.site==location]
        K2_amp_obs[stn]=harm_other.K2_amp[harm_other.site==location]/100
        K2_pha_obs[stn]=harm_other.K2_pha[harm_other.site==location]
    #Mask the arrays so that we can do statistics without the 0's throwing thigns off.
    O1_amp_obs =np.ma.masked_values(O1_amp_obs, 0)
    O1_pha_obs =np.ma.masked_values(O1_pha_obs, 0)
    S2_amp_obs =np.ma.masked_values(S2_amp_obs, 0)
    S2_pha_obs =np.ma.masked_values(S2_pha_obs, 0)
    P1_amp_obs =np.ma.masked_values(P1_amp_obs, 0)
    P1_pha_obs =np.ma.masked_values(P1_pha_obs, 0)
    N2_amp_obs =np.ma.masked_values(N2_amp_obs, 0)
    N2_pha_obs =np.ma.masked_values(N2_pha_obs, 0)
    Q1_amp_obs =np.ma.masked_values(Q1_amp_obs, 0)
    Q1_pha_obs =np.ma.masked_values(Q1_pha_obs, 0)
    K2_amp_obs =np.ma.masked_values(K2_amp_obs, 0)
    K2_pha_obs =np.ma.masked_values(K2_pha_obs, 0)

359.453826095

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 [47]:

#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[47]:

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

In [48]:

#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[48]:

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

In [49]:

#Plotting - O1

fig=tidetools.plot_scatter_pha_amp(O1_amp,O1_amp_obs,O1_pha,O1_pha_obs,'O1',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[49]:

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

In [50]:

#Plotting - S2

fig=tidetools.plot_scatter_pha_amp(S2_amp,S2_amp_obs,S2_pha,S2_pha_obs,'S2',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[50]:

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

In [51]:

#Plotting - P1

fig=tidetools.plot_scatter_pha_amp(P1_amp,P1_amp_obs,P1_pha,P1_pha_obs,'P1',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[51]:

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

In [52]:

#Plotting - N2

fig=tidetools.plot_scatter_pha_amp(N2_amp,N2_amp_obs,N2_pha,N2_pha_obs,'N2',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[52]:

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

In [53]:

#Plotting - Q1

fig=tidetools.plot_scatter_pha_amp(Q1_amp,Q1_amp_obs,Q1_pha,Q1_pha_obs,'Q1',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[53]:

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

In [54]:

#Plotting - K2

fig=tidetools.plot_scatter_pha_amp(K2_amp,K2_amp_obs,K2_pha,K2_pha_obs,'K2',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[54]:

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

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 [55]:

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

RMS Error¶

In [56]:

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 [57]:

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 [58]:

#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)*1.03,np.array(M2_pha)+2.3)
D_K1= complex_diff(np.array(K1_amp_obs),np.array(K1_pha_obs), np.array(K1_amp)*0.99,np.array(K1_pha)-0.5)
D_O1= complex_diff(np.ma.array(O1_amp_obs),np.ma.array(O1_pha_obs), np.ma.array(O1_amp),np.ma.array(O1_pha))
D_S2= complex_diff(np.ma.array(S2_amp_obs),np.ma.array(S2_pha_obs), np.ma.array(S2_amp),np.ma.array(S2_pha))
D_P1= complex_diff(np.ma.array(P1_amp_obs),np.ma.array(P1_pha_obs), np.ma.array(P1_amp),np.ma.array(P1_pha))
D_N2= complex_diff(np.ma.array(N2_amp_obs),np.ma.array(N2_pha_obs), np.ma.array(N2_amp),np.ma.array(N2_pha))
D_Q1= complex_diff(np.ma.array(Q1_amp_obs),np.ma.array(Q1_pha_obs), np.ma.array(Q1_amp),np.ma.array(Q1_pha))
D_K2= complex_diff(np.ma.array(K2_amp_obs),np.ma.array(K2_pha_obs), np.ma.array(K2_amp),np.ma.array(K2_pha))

Saving the results¶

We will now save these statistics in a spreadsheet

In [59]:

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 no PR',
                    mean(M2_amp[1:split2]-M2_amp_obs[1:split2]),mean(Dphi_M2[1:split2]),mean(D_M2[1:split2]), 
                    mean(K1_amp[1:split2]-K1_amp_obs[1:split2]),mean(Dphi_K1[1:split2]),mean(D_K1[1:split2])])
    writer.writerow(['RMS Difference no North no PR',
                    rms(M2_amp[1:split2]-M2_amp_obs[1:split2]),rms(Dphi_M2[1:split2]),rms(D_M2[1:split2]), 
                    rms(K1_amp[1:split2]-K1_amp_obs[1:split2]),rms(Dphi_K1[1:split2]),rms(D_K1[1: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 [31]:

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

plt.subplot(3,2,1)
plt.plot(np.array(M2_amp)*1.03, '-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(np.array(K1_amp)*0.99, '-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
pha_uwm = 180./np.pi * np.unwrap(np.array(M2_pha)*np.pi/180.)
plt.plot(pha_uwm+2.3, '-bo', label = 'model')
pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
plt.plot(pha_uw, 'r-o', label = 'observation')
plt.title('M2 Phase')

plt.subplot(3,2,4)
pha_uw = 180./np.pi * np.unwrap(np.array(K1_pha)*np.pi/180.)
plt.plot(pha_uw-0.5, '-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[31]:

<matplotlib.legend.Legend at 0x7f598b157a10>

In [32]:

M2_amp_topogbottfric = M2_amp
M2_pha_topogbottfric = pha_uwm

In [33]:

M2_amp_rnshlat2 = M2_amp
M2_pha_rnshlat2 = pha_uwm

In [34]:

M2_amp_bot1em2B = M2_amp
M2_pha_bot1em2B = pha_uwm

In [35]:

M2_amp_bot1em2 = M2_amp
M2_pha_bot1em2 = pha_uwm

In [36]:

M2_amp_corr15 = M2_amp
M2_pha_corr15 = pha_uwm

In [37]:

M2_amp_topog = M2_amp
M2_pha_topog = pha_uwm

In [38]:

M2_amp_rnshlat = M2_amp
M2_pha_rnshlat = pha_uwm

In [39]:

plt.figure(figsize=(12,5))
plt.plot(np.array(M2_amp_topog)*0.97, '-bo', label = 'topog')
plt.plot(M2_amp_obs, 'r-s', label = 'observation')
plt.plot(np.array(M2_amp_rnshlat2)*1.08, '-m^')
plt.plot(np.array(M2_amp_bot1em2B)*1.07, '-c^')
plt.plot(M2_amp_corr15, '-go', label='corr15')
plt.plot(np.array(M2_amp_topogbottfric)*1.03, '-yo')
makeit = np.array(M2_amp_corr15) + (np.array(M2_amp_bot1em2B)*1.07-np.array(M2_amp_corr15)) + (np.array(M2_amp_topog)*0.97-np.array(M2_amp_corr15))
plt.plot(makeit, 'k*-')

Out[39]:

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

In [40]:

plt.figure(figsize=(12,5))
plt.plot(np.array(M2_pha_topog)+2.9, '-bo', label = 'topog')
pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
plt.plot(pha_uw, 'r-s', label = 'observation')
plt.plot(np.array(M2_pha_rnshlat2)-0.8, '-m^')
plt.plot(np.array(M2_pha_bot1em2B)-0.5, '-c^')
plt.plot(M2_pha_corr15, '-go', label = 'corr15')
plt.plot(np.array(M2_pha_topogbottfric)+2, '-yo')
makeit = np.array(M2_pha_corr15) + (np.array(M2_pha_bot1em2B)-0.8-np.array(M2_pha_corr15)) + (np.array(M2_pha_topog)+2.9-np.array(M2_pha_corr15))
plt.plot(makeit, 'k*-')

Out[40]:

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

In [41]:

cmap = plt.get_cmap('PuBu')
cmap.set_bad('burlywood')

fig,axs=plt.subplots(4, 2, figsize=(8,20))

constituent = ('M2', 'K1', 'O1', 'S2', 'P1', 'N2', 'Q1', 'K2')
error_D = (D_M2, D_K1, D_O1, D_S2, D_P1, D_N2, D_Q1, D_K2)


for row in range(4):

    for ax, error_D1, const in zip(axs[row], error_D[row*2:row*2+2], constituent[row*2:row*2+2]):
        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.set_title(const)
        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 [42]:

fig, axs = plt.subplots(6,2,figsize=(10,15))
axs[0,0].plot(np.array(O1_amp)/np.array(K1_amp), '-bo', label = 'model')
axs[0,0].plot((0,28),(0.560,0.560), 'r-', label = 'observation')
axs[0,0].set_title('O1/K1 Amplitude')
axs[0,1].plot(np.array(O1_pha)-np.array(K1_pha), '-bo', label = 'model')
axs[0,1].plot((0,28),(-22.9,-22.9), 'r-', label = 'observation')
axs[0,1].set_title('O1-K1 Phase')

axs[1,0].plot(np.array(S2_amp)/np.array(M2_amp), '-bo', label = 'model')
axs[1,0].plot((0,28),(0.249,0.249), 'r-', label = 'observation')
axs[1,0].set_title('S2/M2 Amplitude')
pha_uw = 180./np.pi * np.unwrap((np.array(S2_pha)-np.array(M2_pha))*np.pi/180.)
axs[1,1].plot(pha_uw, '-bo', label = 'model')
axs[1,1].plot((0,28),( 28.7, 28.7), 'r-', label = 'observation')
axs[1,1].set_title('S2-M2 Phase')

axs[2,0].plot(np.array(P1_amp)/np.array(K1_amp), '-bo', label = 'model')
axs[2,0].plot((0,28),(0.311,0.311), 'r-', label = 'observation')
axs[2,0].set_title('P1/K1 Amplitude')

pha_uw = 180./np.pi * np.unwrap((np.array(P1_pha)-np.array(K1_pha))*np.pi/180.)
axs[2,1].plot(pha_uw, '-bo', label = 'model')
axs[2,1].plot((0,28),(-3,-3), 'r-', label = 'observation')
axs[2,1].set_title('P1-K1 Phase')

axs[3,0].plot(np.array(N2_amp)/np.array(M2_amp), '-bo', label = 'model')
axs[3,0].plot((0,28),(0.200,0.200), 'r-', label = 'observation')
axs[3,0].set_title('N2/M2 Amplitude')

pha_uw = 180./np.pi * np.unwrap((np.array(N2_pha)-np.array(M2_pha))*np.pi/180.)
axs[3,1].plot(pha_uw, '-bo', label = 'model')
axs[3,1].plot((0,28),(-28.3, -28.3), 'r-', label = 'observation')
axs[3,1].set_title('N2-M2 Phase')

axs[4,0].plot(np.array(Q1_amp)/np.array(K1_amp), '-bo', label = 'model')
axs[4,0].plot((0,28),(0.089,0.089), 'r-', label = 'observation')
axs[4,0].set_title('Q1/K1 Amplitude')

pha_uw = 180./np.pi * np.unwrap((np.array(Q1_pha)-np.array(K1_pha))*np.pi/180.)
axs[4,1].plot(pha_uw+360., '-bo', label = 'model')
axs[4,1].plot((0,28),(-27.3,-27.3), 'r-', label = 'observation')
axs[4,1].set_title('Q1-K1 Phase')

axs[5,0].plot(np.array(K2_amp)/np.array(M2_amp), '-bo', label = 'model')
axs[5,0].plot((0,28),(0.068,0.068), 'r-', label = 'observation')
axs[5,0].set_title('K2/M2 Amplitude')

pha_uw = 180./np.pi * np.unwrap((np.array(K2_pha)-np.array(M2_pha))*np.pi/180.)
axs[5,1].plot(pha_uw, '-bo', label = 'model')
axs[5,1].plot((0,28),(28.7, 28.7), 'r-', label = 'observation')
axs[5,1].set_title('K2-M2 Phase')

Out[42]:

<matplotlib.text.Text at 0x7f5935c30850>

In [43]:

print te

In [44]:

sample = 17
start = np.zeros(sample)
tend = np.zeros(sample)
for i in range(sample):
    start[i] = 196+(480-196)*np.random.rand()
    tend[i] = te-(480-196)*np.random.rand()
print start
print tend
timelength = (tend-start)/96.
print np.mean(timelength),2*np.std(timelength)
print time[start[1]:tend[1]]

[ 477.58594177  251.70925818  302.80288928  462.29204484  474.77495378
  424.60064219  293.30229037  207.04427899  387.72073575  240.58269087
  413.30969544  226.68384062  341.59429178  198.60234818  407.94633113
  234.94516407  406.83508316]
[ 3638.45825492  3828.93243835  3693.41964051  3699.52463167  3735.29513176
  3637.43650941  3591.68710002  3668.89524392  3713.6094466   3573.6362958
  3707.77812636  3798.52219798  3753.69133443  3795.20079988  3590.21827328
  3793.96900934  3703.26582262]
35.031377314 2.9175750749
[  62.875   63.125   63.375 ...,  956.375  956.625  956.875]

In [45]:

#allocate space for our arrays
M2_amp=np.zeros((numsta,sample)); M2_pha=np.zeros((numsta,sample))
K1_amp=np.zeros((numsta,sample)); K1_pha=np.zeros((numsta,sample))
O1_amp=np.zeros((numsta,sample)); O1_pha=np.zeros((numsta,sample))
S2_amp=np.zeros((numsta,sample)); S2_pha=np.zeros((numsta,sample))
P1_amp=np.zeros((numsta,sample)); P1_pha=np.zeros((numsta,sample))
N2_amp=np.zeros((numsta,sample)); N2_pha=np.zeros((numsta,sample))
Q1_amp=np.zeros((numsta,sample)); Q1_pha=np.zeros((numsta,sample))
K2_amp=np.zeros((numsta,sample)); K2_pha=np.zeros((numsta,sample))



for it,tst,tet in zip(range(sample),start,tend):
  
    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(octuple,time[tst:tet],ssh[tst:tet]) 
        if fitted[0] < 0:
            fitted[0] = -fitted[0]
            fitted[1] = fitted[1]+180

        M2_amp[stn,it] = fitted[0]*M2ft
        pha = fitted[1]+M2uvt
        if  pha > 360:
            pha=pha-360
        M2_pha[stn,it] = pha
        
        if fitted[2] < 0:
            fitted[2] = -fitted[2]
            fitted[3] = fitted[3]+180

        K1_amp[stn,it] = fitted[2]*K1ft
        pha= fitted[3]+K1uvt
        if  pha > 360:
            pha=pha-360
        K1_pha[stn,it]= pha   
                
        if fitted[4] < 0:
            fitted[4] = -fitted[4]
            fitted[5] = fitted[5]+180
        O1_amp[stn,it] =fitted[4]*O1ft
        pha= fitted[5]+O1uvt
        if  pha > 360:
            pha=pha-360
        O1_pha[stn,it]= pha 
        
        if fitted[6] < 0:
            fitted[6] = -fitted[6]
            fitted[7] = fitted[7]+180
        S2_amp[stn,it] =fitted[6]*S2ft
        pha= fitted[7]+S2uvt
        if  pha > 360:
            pha=pha-360
        S2_pha[stn,it]= pha 
        
        if fitted[8] < 0:
            fitted[8] = -fitted[8]
            fitted[9] = fitted[9]+180
        P1_amp[stn,it] = fitted[8]*P1ft
        pha= fitted[9]+P1uvt
        if  pha > 360:
            pha=pha-360
        P1_pha[stn,it] =pha 
    
        if fitted[10] < 0:
            fitted[10] = -fitted[10]
            fitted[11] = fitted[11]+180
        N2_amp[stn,it] = fitted[10]*N2ft
        pha= fitted[11]+N2uvt
        if  pha > 360:
            pha=pha-360
        N2_pha[stn,it] = pha
        
        if fitted[12] < 0:
            fitted[12] = -fitted[12]
            fitted[13] = fitted[13]+180
        Q1_amp[stn,it] = fitted[12]*Q1ft
        pha= fitted[13]+Q1uvt
        if  pha > 360:
            pha=pha-360
        Q1_pha[stn,it] = pha
        
        if fitted[14] < 0:
            fitted[14] = -fitted[14]
            fitted[15] = fitted[15]+180
        K2_amp[stn,it] = fitted[14]*K2ft
        pha= fitted[15]+K2uvt
        if  pha > 360:
            pha=pha-360
        K2_pha[stn,it] = pha
    

In [46]:

jdef = range(3)
south = range(14,18)
north = range(29,31)
print 'M2'
print '     JdeFuca'
print np.mean(M2_amp[jdef]),2*np.std(np.mean(M2_amp[jdef],axis=0))
print np.mean(M2_amp_obs[jdef]), np.mean(M2_amp_obs[jdef])-np.mean(M2_amp[jdef])
print np.mean(M2_amp_obs[jdef])/np.mean(M2_amp[jdef])
print np.mean(M2_pha[jdef]),2*np.std(np.mean(M2_pha[jdef],axis=0))
print np.mean(M2_pha_obs[jdef]), np.mean(M2_pha_obs[jdef])-np.mean(M2_pha[jdef])
print '     South'
print np.mean(M2_amp[south]),2*np.std(np.mean(M2_amp[south],axis=0))
print np.mean(M2_amp_obs[south]), np.mean(M2_amp_obs[south])-np.mean(M2_amp[south])
print np.mean(M2_amp_obs[south])/np.mean(M2_amp[south])
print np.mean(M2_pha[south]),2*np.std(np.mean(M2_pha[south],axis=0))
print np.mean(M2_pha_obs[south]), np.mean(M2_pha_obs[south])-np.mean(M2_pha[south])
print '     North'
print np.mean(M2_amp[north]),2*np.std(np.mean(M2_amp[north],axis=0))
print np.mean(M2_amp_obs[north]), np.mean(M2_amp_obs[north])-np.mean(M2_amp[north])
print np.mean(M2_amp_obs[north])/np.mean(M2_amp[north])
print np.mean(M2_pha[north]),2*np.std(np.mean(M2_pha[north],axis=0))
print np.mean(M2_pha_obs[north]), np.mean(M2_pha_obs[north])-np.mean(M2_pha[north])
print '==============================================='
print 'K1'
print '     JdeFuca'
print np.mean(K1_amp[jdef]),2*np.std(np.mean(K1_amp[jdef],axis=0))
print np.mean(K1_amp_obs[jdef]), np.mean(K1_amp_obs[jdef])-np.mean(K1_amp[jdef])
print np.mean(K1_amp_obs[jdef])/np.mean(K1_amp[jdef])
print np.mean(K1_pha[jdef]),2*np.std(np.mean(K1_pha[jdef],axis=0))
print np.mean(K1_pha_obs[jdef]), np.mean(K1_pha_obs[jdef])-np.mean(K1_pha[jdef])
print '     South'
print np.mean(K1_amp[south]),2*np.std(np.mean(K1_amp[south],axis=0))
print np.mean(K1_amp_obs[south]), np.mean(K1_amp_obs[south])-np.mean(K1_amp[south])
print np.mean(K1_amp_obs[south])/np.mean(K1_amp[south])
print np.mean(K1_pha[south]),2*np.std(np.mean(K1_pha[south],axis=0))
print np.mean(K1_pha_obs[south]), np.mean(K1_pha_obs[south])-np.mean(K1_pha[south])
print '     North'
print np.mean(K1_amp[north]),2*np.std(np.mean(K1_amp[north],axis=0))
print np.mean(K1_amp_obs[north]), np.mean(K1_amp_obs[north])-np.mean(K1_amp[north])
print np.mean(K1_amp_obs[north])/np.mean(K1_amp[north])
print np.mean(K1_pha[north]),2*np.std(np.mean(K1_pha[north],axis=0))
print np.mean(K1_pha_obs[north]), np.mean(K1_pha_obs[north])-np.mean(K1_pha[north])
print '==============================================='

M2
     JdeFuca
0.57027031692 0.00103796378824
0.512333333333 -0.0579369835863
0.898404349889
253.096350381 0.141615031206
270.8 17.7036496191
     South
0.934552595668 0.00146821513786
0.94525 0.0106974043319
1.01144655141
31.3392568628 0.0316530577237
31.35 0.0107431372001
     North
1.15693865159 0.00097282470679
1.1705 0.0135613484117
1.0117217524
274.022116152 0.094298801686
274.05 0.0278838481535
===============================================
K1
     JdeFuca
0.561181482534 0.00482383285143
0.542666666667 -0.0185148158668
0.967007436198
257.223937011 0.60637617942
261.533333333 4.30939632272
     South
0.898845325552 0.0091415382251
0.87225 -0.0265953255516
0.970411677298
286.237346809 0.52276691814
285.95 -0.287346808551
     North
0.581489614774 0.00307231793669
0.5675 -0.0139896147741
0.975941763329
260.911889478 0.431113106731
260.7 -0.211889478139
===============================================

In [47]:

print 'O1'
print '     South'
print np.mean(O1_amp[south]/K1_amp[south]),2*np.std(np.mean(O1_amp[south]/K1_amp[south],axis=0))
print np.mean(O1_amp_obs[south]/K1_amp_obs[south]), (np.mean(O1_amp_obs[south]/K1_amp_obs[south])
                                                    -np.mean(O1_amp[south]/K1_amp[south]))
print np.mean(O1_amp_obs[south]/K1_amp_obs[south])/np.mean(O1_amp[south]/K1_amp[south])
print np.mean(O1_pha[south]-K1_pha[south]),2*np.std(np.mean(O1_pha[south]-K1_pha[south],axis=0))
print np.mean(O1_pha_obs[south]-K1_pha_obs[south]), (np.mean(O1_pha_obs[south]-K1_pha_obs[south])
                                                     -np.mean(O1_pha[south]-K1_pha[south]))
print '     North'
print np.mean(O1_amp[north]/K1_amp[north]),2*np.std(np.mean(O1_amp[north]/K1_amp[north],axis=0))
print np.mean(O1_amp_obs[north]/K1_amp_obs[north]), (np.mean(O1_amp_obs[north]/K1_amp_obs[north])
                                                    -np.mean(O1_amp[north]/K1_amp[north]))
print np.mean(O1_amp_obs[north]/K1_amp_obs[north])/np.mean(O1_amp[north]/K1_amp[north])
print np.mean(O1_pha[north]-K1_pha[north]),2*np.std(np.mean(O1_pha[north]-K1_pha[north],axis=0))
print np.mean(O1_pha_obs[north]-K1_pha_obs[north]), (np.mean(O1_pha_obs[north]-K1_pha_obs[north])
                                                     -np.mean(O1_pha[north]-K1_pha[north]))
print '==============================================='
print 'S2'
code = ('south','north')
for dir,dire in zip(code,(south,north)):
    print dir
    print np.mean(S2_amp[dire]/M2_amp[dire]),2*np.std(np.mean(S2_amp[dire]/M2_amp[dire],axis=0))
    print np.mean(S2_amp_obs[dire]/M2_amp_obs[dire]), (np.mean(S2_amp_obs[dire]/M2_amp_obs[dire])
                                                    -np.mean(S2_amp[dire]/M2_amp[dire]))
    print np.mean(S2_amp_obs[dire]/M2_amp_obs[dire])/np.mean(S2_amp[dire]/M2_amp[dire])
    unwrap = np.unwrap(np.array(S2_pha)*np.pi/180.)*180./np.pi
    M2_un = np.unwrap(np.array(M2_pha)*np.pi/180.)*180./np.pi
    plt.plot (unwrap[dire],'r',M2_un[dire],'b')
    print np.mean(unwrap[dire]-M2_un[dire])+360.,2*np.std(np.mean(unwrap[dire]-M2_un[dire],axis=0))
    print np.mean(S2_pha_obs[dire]-M2_pha_obs[dire]), (np.mean(S2_pha_obs[dire]-M2_pha_obs[dire])
                                                     -np.mean(unwrap[dire]-M2_un[dire]))-360.

O1
     South
0.568224038732 0.00405451773849
0.560324825986 -0.00789921274584
0.986098418568
-22.6197245768 0.434706989358
-22.9 -0.2802754232
     North
0.579514124633 0.00240014358402
0.576344086022 -0.00317003861139
0.994529833741
-378.754622394 0.319444808762
-18.96 359.794622394
===============================================
S2
south
0.249158655972 0.00534431042522
0.249455337691 0.000296681718894
1.00119073414
388.5324582 1.59070743891
28.7 -359.8324582
north
0.329138505829 0.00607781520503
0.329291204099 0.000152698270256
1.00046393317
29.3054642426 1.24059875875
29.4 0.0945357573903

In [48]:

const = ('P1', 'Q1')
model_amp = (P1_amp, Q1_amp)
model_pha = ()
for const, model_amp, model_pha, obs_amp, obs_pha in zip(('P1','Q1'),
                                                    (P1_amp, Q1_amp),(P1_pha, Q1_pha), 
                                                  (P1_amp_obs, Q1_amp_obs), (P1_pha_obs, Q1_pha_obs)):
    print const
    for dir,dire in zip(code,(south,north)):
        print dir
        print np.mean(model_amp[dire]/K1_amp[dire]),2*np.std(np.mean(model_amp[dire]/K1_amp[dire],axis=0))
        print np.mean(obs_amp[dire]/K1_amp_obs[dire]), (np.mean(obs_amp[dire]/K1_amp_obs[dire])
                                                    -np.mean(model_amp[dire]/K1_amp[dire]))
        print np.mean(obs_amp[dire]/K1_amp_obs[dire])/np.mean(model_amp[dire]/K1_amp[dire])
        unwrap = np.unwrap(np.array(model_pha)*np.pi/180.)*180./np.pi
        K1_un = np.unwrap(np.array(K1_pha)*np.pi/180.)*180./np.pi
        print np.mean(unwrap[dire]-K1_un[dire]),2*np.std(np.mean(unwrap[dire]-K1_un[dire],axis=0))
        print np.mean(obs_pha[dire]-K1_pha_obs[dire]), (np.mean(obs_pha[dire]-K1_pha_obs[dire])
                                                     -np.mean(unwrap[dire]-K1_un[dire]))

P1
south
0.292964133189 0.0103204373458
0.31090487239 0.0179407392005
1.06123868818
-3.01904158282 1.46755505989
-3.0 0.019041582819
north
0.29503443161 0.00594792098505
0.306451612903 0.0114171812931
1.03869779277
-2.45030431581 1.21560920822
-2.33 0.120304315812
Q1
south
0.0897486251037 0.00184512033631
0.0893271461717 -0.000421478932003
0.995303783969
-26.6620584731 1.35683936235
-27.3 -0.637941526856
north
0.0977096478675 0.00128535506812
0.0955197132616 -0.00218993460582
0.977587324756
-25.0880781336 0.563972614811
-25.11 -0.0219218664093

In [49]:

const = ('N2', 'K2')
model_amp = (N2_amp, K2_amp)
model_pha = ()
for const, model_amp, model_pha, obs_amp, obs_pha in zip(('N2','K2'),
                                                    (N2_amp, K2_amp),(N2_pha, K2_pha), 
                                                  (N2_amp_obs, K2_amp_obs), (N2_pha_obs, K2_pha_obs)):
    print const
    for dir,dire in zip(code,(south,north)):
        print dir
        print np.mean(model_amp[dire]/M2_amp[dire]),2*np.std(np.mean(model_amp[dire]/M2_amp[dire],axis=0))
        print np.mean(obs_amp[dire]/M2_amp_obs[dire]), (np.mean(obs_amp[dire]/M2_amp_obs[dire])
                                                    -np.mean(model_amp[dire]/M2_amp[dire]))
        print np.mean(obs_amp[dire]/M2_amp_obs[dire])/np.mean(model_amp[dire]/M2_amp[dire])
        unwrap = model_pha#np.unwrap(np.array(model_pha)*np.pi/180.)*180./np.pi
        M2_un = np.unwrap(np.array(M2_pha)*np.pi/180.)*180./np.pi
        print unwrap[dire]
        print np.mean(unwrap[dire]-M2_un[dire]),2*np.std(np.mean(unwrap[dire]-M2_un[dire],axis=0))
        print np.mean(obs_pha[dire]-M2_pha_obs[dire]), (np.mean(obs_pha[dire]-M2_pha_obs[dire])
                                                     -np.mean(unwrap[dire]-M2_un[dire]))

N2
south
0.200668908813 0.000766960875523
0.200435729847 -0.000233178965392
0.998837991561
[[ 2.47050896  1.67152851  2.05902468  2.27650477  2.32695203  2.42408517
   2.16698289  1.74157004  2.07474311  1.8259618   2.24565942  1.99754907
   1.99218797  1.85300753  2.37859361  1.89758021  2.21695007]
 [ 2.64015103  1.86008     2.21571661  2.4356851   2.51146337  2.56839588
   2.31199055  1.94584332  2.2436573   2.00200611  2.39776706  2.16075187
   2.15872126  2.03969365  2.51093725  2.06935958  2.36733982]
 [ 4.33989174  3.53837205  3.87874201  4.15728353  4.21689654  4.25632029
   3.99149269  3.6227554   3.98834067  3.72229019  4.12787525  3.69542673
   3.85161418  3.61532551  4.24217768  3.64432485  4.09659889]
 [ 3.25225908  2.45004268  2.80812583  3.04167735  3.10261048  3.17319783
   2.93317543  2.57866554  2.85955964  2.66374531  3.00769282  2.6840446
   2.74660455  2.60547405  3.14991152  2.61235508  2.97977409]]
-28.549645192 0.474413458632
-28.3 0.249645191984
north
0.220278013553 0.0012382523029
0.219726729291 -0.000551284261964
0.997497325071
[[ 250.67598328  250.93328407  250.73248814  250.82878225  250.85223144
   250.65400404  250.61265863  250.71740972  250.76947646  250.61424709
   250.74224802  250.90318569  250.81750521  250.9095643   250.55672696
   250.90904825  250.74081281]
 [ 249.68111681  249.87005543  249.69041736  249.8305825   249.85099172
   249.6391468   249.56960222  249.64810204  249.7655167   249.54802581
   249.73626248  249.87138211  249.80241083  249.85536413  249.53337337
   249.86772225  249.73447295]]
-23.7731690374 0.291602526719
-23.68 0.0931690373623
K2
south
0.067660414654 0.00898869463758
0.0675381263617 -0.000122288292351
0.998192616865
[[ 57.78786901  60.48710246  59.60997698  57.71171322  56.13755834
   55.54621316  59.25861696  61.8947436   57.4739308   62.64026547
   56.32489285  60.69038447  60.06472957  62.48201792  55.47737134
   60.79414971  57.68793348]
 [ 59.40132468  61.42368389  60.8389906   58.90073687  57.43963515
   57.05801225  60.93757341  63.51230107  58.8431997   64.42639591
   57.67089456  61.67232391  61.09991259  63.40301951  57.41822419
   61.75889819  58.94173688]
 [ 57.99751481  60.47397688  59.15971992  56.45677742  55.76637384
   56.26077236  60.5633273   62.99385245  56.53603397  64.71348375
   55.40856288  60.83172836  58.39979539  61.84373159  56.36835607
   60.9265417   56.20682641]
 [ 63.22971193  65.49928062  64.85953206  62.59640549  61.29989632
   61.28080393  65.260183    67.90230816  62.89565462  68.97182562
   61.51820557  66.03394692  64.87738272  67.1095029   61.32103937
   65.99208807  62.54211017]]
29.2624135085 4.78212562543
28.7 -0.562413508518
north
0.100832445107 0.00809837349216
0.0913748932536 -0.00945755185318
0.906205271099
[[ 291.92498078  298.38337579  299.18516061  295.28105093  294.64901558
   294.25132921  296.45040936  297.21338306  297.0568058   295.54844398
   297.13762313  298.8761011   299.00107577  299.34438136  295.00841367
   299.2246711   297.1857711 ]
 [ 289.81010528  297.25117552  297.96379497  293.7573682   292.74760028
   292.22534261  294.95723009  296.39155301  295.65707967  294.36051596
   295.75684539  297.97366289  297.83904231  298.44407558  293.18569801
   298.16241555  295.88084703]]
22.0980704567 4.55170967795
21.96 -0.138070456683

In [49]:

In [50]:

fig,axs = plt.subplots(8,2,figsize=(15,25))
for i in range(sample):
    pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,i])*np.pi/180.)
    axs[0,1].plot(pha_uw ,'-ob', label = 'model')
pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
axs[0,1].plot(pha_uw, 'r-*', label = 'observation')
axs[0,1].set_title('M2 Phase')
for i in range(sample):
    axs[0,0].plot(M2_amp[:,i], '-bo', label = 'model')
axs[0,0].plot(M2_amp_obs, 'r-*', label = 'observation')
axs[0,0].set_title('M2 Amp')

for i in range(sample):
    pha_uw = 180./np.pi * np.unwrap(np.array(K1_pha[:,i])*np.pi/180.)
    axs[1,1].plot(pha_uw, '-bo', label = 'model')
axs[1,1].plot(K1_pha_obs, 'r-*', label = 'observation')
axs[1,1].set_title('K1 Phase')
for i in range(sample):
    axs[1,0].plot(K1_amp[:,i], '-bo', label = 'model')
axs[1,0].plot(K1_amp_obs, 'r-*', label = 'observation')
axs[1,0].set_title('K1 Amp')

for i in range(sample):
    if O1_pha[0,i] < 0:
        O1_pha[0,i] = O1_pha[0,i] + 360
    pha_uw = 180./np.pi * np.unwrap(np.array(O1_pha[:,i])*np.pi/180.)
    axs[2,1].plot(pha_uw, '-bo', label = 'model')
axs[2,1].plot(O1_pha_obs, 'r-*', label = 'observation', markersize = 15)
axs[2,1].set_title('O1 Phase')
for i in range(sample):
    axs[2,0].plot(O1_amp[:,i], '-bo', label = 'model')
axs[2,0].plot(O1_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[2,0].set_title('O1 Amp')

for i in range(sample):
    if S2_pha[0,i] < 0:
        S2_pha[0,i] = S2_pha[0,i] + 360
    pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,i])*np.pi/180.)
    axs[3,1].plot(pha_uw, '-bo', label = 'model')
pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha_obs)*np.pi/180.)
vsmall = 1e-6
pha_uwm = np.ma.masked_array(pha_uw, mask=(abs(pha_uw-360)<vsmall))
axs[3,1].plot(pha_uwm, 'r-*', label = 'observation', markersize = 15)
axs[3,1].set_title('S2 Phase')
for i in range(sample):
    axs[3,0].plot(S2_amp[:,i], '-bo', label = 'model')
axs[3,0].plot(S2_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[3,0].set_title('S2 Amp')

for i in range(sample):
    axs[4,0].plot(P1_amp[:,i], '-bo', label = 'model')
axs[4,0].plot(P1_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[4,0].set_title('P1 Amp')
for i in range(sample):
    pha_uw = 180./np.pi * np.unwrap(np.array(P1_pha[:,i])*np.pi/180.)
    axs[4,1].plot(pha_uw, '-bo', label = 'model')
axs[4,1].plot(P1_pha_obs, 'r-*', label = 'observation', markersize = 15)
axs[4,1].set_title('P1 Phase')

for i in range(sample):
    axs[5,0].plot(N2_amp[:,i], '-bo', label = 'model')
axs[5,0].plot(N2_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[5,0].set_title('N2 Amp')
for i in range(sample):
    pha_uw = 180./np.pi * np.unwrap(np.array(N2_pha[:,i])*np.pi/180.)
    axs[5,1].plot(pha_uw, '-bo', label = 'model')
pha_uw = 180./np.pi * np.unwrap(np.array(N2_pha_obs)*np.pi/180.)
pha_uwm = np.ma.masked_array(pha_uw, mask=(abs(pha_uw-360)<vsmall))
axs[5,1].plot(pha_uwm, 'r-*', label = 'observation', markersize = 15)  
axs[5,1].set_title('N2 Phase')

for i in range(sample):
    axs[6,0].plot(Q1_amp[:,i], '-bo', label = 'model')
axs[6,0].plot(Q1_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[6,0].set_title('Q1 Amp')
for i in range(sample):
    pha_uw = 180./np.pi * np.unwrap(np.array(Q1_pha[:,i])*np.pi/180.)
    for j in range(numsta):
        if pha_uw[j] < 0:
            pha_uw[j] += 360
    axs[6,1].plot(pha_uw, '-bo', label = 'model')
axs[6,1].plot(Q1_pha_obs, 'r-*', label = 'observation', markersize = 15)
axs[6,1].set_title('Q1 Phase')

for i in range(sample):
    axs[7,0].plot(K2_amp[:,i], '-bo', label = 'model')
axs[7,0].plot(K2_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[7,0].set_title('K2 Amp')
for i in range(sample):
    pha_uw = 180./np.pi * np.unwrap(np.array(K2_pha[:,i])*np.pi/180.)
    for j in range(numsta):
        if pha_uw[j] < 120:
            pha_uw[j] += 360
    axs[7,1].plot(pha_uw ,'-bo', label = 'model')
for j in range(numsta):
        if K2_pha_obs[j] < 120:
            K2_pha_obs[j] += 360
axs[7,1].plot(K2_pha_obs, 'r-*', label = 'observation', markersize = 15)
axs[7,1].set_title('K2 Phase')

Out[50]:

<matplotlib.text.Text at 0x7f5934e3f9d0>

In [51]:

fig,axs = plt.subplots(1,2,figsize=(15,10))
K1mean=np.average(K1_pha, axis=1)
K1max = np.max(K1_pha,axis=1)
K1min = np.min(K1_pha,axis=1)
asymmetric_error = [ K1mean-K1min, K1max-K1mean]
axs[0].errorbar(range(31),K1mean, yerr = asymmetric_error)
axs[0].plot(K1_pha_obs, 'r-*', label = 'observation')
K1mean=np.average(K1_amp, axis=1)
K1max = np.max(K1_amp,axis=1)
K1min = np.min(K1_amp,axis=1)
asymmetric_error = [ K1mean-K1min, K1max-K1mean]
axs[1].errorbar(range(31),K1mean, yerr = asymmetric_error)
axs[1].plot(K1_amp_obs, 'r-*', label = 'observation')

Out[51]:

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

Comparing corr15 and tides_test¶

In [61]:

tides_test = pd.read_csv('/ocean/imachuca/MEOPAR/analysis/Idalia/tides_test.csv')
tides_test2 = pd.read_csv('/ocean/imachuca/MEOPAR/analysis/Idalia/tides_test2.csv')
corr15 = pd.read_csv('/ocean/imachuca/MEOPAR/analysis/Idalia/corr15.csv')
oldtopog = pd.read_csv('/ocean/imachuca/MEOPAR/analysis/Idalia/oldtopog.csv')

In [90]:

tides_test['Station Name']

Out[90]:

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                        Winchelsea
17                      Halfmoon Bay
18                   Irvines Landing
19                      Powell River
20                      Little River
21                              Lund
22                       Twin Islets
23                    Campbell River
24                    Maude Island E
25                       Nymphe Cove
26                   Seymour Narrows
27                         Brown Bay
28                     Chatham Point
29                        Kelsey Bay
30                      Yorke Island
31                   Mean Difference
32                    RMS Difference
33    Mean Difference no North no PR
34     RMS Difference no North no PR
Name: Station Name, dtype: object

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.

D is complex difference

In [81]:

vars_comp=['R (M2)','Delta phi (M2)','D (M2)','R (K1)','Delta phi (K1)','D (K1)']
cols=np.arange(len(vars_comp))
fig,axs = plt.subplots(6,1,figsize=(20,24))
for col, var in zip(cols, vars_comp):
    axs[col].plot(oldtopog[var], 'g', label='oldtopog', linewidth=2)
    axs[col].plot(corr15[var], 'Magenta', label='corr15', linewidth=1.5)
    axs[col].plot(tides_test[var], 'b', label='tides_test', linewidth=2)
    axs[col].plot(tides_test2[var], 'OrangeRed', label='tides_test2', linewidth=1.5)
    axs[col].legend()
    axs[col].set_title(var)
    axs[col].xaxis.set_ticks(np.arange(0, 35, 1))
    axs[col].grid()

In [ ]: