Seasonal & Cyclical Alignment for Movement Data

Seasonal and cyclical alignment is the process of mapping multi-year mobility time series onto unified phase-aware coordinates so that behavioral rhythms — commute peaks, freight cycles, tourism surges — can be compared directly across calendar years, despite leap years, shifting event schedules, and long-term structural drift.

Without it, year-over-year comparisons in Temporal Aggregation & Window Mapping pipelines suffer from phase misalignment, rolling windows capture artificial discontinuities, and forecasting models overfit to calendar artifacts rather than true demand patterns.

Prerequisites Checklist

Complete these before implementing the alignment pipeline. This section sits downstream of ingestion and timezone normalisation inside Temporal Aggregation & Window Mapping.

• Python packages: pandas >= 2.0, numpy, scipy, statsmodels >= 0.14, and optionally tslearn >= 0.6 for DTW.
• Data schema: A DatetimeIndex-aligned Series or DataFrame with columns entity_id, timestamp (UTC, timezone-aware), and a numeric demand signal (counts, volume, or dwell_seconds). At least two full annual cycles are required — three or more improve STL robustness.
• Sampling consistency: Uniform intervals (15 min, hourly, or daily) are required before decomposition. Irregular sampling must be resolved first via gap-filling in sparse trajectories, including Kalman-based interpolation for sub-minute GPS traces.
• Spatial pre-aggregation: Trajectories should already be aggregated to consistent spatial zones (H3 hexagons, census tracts, or transit corridors) via coordinate reference system mapping before temporal alignment begins — mixing raw coordinate precision with temporal decomposition inflates residual noise.
• Upstream validation: Timestamps must be free of DST boundary artifacts. Apply pd.date_range(freq='h', tz='UTC') normalization and verify with time-series synchronization strategies before proceeding.

Failure Mode Taxonomy

Deterministic Pipeline Overview

The alignment pipeline has five ordered stages. Each stage produces a validated intermediate that the next stage consumes.

1. Timezone canonicalisation — convert all timestamps to UTC; verify no naive datetimes remain.
2. Uniform resampling — resample or interpolate to a consistent interval; record the fill-rate per zone.
3. Circular temporal encoding — generate sine/cosine pairs for day-of-year, day-of-week, and hour-of-day.
4. STL decomposition — decompose each zone’s series into trend, seasonal, and residual components; validate stationarity of residuals.
5. DTW phase alignment — align the isolated seasonal component of each target year to a reference year using constrained dynamic time warping; output a phase-corrected series.

Implementation Walkthrough

from __future__ import annotations

import numpy as np
import pandas as pd
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.stattools import adfuller
from typing import Optional


def encode_circular(
    series: pd.Series,
    period: float,
) -> tuple[pd.Series, pd.Series]:
    """Return sine and cosine encodings of `series` normalised to `period`.

    Args:
        series: Numeric values (e.g. day-of-year 1..366, hour 0..23).
        period: Full cycle length (365.25 for annual, 7 for weekly, 24 for daily).

    Returns:
        Tuple of (sin_enc, cos_enc) Series with the same index as `series`.
    """
    if series.empty:
        return series.copy(), series.copy()
    angle = 2.0 * np.pi * series / period
    return np.sin(angle), np.cos(angle)


def decompose_stl(
    ts: pd.Series,
    period: int,
    robust: bool = True,
) -> tuple[pd.Series, pd.Series, pd.Series]:
    """Decompose a time series into trend, seasonal, and residual components.

    The series must have a uniform DatetimeIndex and no NaN values.
    Uses statsmodels STL with LOESS smoothing.

    Args:
        ts: Uniformly-sampled numeric Series with DatetimeIndex.
        period: Seasonal period in samples (e.g. 365 for daily annual).
        robust: Whether to use iteratively-reweighted least-squares (recommended
                for mobility data with outlier events such as public holidays).

    Returns:
        Tuple of (trend, seasonal, residual) Series.

    Raises:
        ValueError: If `ts` is shorter than 2 × period (insufficient data).
    """
    if len(ts) < 2 * period:
        raise ValueError(
            f"Series length {len(ts)} is less than 2 × period ({2 * period}). "
            "Provide at least two full seasonal cycles."
        )
    if ts.isna().any():
        raise ValueError("ts must not contain NaN values. Run gap-filling first.")

    result = STL(ts, period=period, robust=robust).fit()
    return result.trend, result.seasonal, result.resid


