Tides and the sea level - Part 2
================================

Two weeks ago we calculated the sea level elevation in Stavanger related to
tidal motions given the amplitudes, periods and phases of 6 components. We are
getting back to that example again, but now we will take a look at some observed
data.

**Task 1:** First things first, we read from the dataset
:download:`stavanger.txt`. For this, we use the function ``read_csv`` from
Pandas:

.. code:: python

    import pandas as pd
    
    df = pd.read_csv('stavanger.txt', sep=r"\s+") 
    
    df.keys()

This gives us the following output:

.. code-block:: text

    Index(['time', 'ssh', 'anomaly', 'M2', 'S2', 'N2', 'K2', 'O1', 'K1', 'tide'], dtype='object')

The ``df`` Dataframe contains hourly dates between 09/2018 and 09/2019 of the
observed sea surface height (cm) in Stavanger, its anomaly (ssh - mean) and the
tidal elevations (of each component and their summation) that we calculated
previously.

**Task 2:** To proceed with the comparison between the tidal components
and the observed sea surface height, let’s select the first *45 days* of
the ``df`` dataset.

**Tip**: use the function ``iloc`` from Pandas.

.. code:: python

    import numpy as np 
    
    t = np.arange(0, 45*24+1, 1) # 45 days
    df_sm = df.iloc[0:len(t)]
    
    df_sm

We get the following output:

.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>time</th>
          <th>ssh</th>
          <th>anomaly</th>
          <th>M2</th>
          <th>S2</th>
          <th>N2</th>
          <th>K2</th>
          <th>O1</th>
          <th>K1</th>
          <th>tide</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2018-09-13</th>
          <td>00:00:00</td>
          <td>106.2</td>
          <td>40.428946</td>
          <td>14.966470</td>
          <td>6.650923</td>
          <td>1.380173</td>
          <td>1.128245</td>
          <td>1.367263</td>
          <td>1.099229</td>
          <td>26.592303</td>
        </tr>
        <tr>
          <th>2018-09-13</th>
          <td>01:00:00</td>
          <td>105.3</td>
          <td>39.528946</td>
          <td>15.634982</td>
          <td>5.426019</td>
          <td>2.258865</td>
          <td>1.173072</td>
          <td>1.174997</td>
          <td>0.745735</td>
          <td>26.413670</td>
        </tr>
        <tr>
          <th>2018-09-13</th>
          <td>02:00:00</td>
          <td>92.5</td>
          <td>26.728946</td>
          <td>12.386691</td>
          <td>2.747218</td>
          <td>2.592491</td>
          <td>0.902034</td>
          <td>0.913493</td>
          <td>0.341124</td>
          <td>19.883051</td>
        </tr>
        <tr>
          <th>2018-09-13</th>
          <td>03:00:00</td>
          <td>88.4</td>
          <td>22.628946</td>
          <td>6.035344</td>
          <td>-0.667698</td>
          <td>2.300545</td>
          <td>0.388111</td>
          <td>0.598161</td>
          <td>-0.086869</td>
          <td>8.567594</td>
        </tr>
        <tr>
          <th>2018-09-13</th>
          <td>04:00:00</td>
          <td>79.9</td>
          <td>14.128946</td>
          <td>-1.827951</td>
          <td>-3.903705</td>
          <td>1.453475</td>
          <td>-0.230316</td>
          <td>0.247583</td>
          <td>-0.508908</td>
          <td>-4.769821</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>2018-10-27</th>
          <td>20:00:00</td>
          <td>49.2</td>
          <td>-16.571054</td>
          <td>-13.860386</td>
          <td>-2.747218</td>
          <td>1.269579</td>
          <td>0.628035</td>
          <td>-0.076464</td>
          <td>1.341275</td>
          <td>-13.445178</td>
        </tr>
        <tr>
          <th>2018-10-27</th>
          <td>21:00:00</td>
          <td>50.4</td>
          <td>-15.371054</td>
          <td>-8.388605</td>
          <td>0.667698</td>
          <td>2.192932</td>
          <td>1.052723</td>
          <td>0.288132</td>
          <td>1.050467</td>
          <td>-3.136652</td>
        </tr>
        <tr>
          <th>2018-10-27</th>
          <td>22:00:00</td>
          <td>54.5</td>
          <td>-11.271054</td>
          <td>-0.815350</td>
          <td>3.903705</td>
          <td>2.587127</td>
          <td>1.193952</td>
          <td>0.635750</td>
          <td>0.687655</td>
          <td>8.192839</td>
        </tr>
        <tr>
          <th>2018-10-27</th>
          <td>23:00:00</td>
          <td>66.2</td>
          <td>0.428946</td>
          <td>6.962163</td>
          <td>6.093717</td>
          <td>2.357046</td>
          <td>1.013693</td>
          <td>0.945905</td>
          <td>0.277707</td>
          <td>17.650232</td>
        </tr>
        <tr>
          <th>2018-10-28</th>
          <td>00:00:00</td>
          <td>80.4</td>
          <td>14.628946</td>
          <td>12.995547</td>
          <td>6.650923</td>
          <td>1.558206</td>
          <td>0.560485</td>
          <td>1.200323</td>
          <td>-0.151276</td>
          <td>22.814208</td>
        </tr>
      </tbody>
    </table>
    <p>1081 rows × 10 columns</p>
    </div>

