#!/usr/bin/env python
# coding: utf-8

# In[ ]:


get_ipython().run_line_magic('matplotlib', 'inline')

import datetime

import numpy as np
import matplotlib.pyplot as plt

import pandas as pd

import utide

print(utide.__version__)


# Look at the data file to see what structure it has.

# In[ ]:


with open("can1998.dtf") as f:
    lines = f.readlines()

print("".join(lines[:5]))


# It looks like the fields are seconds, year, month, day, hour, elevation, flag.  We need a date parser function to combine the date and time fields into a single value to be used as the datetime index.

# In[ ]:


names = ["seconds", "year", "month", "day", "hour", "elev", "flag"]

obs = pd.read_csv(
    "can1998.dtf",
    names=names,
    skipinitialspace=True,
    delim_whitespace=True,
    na_values="9.990",
)


date_cols = ["year", "month", "day", "hour"]
index = pd.to_datetime(obs[date_cols])
obs = obs.drop(date_cols, axis=1)
obs.index = index

obs.head(5)


# Although there are no elevations marked bad via special value, which should be `nan` after reading the file, the flag value of 2 indicates the values are unreliable, so we will mark them with `nan`, calculate the deviations of the elevations from their mean (stored in a new column called "anomaly"), and then interpolate to fill in the `nan` values in the anomaly.

# In[ ]:


bad = obs["flag"] == 2
corrected = obs["flag"] == 1

obs.loc[bad, "elev"] = np.nan
obs["anomaly"] = obs["elev"] - obs["elev"].mean()
obs["anomaly"] = obs["anomaly"].interpolate()
print(f"{bad.sum()} points were flagged 'bad' and interpolated")
print(f"{corrected.sum()} points were flagged 'corrected' and left unchanged")


# Now we can call solve to obtain the coefficients.

# In[ ]:


coef = utide.solve(
    obs.index,
    obs["anomaly"],
    lat=-25,
    method="ols",
    conf_int="MC",
    verbose=False,
)


# The amplitudes and phases from the fit are now in the `coef` data structure (a Bunch), which can be used directly in the `reconstruct` function to generate a hindcast or forecast of the tides at the times specified in the `time` array.

# In[ ]:


print(coef.keys())


# In[ ]:


tide = utide.reconstruct(obs.index, coef, verbose=False)


# The output from the reconstruction is also a Bunch:

# In[ ]:


print(tide.keys())


# In[ ]:


t = obs.index.to_pydatetime()

fig, (ax0, ax1, ax2) = plt.subplots(figsize=(17, 5), nrows=3, sharey=True, sharex=True)

ax0.plot(t, obs.anomaly, label="Observations", color="C0")
ax1.plot(t, tide.h, label="Prediction", color="C1")
ax2.plot(t, obs.anomaly - tide.h, label="Residual", color="C2")
fig.legend(ncol=3, loc="upper center");