Home › tutorial

From a gauge file to a decomposition

This walks one real-shaped record all the way through: load, inspect, flag, fit, reconstruct, plot, save. Every number printed below comes from the example file that ships with the test data, so you can run the same lines and compare.

1. The record

The example is a single gauge, referred to as station-7 throughout this site: 371 days of six-minute samples starting 1 October 2025, which is 89 040 rows. Water level is in metres, positive upwards, against whatever vertical reference the file declares. The first lines look like this:

# station-7, level in m, reference: station zero, time in UTC
2025-10-01T00:00:00Z,2.014,1
2025-10-01T00:06:00Z,2.067,1
2025-10-01T00:12:00Z,2.118,1
2025-10-01T00:18:00Z,,9

The third column is a quality flag and the fourth row is already a gap: an empty value with flag 9. The reader keeps that row as a NaN rather than dropping it, because a gap at a known time is information and a silently shortened array is not.

2. Loading and looking

import tidalcore as tc

rec = tc.read_gauge("station-7.csv", format="csv-flagged", tz="UTC")
print(rec)

# Record(station='station-7', n=89040, step=6min, span=371.0d,
#        start=2025-10-01T00:00Z, end=2026-10-07T00:00Z,
#        units='m', reference='station zero', missing=978)

Nothing has been computed yet. The reader only parses, checks that the sampling step is consistent, and reports what it found. A record with an irregular step is accepted; step then prints the median interval and the fit uses the actual timestamps.

3. Quality control before anything else

Run the flaggers before the fit, not after. A single 4-metre spike will tilt a least-squares solution quietly and by an amount nobody can guess from the output.

rec = tc.qc.flag_gaps(rec, max_gap="3h")
rec = tc.qc.flag_spikes(rec, window="25h", k=6)
print(tc.qc.summary(rec))

# flag 1  good             87402   98.16 %
# flag 2  probably good      648    0.73 %
# flag 3  probably bad        12    0.01 %
# flag 4  bad                  0    0.00 %
# flag 9  missing            978    1.10 %

The spike detector compares each sample against a running median over a 25-hour window, which is slightly longer than a full diurnal cycle, and flags anything more than k robust deviations away. Twelve samples were caught here; all twelve sit inside a two-hour span on 14 January 2026, which is the kind of cluster that means an instrument event rather than weather.

4. Choosing a constituent set

With 371 days of data, 37 of the 68 catalogued constituents are separable; asking for more does not fail, it just lets pairs of terms trade amplitude against each other. The default picks the resolvable set for you:

print(tc.resolvable(rec.span, catalogue="default"))
# 37 constituents: M2 S2 N2 K2 K1 O1 P1 Q1 ... Sa Ssa Mm Mf MSf

Why 37 and not 68 is a question about record length rather than about software, and it has a clean answer in the Rayleigh table. The short version: separating K1 from P1 needs 182.6 days, so a year-long record gets both, while a 29-day record must infer one from the other.

5. The fit

sol = tc.analyse(rec, constituents="auto", nodal=True)
print(sol.table(top=8))

# name  amplitude_m  phase_deg     snr
#   M2        1.412      211.6   318.4
#   S2        0.463      238.2   104.1
#   N2        0.281      190.8    63.2
#   K2        0.126      235.9    28.4
#   K1        0.109       77.3    24.1
#   O1        0.082       54.0    18.6
#   P1        0.036       74.9     8.2
#   Q1        0.017       31.5     3.9

Amplitudes are in the units of the input, here metres. Phase is the lag in degrees, taken at the reference time of the solution, which defaults to the midpoint of the record and is printed by sol.ref_time. If you compare phases with another analysis, check that reference first; it is the usual cause of a constant offset between two otherwise identical results.

The signal-to-noise column is the ratio of fitted amplitude to its standard error. Terms below about 2 are not evidence of anything, and sol.table() marks them with a trailing asterisk so they do not get copied into a report by accident.

6. Reconstruction and residual

tide = tc.reconstruct(sol, rec.time)
rest = tc.residual(rec, sol)

print(f"observed sd {rec.level.std():.3f} m")
print(f"residual rms {tc.rms(rest):.3f} m")

# observed sd 1.320 m
# residual rms 0.186 m

Those two numbers are the honest summary of the whole exercise: the fit accounts for 98.0 % of the variance in the record, and the remaining 0.186 m root-mean-square is the part no harmonic analysis will ever give you. The largest single excursion in the residual is 0.53 m, on 11 March 2026, during the storm visible in the figure below.

The figure was produced with the plot helper, which needs the plots extra and does nothing clever:

fig = tc.plots.decomposition(rec, sol, start="2026-03-06", days=7)
fig.savefig("station-7-march.png", dpi=150)

7. Filling short gaps

If a downstream tool needs an unbroken series, fill from the solution rather than by interpolation, and keep the mark:

filled = tc.fill_gaps(rec, sol, max_gap="72h")
print((filled.flags == 5).sum())
# 978

Every filled sample carries flag 5, so a later analysis of the filled file will exclude exactly the values the tide invented. Gaps longer than max_gap are left alone on purpose: three days is already long enough that surge and river stage have changed, and a harmonic fill across a week is a forecast wearing a measurement's clothes.

8. Saving the result

sol.to_json("station-7.json", include=("amplitudes", "phases", "errors", "meta"))

The same run from a shell, with the coverage bar made explicit:

tidalcore analyse station-7.csv --format csv-flagged \
    --nodal --min-coverage 0.9 -o station-7.json

9. Before you quote any of this

1. State the record span and the constituent count together. An amplitude without a record length is not a result.
2. State which flags were admitted to the fit, and the coverage that remained after exclusions.
3. Check the reference time before comparing phases with anyone else's numbers.
4. Look at the residual series, not only its root-mean-square. A fit can be excellent on average and wrong for one fortnight.
5. If any term you care about has a signal-to-noise ratio under 2, drop it from the discussion or say plainly that it is at the noise floor.

Going further

Three side pages go deeper than this walk-through on points it only touched:

Next: the API overview lists the return types and exceptions used above, and formats and QC explains the flag column that this record came with.