%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.
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.
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.
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.
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.
print(coef.keys())
tide = utide.reconstruct(obs.index, coef, verbose=False)
The output from the reconstruction is also a Bunch:
print(tide.keys())
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");