def align_with_dtw(
    reference: np.ndarray,
    target: np.ndarray,
    sakoe_chiba_radius: Optional[int] = 14,
) -> np.ndarray:
    """Align `target` seasonal profile to `reference` using constrained DTW.

    Requires tslearn. The Sakoe-Chiba band prevents unrealistic temporal
    stretching — keep the radius at ≤10% of the series length.

    Args:
        reference: 1-D array of the reference seasonal profile (e.g. year N).
        target: 1-D array to be aligned (e.g. year N+1). Same length as reference.
        sakoe_chiba_radius: Maximum warping distance in samples. None = unconstrained
                            (not recommended for production use).

    Returns:
        Aligned version of `target` resampled onto `reference`'s time axis.

    Raises:
        ImportError: If tslearn is not installed.
        ValueError: If reference and target have different lengths.
    """
    try:
        from tslearn.metrics import dtw_path
    except ImportError as exc:
        raise ImportError("Install tslearn>=0.6: pip install tslearn") from exc

    if len(reference) != len(target):
        raise ValueError(
            f"reference length ({len(reference)}) != target length ({len(target)}). "
            "Both series must span the same number of time steps."
        )

    r = reference.reshape(-1, 1)
    t = target.reshape(-1, 1)
    path, _ = dtw_path(r, t, sakoe_chiba_radius=sakoe_chiba_radius)

    # Resample target onto the reference time axis via the optimal warp path
    aligned = np.zeros(len(reference))
    counts = np.zeros(len(reference), dtype=int)
    for ref_idx, tgt_idx in path:
        aligned[ref_idx] += target[tgt_idx]
        counts[ref_idx] += 1
    counts = np.where(counts == 0, 1, counts)  # guard against unmapped positions
    return aligned / counts


def validate_alignment(
    reference_seasonal: pd.Series,
    aligned_seasonal: pd.Series,
    original_ts: pd.Series,
) -> dict[str, float]:
    """Run post-alignment sanity checks.

    Args:
        reference_seasonal: Seasonal component of the reference year.
        aligned_seasonal: Seasonal component after DTW alignment.
        original_ts: The full original time series (for variance partitioning).

    Returns:
        Dict with keys: max_lag_offset, seasonal_variance_ratio, adf_pvalue.
    """
    cross_corr = np.correlate(
        reference_seasonal - reference_seasonal.mean(),
        aligned_seasonal - aligned_seasonal.mean(),
        mode="full",
    )
    lags = np.arange(-len(reference_seasonal) + 1, len(reference_seasonal))
    max_lag = int(lags[np.argmax(cross_corr)])

    seasonal_var = float(np.var(aligned_seasonal))
    total_var = float(np.var(original_ts))
    seasonal_ratio = seasonal_var / total_var if total_var > 0 else 0.0

    # Stationarity test on residual (lower p-value = more stationary = better)
    _, adf_p, *_ = adfuller(original_ts - aligned_seasonal, autolag="AIC")

    return {
        "max_lag_offset": max_lag,
        "seasonal_variance_ratio": seasonal_ratio,
        "adf_pvalue": float(adf_p),
    }

End-to-end usage

import pandas as pd
import numpy as np

# 1. Load and canonicalise timezone
raw = pd.read_parquet("fleet_daily_counts.parquet")
raw["timestamp"] = pd.to_datetime(raw["timestamp"], utc=True)
raw = raw.sort_values("timestamp").set_index("timestamp")

ts = raw["trip_count"].resample("D").sum().fillna(0)  # daily totals, no gaps

# 2. Circular encoding (store alongside features for ML downstream)
doy = ts.index.day_of_year.astype(float)
sin_doy, cos_doy = encode_circular(pd.Series(doy, index=ts.index), period=365.25)

# 3. STL decomposition — annual seasonality on a daily series
trend, seasonal, resid = decompose_stl(ts, period=365, robust=True)

# 4. Split into reference year (2022) and target year (2023)
ref_seasonal = seasonal["2022"].values
tgt_seasonal = seasonal["2023"].values

# Pad/trim to equal length if years differ by a day
min_len = min(len(ref_seasonal), len(tgt_seasonal))
ref_seasonal = ref_seasonal[:min_len]
tgt_seasonal = tgt_seasonal[:min_len]

# 5. DTW alignment with a ±14-day Sakoe-Chiba band
aligned = align_with_dtw(ref_seasonal, tgt_seasonal, sakoe_chiba_radius=14)

# 6. Validate
aligned_series = pd.Series(aligned, index=seasonal["2023"].index[:min_len])
metrics = validate_alignment(
    pd.Series(ref_seasonal), aligned_series, ts["2023"].iloc[:min_len]
)
print(metrics)
# Expected: max_lag_offset ≈ 0, seasonal_variance_ratio > 0.6, adf_pvalue < 0.05

Mathematical Grounding

Circular encoding. A scalar day-of-year $d \in [1, 366]$ is embedded on the unit circle as:

$$\phi_s = \sin!\left(\frac{2\pi d}{365.25}\right), \quad \phi_c = \cos!\left(\frac{2\pi d}{365.25}\right)$$

Both components are required. Using only $\phi_s$ loses quadrant information — day 90 and day 275 (symmetric about the summer solstice) become indistinguishable. The pair $(\phi_s, \phi_c)$ uniquely identifies any day within the year and preserves angular proximity, so that the Euclidean distance $|\phi_{\text{Dec 30}} - \phi_{\text{Jan 2}}|$ is small rather than near-maximal.

STL decomposition. STL (Seasonal-Trend decomposition using Loess) iterates between LOESS smoothing of the seasonal subseries and a robust trend estimator. The decomposition is:

$$y_t = T_t + S_t + R_t$$