**Task 3:** Let’s visualize how things look like. Plot the *anomaly* and
*tide* variables.

.. code:: python

    import matplotlib.pyplot as plt
    
    plt.plot(df_sm.anomaly.values, 'k--', label='Anomaly');
    plt.plot(df_sm.tide.values, 'r-', label='Tides');
    
    plt.legend(loc='upper right')
    plt.show()

We get the following plot:

.. image:: plot1.png


Very good! As expected (why?), the observed sea surface height is not
equal to the pure tide signal, but if you zoom in at the peaks (or
throughs), you will see that they are in phase and this is an indicator
that we are on the right track.

**Question**: How much of the observed signal is composed by tides?

To answer this, we will estimate the energy (*J*) contained in both
signals by using the following equation:

.. math:: E = \sum_{n=0}^{N} |F(t)|^2,

where *N* is the number of points of your signal (length) and *F(t)* is
your time series.

**Task 4:** Given the previous equation, we calculate the energy of the
*anomaly* and *tide* signals and the ratio between them. We implement the
function ``energy`` that takes as arguments both signals and *returns* the
:math:`\frac{E_{tide}}{E_{anomaly}}` ratio as percentage.

.. code:: python

    def energy(anomaly, tide):
        e_anom = ((abs(anomaly))**2)
        e_tide = ((abs(tide))**2)
        
        # plot energies
        plt.plot(e_anom, 'k', e_tide,'r')
        
        ratio = np.sum(e_tide) / np.sum(e_anom) # calculate the ratio between the energies
        
        return ratio*100 # in percentage
    
    
    # calculate the ratio
    ratio = energy(df_sm.anomaly.values, df_sm.tide.values)

    plt.show()
    
    print('The ratio is: %f%%\n'%(ratio))

We get the following output:

.. code-block:: text

    The ratio is: 32.158691%

and the following plot:

.. image:: plot2.png


As any other wave, tides also travel (propagate) in the ocean. Their
phase velocity is given by:

.. math:: c = \sqrt{(g*D)}

where *g* is the acceleration of gravity and *D* is the ocean depth.

**Task 5:** Considering *g* = 10 m s\ :math:`^{-2}` and *D* = 4000 m,
how fast (in m s\ :math:`^{-1}`) does this wave travel?

How much time (in hours) does it take to travel from Stavanger to Bodø?
Consider the distance between the two cities as 1200 km.

.. code:: python

    g = 10 # gravity
    D = 4000 # water depth
    L = 1200 # distance, in meters
    
    c = (g*D)**0.5
    print('The phase velocity is %f m/s'%(c))
    
    dt = L/(c*3.6)
    print('The phase velocity is %f hours'%(dt))

We get the following output:

.. code-block:: text

    The phase velocity is 200.000000 m/s
    The phase velocity is 1.666667 hours


This is impressively fast, right? Putting these numbers in words, a
tidal crest observed in Stavanger is also observed in Bodø, 1200 km
away, a bit more than one and a half hours later!