where $T_t$ is the trend, $S_t$ the seasonal component (repeating with period $p$), and $R_t$ the residual. The robust=True flag down-weights large residuals, preventing public holidays or pandemic-era anomalies from distorting the extracted seasonal shape.

Dynamic Time Warping. DTW finds the monotone index mapping $\pi: [1, N] \to [1, M]$ that minimises cumulative pointwise distance $\sum_{(i,j) \in \pi} |x_i - y_j|^2$. The Sakoe-Chiba constraint $|i - j| \leq r$ limits the warping to a diagonal band of half-width $r$, preventing the algorithm from matching, say, a February demand trough to an August one.

Calibration & Parameter Tuning

For daily resolution data with an annual period, the seasonal_deg (polynomial degree of seasonal LOESS) should stay at 0 or 1. Higher degrees allow the seasonal shape to evolve over time, which is useful for decade-long series but can absorb genuine year-specific shocks you want to preserve in the residual.

Integration & Compatibility

The aligned seasonal series slot directly into adjacent pipeline stages within Temporal Aggregation & Window Mapping:

• Rolling Statistics for Mobility Metrics: Feed phase-aligned series as the input to sliding-window computations. Because alignment removes calendar-induced variance, rolling averages and volatility bands now capture genuine behavioral shifts rather than seasonal artifacts.
• Dynamic Time-Binning Strategies: Use the extracted seasonal profile to set adaptive bin boundaries — narrow bins near seasonal peaks, wider bins during off-season troughs — rather than applying uniform intervals.
• Gap-filling upstream: Any NaN values from sparse sensor coverage must be resolved via gap-filling in sparse trajectories before STL decomposition. STL will raise a convergence warning on series with more than ~5% missing values even when missing='drop' is used.
• Downstream forecasting models: The (trend, sin_doy, cos_doy) feature set produced here integrates cleanly with LightGBM, XGBoost, or Facebook Prophet as pre-engineered temporal features. Pass the aligned seasonal component as a regressor to Prophet via add_regressor().
• Stay-point analysis: When aligning dwell-time series (e.g., average stop duration at transit nodes), combine this pipeline with stay-point detection algorithms to ensure only genuine stops contribute to the seasonal signal.

Production Implementation Notes

1. Vectorize over iteration. Avoid apply() or Python loops for circular encoding and resampling. Use numpy broadcasting and pandas vectorized datetime accessors (.dt.day_of_year, .dt.isocalendar()).
2. Timezone discipline. Store all raw timestamps in UTC. Convert to local time only for visualization or reporting. Mixing naive and aware timestamps causes silent misalignment during DST transitions.
3. Memory management. Multi-year, high-frequency mobility datasets easily exceed RAM limits. Process data in spatially partitioned chunks (e.g. by H3 parent cell or transit corridor) and use dask.dataframe for out-of-core STL decomposition.
4. Reproducibility. STL and DTW implementations may use stochastic smoothing. Set np.random.seed() and pin statsmodels and tslearn versions in your environment lockfile.
5. Pandas 2.x frequency strings. The uppercase 'H' hourly alias was deprecated in pandas 2.2; use lowercase 'h'. Similarly, 'T' for minutes is replaced by 'min'. Run pd.tseries.frequencies.to_offset('h') to verify your environment before deploying.

FAQ

Why use both sin and cos for encoding — can I just use sin? Using only sin makes day 90 and day 275 equidistant from the origin, giving them identical encoded values despite being six months apart. The cosine component breaks this symmetry: the two-component pair $(\sin\theta, \cos\theta)$ is a bijection from $[0, 2\pi)$ to the unit circle, so each day receives a unique encoding.

What does an adf_pvalue above 0.05 in the validation output mean? A p-value above 0.05 from the Augmented Dickey-Fuller test indicates the residual series is non-stationary — there is likely a remaining trend or low-frequency cycle that STL did not capture. Increase the seasonal_deg parameter, extend your series to cover more cycles, or apply a second STL pass at the weekly frequency.

Can the same pipeline handle both daily and sub-hourly data? Yes, but adjust period to match the resolution. For 15-minute data with a weekly cycle, use period=672 (7 days × 24 h × 4 intervals/h). Memory usage scales linearly with series length, so partition by spatial unit before running STL on sub-hourly series.

How do I align more than two years at once? Designate one year as the reference (typically the most complete or representative year). Align all other years pairwise against that reference using the same Sakoe-Chiba radius. Do not chain alignments (year 1→2, year 2→3) as cumulative warp errors compound.

Does this approach work for non-annual cycles like tidal or event-driven rhythms? Yes. Set period to the cycle length in samples. For tidal rhythms at hourly resolution, period=25 (approximate M2 tidal period ≈ 12.4 h, so 24.8 samples per cycle rounded to 25). For event-driven cycles with irregular spacing, DTW alignment is more appropriate than STL; use it directly on normalized demand profiles around event windows.

Related

• Aligning Seasonal Travel Patterns Across Multiple Years
• Rolling Statistics for Mobility Metrics
• Dynamic Time-Binning Strategies
• Gap-Filling in Sparse Trajectories
• Time-Series Synchronization Strategies

Back to Temporal Aggregation & Window Mapping