Skip to content

rivabar.prediction

rivabar.prediction

Migration prediction using the Howard and Knutson (1984) meander migration model.

Implements the convolution-based approach described in Sylvester et al. (2019, Geology), where the predicted migration rate is computed from curvature using an exponential downstream weighting function.

Key equations

Nominal migration rate (Eq. 2 in the supplementary material): R0(s) = kl * W * C(s)

Adjusted (predicted) migration rate (Eq. 3): R1(s) = Omega * R0(s) + Gamma * (integral of R0(s-xi)*G(xi) dxi) / (integral of G(xi) dxi)

where G(xi) = exp(-2 * Cf * xi / D) is the exponential weighting function, D is river depth (estimated from width), and Cf is a friction factor.

Depth from width (Eq. 5, from Konsoer et al., 2013): W = 18.8 * D^1.41 => D = (W / 18.8)^(1/1.41)

References

Howard, A.D. and Knutson, T.R., 1984. Sufficient conditions for river meandering: A simulation approach. Water Resources Research, 20, p.1659-1667. Sylvester, Z., Durkin, P. and Covault, J.A., 2019. High curvatures drive river meandering. Geology, 47(3), p.263-266.

depth_from_width(W)

Estimate river depth from channel width using the Konsoer et al. (2013) regression.

W = 18.8 * D^1.41 => D = (W / 18.8)^(1/1.41)

Parameters:

Name Type Description Default
W float or array_like

Channel width in metres.

 required

Returns:

Name Type Description
D float or ndarray

Estimated depth in metres.
Source code in rivabar/prediction.py 
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def depth_from_width(W):
    """
    Estimate river depth from channel width using the Konsoer et al. (2013)
    regression.

    W = 18.8 * D^1.41  =>  D = (W / 18.8)^(1/1.41)

    Parameters
    ----------
    W : float or array_like
        Channel width in metres.

    Returns
    -------
    D : float or ndarray
        Estimated depth in metres.
    """
    return (np.asarray(W) / 18.8) ** (1.0 / 1.41)

nominal_migration_rate(curvature, W, kl)

Compute the nominal migration rate R0 = kl * W * C.

Parameters:

Name Type Description Default
curvature array_like

Curvature along the channel (1/m).

 required
W float

Mean channel width in metres.

 required
kl float

Migration rate constant (erodibility), in m/yr.

 required

Returns:

Name Type Description
R0 ndarray

Nominal migration rate (m/yr).
Source code in rivabar/prediction.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def nominal_migration_rate(curvature, W, kl):
    """
    Compute the nominal migration rate R0 = kl * W * C.

    Parameters
    ----------
    curvature : array_like
        Curvature along the channel (1/m).
    W : float
        Mean channel width in metres.
    kl : float
        Migration rate constant (erodibility), in m/yr.

    Returns
    -------
    R0 : ndarray
        Nominal migration rate (m/yr).
    """
    return -kl * W * np.asarray(curvature)

predicted_migration_rate(R0, s, D, Cf, Omega=-1.0, Gamma=2.5)

Compute the predicted migration rate using the Howard and Knutson (1984) convolution model.

R1(s) = Omega * R0(s) + Gamma * conv(R0, G) / integral(G)

where G(xi) = exp(-2 * Cf * xi / D) and the convolution integrates upstream influence (causal, one-sided).

Parameters:

Name Type Description Default
R0 array_like

Nominal migration rate along the channel.

 required
s array_like

Along-channel distance (m), must be uniformly spaced.

 required
D float

River depth in metres.

 required
Cf float

Friction factor (dimensionless); typical range 0.003-0.009.

 required
Omega float

Weight for the local contribution (default -1.0).

 -1.0
Gamma float

Weight for the upstream-integrated contribution (default 2.5). Note: Omega + Gamma = 1.5 for the standard HK model, so that R1 = 1.5 * R0 for a bend with constant curvature.

 2.5

Returns:

Name Type Description
R1 ndarray

Predicted migration rate (m/yr), same length as R0.
Source code in rivabar/prediction.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def predicted_migration_rate(R0, s, D, Cf, Omega=-1.0, Gamma=2.5):
    """
    Compute the predicted migration rate using the Howard and Knutson (1984)
    convolution model.

    R1(s) = Omega * R0(s) + Gamma * conv(R0, G) / integral(G)

    where G(xi) = exp(-2 * Cf * xi / D) and the convolution integrates
    upstream influence (causal, one-sided).

    Parameters
    ----------
    R0 : array_like
        Nominal migration rate along the channel.
    s : array_like
        Along-channel distance (m), must be uniformly spaced.
    D : float
        River depth in metres.
    Cf : float
        Friction factor (dimensionless); typical range 0.003-0.009.
    Omega : float, optional
        Weight for the local contribution (default -1.0).
    Gamma : float, optional
        Weight for the upstream-integrated contribution (default 2.5).
        Note: Omega + Gamma = 1.5 for the standard HK model, so that
        R1 = 1.5 * R0 for a bend with constant curvature.

    Returns
    -------
    R1 : ndarray
        Predicted migration rate (m/yr), same length as R0.
    """
    R0 = np.asarray(R0, dtype=float)
    s = np.asarray(s, dtype=float)
    ds = np.mean(np.diff(s))
    n = len(R0)

    # Build exponential kernel G(xi) = exp(-2*Cf*xi/D)
    # Truncate at 10 * D/(2*Cf) where the kernel is negligible
    if Cf > 0 and D > 0:
        decay_length = D / (2.0 * Cf)
        max_xi = min(10.0 * decay_length, (n - 1) * ds)
        n_kernel = max(int(np.ceil(max_xi / ds)), 1)
        xi = np.arange(n_kernel) * ds
        G = np.exp(-2.0 * Cf * xi / D)
    else:
        # No downstream influence: R1 = (Omega + Gamma) * R0
        return (Omega + Gamma) * R0

    # Normalized convolution: sum(R0 * G) / sum(G), where the denominator at
    # each point is the kernel mass actually available upstream of that point
    # (a true weighted average; near the upstream boundary fewer points exist,
    # so normalizing by the full kernel integral would attenuate the term and
    # flip the sign of R1 there)
    conv_full = np.convolve(R0, G, mode='full')[:n] * ds
    norm = np.convolve(np.ones(n), G, mode='full')[:n] * ds
    upstream_term = conv_full / norm

    R1 = Omega * R0 + Gamma * upstream_term
    return R1

calibrate_Cf(R0, observed_mr, s, D, Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5, mask=None)

Optimize the friction factor Cf to minimize the phase shift between predicted and observed migration rates.

Parameters:

Name Type Description Default
R0 array_like

Nominal migration rate (from kl * W * C).

 required
observed_mr array_like

Observed migration rate (m/yr).

 required
s array_like

Along-channel distance (m).

 required
D float

Estimated river depth (m).

 required
Cf_range tuple

(min, max) bounds for Cf search (default (0.001, 0.02)).

 (0.001, 0.02)
Omega float

Local weight (default -1.0).

 -1.0
Gamma float

Upstream weight (default 2.5).

 2.5
mask array_like of bool

Points to use when scoring the correlation. The prediction is always computed on the full contiguous arrays (the convolution requires uniform spacing); the mask only restricts the statistics.

 None

Returns:

Name Type Description
Cf_opt float

Optimal friction factor.
r_opt float

Pearson correlation at optimal Cf.
Source code in rivabar/prediction.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
def calibrate_Cf(R0, observed_mr, s, D, Cf_range=(0.001, 0.02),
                 Omega=-1.0, Gamma=2.5, mask=None):
    """
    Optimize the friction factor Cf to minimize the phase shift between
    predicted and observed migration rates.

    Parameters
    ----------
    R0 : array_like
        Nominal migration rate (from kl * W * C).
    observed_mr : array_like
        Observed migration rate (m/yr).
    s : array_like
        Along-channel distance (m).
    D : float
        Estimated river depth (m).
    Cf_range : tuple, optional
        (min, max) bounds for Cf search (default (0.001, 0.02)).
    Omega : float, optional
        Local weight (default -1.0).
    Gamma : float, optional
        Upstream weight (default 2.5).
    mask : array_like of bool, optional
        Points to use when scoring the correlation. The prediction is always
        computed on the full contiguous arrays (the convolution requires
        uniform spacing); the mask only restricts the statistics.

    Returns
    -------
    Cf_opt : float
        Optimal friction factor.
    r_opt : float
        Pearson correlation at optimal Cf.
    """
    R0 = np.asarray(R0, dtype=float)
    observed_mr = np.asarray(observed_mr, dtype=float)
    if mask is None:
        mask = np.ones(len(R0), dtype=bool)

    def neg_correlation(Cf):
        R1 = predicted_migration_rate(R0, s, D, Cf, Omega, Gamma)
        r = np.corrcoef(R1[mask], observed_mr[mask])[0, 1]
        return -r

    result = minimize_scalar(neg_correlation, bounds=Cf_range, method='bounded')
    Cf_opt = result.x
    r_opt = -result.fun
    return Cf_opt, r_opt

calibrate_kl(predicted, observed, percentile=75)

Estimate the migration rate constant kl by minimizing the amplitude mismatch between predicted and observed migration rates.

Uses the ratio of the specified percentile of absolute values, following the approach in Sylvester et al. (2019).

Parameters:

Name Type Description Default
predicted array_like

Predicted migration rate (computed with kl=1).

 required
observed array_like

Observed migration rate (m/yr).

 required
percentile float

Percentile of absolute values to match (default 75).

 75

Returns:

Name Type Description
kl_scale float

Scale factor to apply to kl. If the initial kl was kl_init, the refined kl is kl_init * kl_scale.
Source code in rivabar/prediction.py 
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def calibrate_kl(predicted, observed, percentile=75):
    """
    Estimate the migration rate constant kl by minimizing the amplitude
    mismatch between predicted and observed migration rates.

    Uses the ratio of the specified percentile of absolute values, following
    the approach in Sylvester et al. (2019).

    Parameters
    ----------
    predicted : array_like
        Predicted migration rate (computed with kl=1).
    observed : array_like
        Observed migration rate (m/yr).
    percentile : float, optional
        Percentile of absolute values to match (default 75).

    Returns
    -------
    kl_scale : float
        Scale factor to apply to kl. If the initial kl was kl_init,
        the refined kl is kl_init * kl_scale.
    """
    predicted = np.asarray(predicted, dtype=float)
    observed = np.asarray(observed, dtype=float)
    p_pred = np.percentile(np.abs(predicted), percentile)
    p_obs = np.percentile(np.abs(observed), percentile)
    if p_pred > 0:
        return p_obs / p_pred
    return 1.0

calibrate_from_curvature(curvature, observed_mr, s, W, Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5, kl_percentile=75, mask=None)

Full calibration of D, Cf, and kl from curvature and observed migration rate for a single pair.

Procedure (following Sylvester et al., 2019): 1. Estimate D from W using the Konsoer et al. (2013) regression. 2. Compute initial kl from R1 = 1.5 * R0 simplification. 3. Optimize Cf to minimize the phase shift. 4. Refine kl to match the amplitude.

Parameters:

Name Type Description Default
curvature array_like

Curvature along the channel (1/m).

 required
observed_mr array_like

Observed migration rate (m/yr).

 required
s array_like

Along-channel distance (m).

 required
W float

Mean channel width (m).

 required
Cf_range tuple

Bounds for Cf optimization (default (0.001, 0.02)).

 (0.001, 0.02)
Omega float

Local weight (default -1.0).

 -1.0
Gamma float

Upstream weight (default 2.5).

 2.5
kl_percentile float

Percentile for amplitude matching (default 75).

 75
mask array_like of bool

Points to use for the calibration statistics (correlations and percentiles). The prediction itself is always computed on the full contiguous arrays, since the convolution assumes uniform spacing.

 None

Returns:

Name Type Description
result dict

Calibration results with keys: - 'D': estimated depth (m) - 'Cf': optimized friction factor - 'kl': calibrated migration rate constant (m/yr) - 'W': mean channel width (m) - 'Omega': local weight used - 'Gamma': upstream weight used - 'r_nominal': correlation between R0 and observed - 'r_predicted': correlation between R1 and observed - 'R0': nominal migration rate array - 'R1': predicted migration rate array - 'observed_mr': observed migration rate array - 's': along-channel distance array
Source code in rivabar/prediction.py 
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
def calibrate_from_curvature(curvature, observed_mr, s, W,
                             Cf_range=(0.001, 0.02),
                             Omega=-1.0, Gamma=2.5,
                             kl_percentile=75, mask=None):
    """
    Full calibration of D, Cf, and kl from curvature and observed migration
    rate for a single pair.

    Procedure (following Sylvester et al., 2019):
    1. Estimate D from W using the Konsoer et al. (2013) regression.
    2. Compute initial kl from R1 = 1.5 * R0 simplification.
    3. Optimize Cf to minimize the phase shift.
    4. Refine kl to match the amplitude.

    Parameters
    ----------
    curvature : array_like
        Curvature along the channel (1/m).
    observed_mr : array_like
        Observed migration rate (m/yr).
    s : array_like
        Along-channel distance (m).
    W : float
        Mean channel width (m).
    Cf_range : tuple, optional
        Bounds for Cf optimization (default (0.001, 0.02)).
    Omega : float, optional
        Local weight (default -1.0).
    Gamma : float, optional
        Upstream weight (default 2.5).
    kl_percentile : float, optional
        Percentile for amplitude matching (default 75).
    mask : array_like of bool, optional
        Points to use for the calibration statistics (correlations and
        percentiles). The prediction itself is always computed on the full
        contiguous arrays, since the convolution assumes uniform spacing.

    Returns
    -------
    result : dict
        Calibration results with keys:
        - 'D': estimated depth (m)
        - 'Cf': optimized friction factor
        - 'kl': calibrated migration rate constant (m/yr)
        - 'W': mean channel width (m)
        - 'Omega': local weight used
        - 'Gamma': upstream weight used
        - 'r_nominal': correlation between R0 and observed
        - 'r_predicted': correlation between R1 and observed
        - 'R0': nominal migration rate array
        - 'R1': predicted migration rate array
        - 'observed_mr': observed migration rate array
        - 's': along-channel distance array
    """
    curvature = np.asarray(curvature, dtype=float)
    observed_mr = np.asarray(observed_mr, dtype=float)
    s = np.asarray(s, dtype=float)
    if mask is None:
        mask = np.ones(len(curvature), dtype=bool)

    # Step 1: estimate depth
    D = depth_from_width(W)

    # Step 2: initial kl from R1 = 1.5 * R0
    # For constant curvature: R1 = (Omega + Gamma) * R0 = 1.5 * kl * W * C
    # Match the 75th percentile of |observed_mr| to 1.5 * kl * W * |C|_p75
    factor = Omega + Gamma  # should be 1.5
    p_obs = np.percentile(np.abs(observed_mr[mask]), kl_percentile)
    p_curv = np.percentile(np.abs(curvature[mask]), kl_percentile)
    if p_curv > 0 and W > 0 and abs(factor) > 0:
        kl_init = p_obs / (abs(factor) * W * p_curv)
    else:
        kl_init = 1.0

    # Compute nominal migration rate
    R0 = nominal_migration_rate(curvature, W, kl_init)
    r_nominal = np.corrcoef(R0[mask], observed_mr[mask])[0, 1]

    # Step 3: optimize Cf
    Cf_opt, r_predicted = calibrate_Cf(R0, observed_mr, s, D, Cf_range,
                                       Omega, Gamma, mask=mask)

    # Step 4: refine kl
    R1_unscaled = predicted_migration_rate(R0, s, D, Cf_opt, Omega, Gamma)
    kl_scale = calibrate_kl(R1_unscaled[mask], observed_mr[mask], kl_percentile)
    kl_final = kl_init * kl_scale

    # Recompute with final kl
    R0_final = nominal_migration_rate(curvature, W, kl_final)
    R1_final = predicted_migration_rate(R0_final, s, D, Cf_opt, Omega, Gamma)
    r_predicted_final = np.corrcoef(R1_final[mask], observed_mr[mask])[0, 1]

    return {
        'D': D,
        'Cf': Cf_opt,
        'kl': kl_final,
        'W': W,
        'Omega': Omega,
        'Gamma': Gamma,
        'r_nominal': r_nominal,
        'r_predicted': r_predicted_final,
        'R0': R0_final,
        'R1': R1_final,
        'observed_mr': observed_mr,
        's': s,
    }

calibrate_pair(results, pair_idx, Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5, kl_percentile=75, use_stable_segments=True)

Calibrate the HK model for a single pair from analysis results.

Parameters:

Name Type Description Default
results dict

Results dictionary from analyze_river_pairs_filtered or analyze_segment_group.

 required
pair_idx int

Index of the pair to calibrate.

 required
Cf_range tuple

Bounds for Cf optimization (default (0.001, 0.02)).

 (0.001, 0.02)
Omega float

Local weight (default -1.0).

 -1.0
Gamma float

Upstream weight (default 2.5).

 2.5
kl_percentile float

Percentile for amplitude matching (default 75).

 75
use_stable_segments bool

If True, only use stable segments (exclude high-variance regions) for calibration. The prediction is still computed everywhere. Default True.

 True

Returns:

Name Type Description
result dict

Calibration results (see calibrate_from_curvature), plus: - 'date1', 'date2': scene dates - 'time_gap_years': time gap between scenes - 'pair_idx': index in the results dict
Source code in rivabar/prediction.py 
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
def calibrate_pair(results, pair_idx, Cf_range=(0.001, 0.02),
                   Omega=-1.0, Gamma=2.5, kl_percentile=75,
                   use_stable_segments=True):
    """
    Calibrate the HK model for a single pair from analysis results.

    Parameters
    ----------
    results : dict
        Results dictionary from ``analyze_river_pairs_filtered`` or
        ``analyze_segment_group``.
    pair_idx : int
        Index of the pair to calibrate.
    Cf_range : tuple, optional
        Bounds for Cf optimization (default (0.001, 0.02)).
    Omega : float, optional
        Local weight (default -1.0).
    Gamma : float, optional
        Upstream weight (default 2.5).
    kl_percentile : float, optional
        Percentile for amplitude matching (default 75).
    use_stable_segments : bool, optional
        If True, only use stable segments (exclude high-variance regions)
        for calibration. The prediction is still computed everywhere.
        Default True.

    Returns
    -------
    result : dict
        Calibration results (see ``calibrate_from_curvature``), plus:
        - 'date1', 'date2': scene dates
        - 'time_gap_years': time gap between scenes
        - 'pair_idx': index in the results dict
    """
    pair_info = results['pair_info'][pair_idx]
    curvature = results['curvatures'][pair_idx]
    distances = results['migration_distances'][pair_idx]
    s = results['along_channel_distances'][pair_idx]
    time_gap_years = pair_info['time_gap_years']
    W = np.mean(pair_info['width1'])

    # Convert distances to migration rate
    observed_mr = distances / time_gap_years

    # Determine which points to use for the calibration statistics. The
    # prediction is always computed on the full contiguous arrays: slicing
    # out high-variance regions before the convolution would concatenate
    # spatially disjoint points and distort the upstream-influence kernel.
    mask = _stable_segment_mask(pair_info, len(curvature)) if use_stable_segments else None

    # Calibrate; statistics restricted to stable segments via the mask
    cal = calibrate_from_curvature(curvature, observed_mr, s, W,
                                   Cf_range, Omega, Gamma, kl_percentile,
                                   mask=mask)

    cal['date1'] = pair_info['date1']
    cal['date2'] = pair_info['date2']
    cal['time_gap_years'] = time_gap_years
    cal['pair_idx'] = pair_idx

    # Residual standard deviation (on stable segments only)
    residuals = observed_mr - cal['R1']
    if mask is not None:
        cal['residual_std'] = np.std(residuals[mask])
    else:
        cal['residual_std'] = np.std(residuals)

    return cal

calibrate_segment(results, df, pair_class_filter='good', Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5, kl_percentile=75, use_stable_segments=True, min_time_gap_years=0, pair_indices=None)

Calibrate the HK model for all qualifying pairs in a segment.

Parameters:

Name Type Description Default
results dict

Results dictionary from analyze_river_pairs_filtered or analyze_segment_group.

 required
df DataFrame

DataFrame from create_dataframe_from_results, used to filter pairs by classification.

 required
pair_class_filter str or list of str

Only calibrate pairs with these classifications (default 'good').

 'good'
Cf_range tuple

Bounds for Cf optimization.

 (0.001, 0.02)
Omega float

Local weight (default -1.0).

 -1.0
Gamma float

Upstream weight (default 2.5).

 2.5
kl_percentile float

Percentile for amplitude matching (default 75).

 75
use_stable_segments bool

Whether to exclude high-variance segments during calibration (default True).

 True
min_time_gap_years float

Minimum time gap (in years) for a pair to be included in calibration (default 0). Short-gap pairs tend to have noisier migration estimates, which can inflate kl.

 0

Returns:

Name Type Description
result dict

Aggregated calibration results: - 'pair_calibrations': list of per-pair calibration dicts - 'D': estimated depth (same for all pairs, depends only on W) - 'Cf_median': median Cf across pairs - 'Cf_values': array of all Cf values - 'kl_median': median kl - 'kl_25': 25th percentile of kl - 'kl_75': 75th percentile of kl - 'kl_values': array of all kl values - 'r_predicted_median': median prediction correlation - 'r_predicted_values': array of all prediction correlations - 'n_pairs_calibrated': number of pairs used - 'W_mean': mean channel width across pairs
Source code in rivabar/prediction.py 
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
def calibrate_segment(results, df, pair_class_filter='good',
                      Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5,
                      kl_percentile=75, use_stable_segments=True,
                      min_time_gap_years=0, pair_indices=None):
    """
    Calibrate the HK model for all qualifying pairs in a segment.

    Parameters
    ----------
    results : dict
        Results dictionary from ``analyze_river_pairs_filtered`` or
        ``analyze_segment_group``.
    df : pandas.DataFrame
        DataFrame from ``create_dataframe_from_results``, used to filter
        pairs by classification.
    pair_class_filter : str or list of str, optional
        Only calibrate pairs with these classifications (default 'good').
    Cf_range : tuple, optional
        Bounds for Cf optimization.
    Omega : float, optional
        Local weight (default -1.0).
    Gamma : float, optional
        Upstream weight (default 2.5).
    kl_percentile : float, optional
        Percentile for amplitude matching (default 75).
    use_stable_segments : bool, optional
        Whether to exclude high-variance segments during calibration
        (default True).
    min_time_gap_years : float, optional
        Minimum time gap (in years) for a pair to be included in
        calibration (default 0). Short-gap pairs tend to have noisier
        migration estimates, which can inflate kl.

    Returns
    -------
    result : dict
        Aggregated calibration results:
        - 'pair_calibrations': list of per-pair calibration dicts
        - 'D': estimated depth (same for all pairs, depends only on W)
        - 'Cf_median': median Cf across pairs
        - 'Cf_values': array of all Cf values
        - 'kl_median': median kl
        - 'kl_25': 25th percentile of kl
        - 'kl_75': 75th percentile of kl
        - 'kl_values': array of all kl values
        - 'r_predicted_median': median prediction correlation
        - 'r_predicted_values': array of all prediction correlations
        - 'n_pairs_calibrated': number of pairs used
        - 'W_mean': mean channel width across pairs
    """
    if pair_indices is not None:
        # Caller supplies an explicit list of pair indices (e.g. the
        # training subset in temporal_cross_validate)
        qualifying = list(pair_indices)
    else:
        if isinstance(pair_class_filter, str):
            pair_class_filter = [pair_class_filter]

        # Find qualifying pairs
        qualifying = df.index[df['pair_class'].isin(pair_class_filter)].tolist()

        # Filter by minimum time gap
        if min_time_gap_years > 0:
            qualifying = [idx for idx in qualifying
                          if results['pair_info'][idx]['time_gap_years']
                          >= min_time_gap_years]

    pair_calibrations = []
    for idx in qualifying:
        try:
            cal = calibrate_pair(results, idx, Cf_range, Omega, Gamma,
                                 kl_percentile, use_stable_segments)
            pair_calibrations.append(cal)
        except Exception as e:
            print(f"  Pair {idx}: calibration failed - {e}")

    if not pair_calibrations:
        return {
            'pair_calibrations': [],
            'n_pairs_calibrated': 0,
        }

    Cf_values = np.array([c['Cf'] for c in pair_calibrations])
    kl_values = np.array([c['kl'] for c in pair_calibrations])
    r_values = np.array([c['r_predicted'] for c in pair_calibrations])
    W_values = np.array([c['W'] for c in pair_calibrations])
    residual_stds = np.array([c['residual_std'] for c in pair_calibrations])

    return {
        'pair_calibrations': pair_calibrations,
        'D': pair_calibrations[0]['D'],  # same for all (depends on mean W)
        'Cf_median': np.median(Cf_values),
        'Cf_values': Cf_values,
        'kl_median': np.median(kl_values),
        'kl_25': np.percentile(kl_values, 25),
        'kl_75': np.percentile(kl_values, 75),
        'kl_values': kl_values,
        'r_predicted_median': np.median(r_values),
        'r_predicted_values': r_values,
        'n_pairs_calibrated': len(pair_calibrations),
        'W_mean': np.mean(W_values),
        'Omega': Omega,
        'Gamma': Gamma,
        'residual_std': np.median(residual_stds),
    }

predict_forward_local(river, local_calibration, segment_calibration=None, path=None, delta_s=50, smoothing_factor=1000000.0, pixel_size=None, prediction_years=None)

Predict migration using spatially-varying kl from local calibration.

Parameters:

Name Type Description Default
river River

Most recent River object.

 required
local_calibration dict

Output of calibrate_pair_local or calibrate_segment_local.

 required
segment_calibration dict or None

If provided (output of calibrate_segment), uses its 'residual_std' for the uncertainty envelope. Otherwise no envelope is computed.

 None
path list of tuples

Sub-path to use.

 None
delta_s float

Resampling interval (default 50 m).

 50
smoothing_factor float

Spline smoothing factor (default 1e6).

 1000000.0
pixel_size float

Pixel size for width conversion.

 None
prediction_years float

Time horizon in years for forward prediction.

 None

Returns:

Name Type Description
prediction dict
  • 'x', 'y': smoothed centerline coordinates
  • 's': along-channel distance
  • 'curvature': curvature values
  • 'width': channel width
  • 'kl_local': interpolated kl profile on the centerline
  • 'predicted_mr': predicted migration rate (local kl)
  • 'predicted_mr_global': predicted MR with global kl If prediction_years is given:
  • 'predicted_x', 'predicted_y': predicted channel position If segment_calibration is given and prediction_years set:
  • 'predicted_x_low/high', 'predicted_y_low/high': envelope
Source code in rivabar/prediction.py 
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
def predict_forward_local(river, local_calibration, segment_calibration=None,
                          path=None, delta_s=50, smoothing_factor=1e6,
                          pixel_size=None, prediction_years=None):
    """
    Predict migration using spatially-varying kl from local calibration.

    Parameters
    ----------
    river : River
        Most recent River object.
    local_calibration : dict
        Output of ``calibrate_pair_local`` or
        ``calibrate_segment_local``.
    segment_calibration : dict or None
        If provided (output of ``calibrate_segment``), uses its
        ``'residual_std'`` for the uncertainty envelope.  Otherwise
        no envelope is computed.
    path : list of tuples, optional
        Sub-path to use.
    delta_s : float, optional
        Resampling interval (default 50 m).
    smoothing_factor : float, optional
        Spline smoothing factor (default 1e6).
    pixel_size : float, optional
        Pixel size for width conversion.
    prediction_years : float, optional
        Time horizon in years for forward prediction.

    Returns
    -------
    prediction : dict
        - 'x', 'y': smoothed centerline coordinates
        - 's': along-channel distance
        - 'curvature': curvature values
        - 'width': channel width
        - 'kl_local': interpolated kl profile on the centerline
        - 'predicted_mr': predicted migration rate (local kl)
        - 'predicted_mr_global': predicted MR with global kl
        If *prediction_years* is given:
        - 'predicted_x', 'predicted_y': predicted channel position
        If *segment_calibration* is given and *prediction_years* set:
        - 'predicted_x_low/high', 'predicted_y_low/high': envelope
    """
    from .utils import get_width_and_curvature

    x, y, s, width, curvature, _ = get_width_and_curvature(
        river, delta_s=delta_s, smoothing_factor=smoothing_factor,
        path=path, pixel_size=pixel_size)

    W = np.mean(width)

    # Extract parameters from local calibration
    # Works with both calibrate_pair_local and calibrate_segment_local output
    if 's_reference' in local_calibration:
        # calibrate_segment_local output
        if 'kl_local_median' in local_calibration:
            kl_s = local_calibration['kl_local_median']
        else:
            # temporally-windowed calibration: use the most recent window
            kl_s = local_calibration['kl_local_series'][-1]['kl_median']
        s_kl = local_calibration['s_reference']
        Cf = (local_calibration['Cf_global'] if 'Cf_global' in local_calibration
              else local_calibration['Cf'])
        kl_global = np.nanmedian(kl_s)
    else:
        # calibrate_pair_local output
        kl_s = local_calibration.get('kl_local_full',
                                     local_calibration['kl_local'])
        s_kl = local_calibration.get('s_full', local_calibration['s'])
        Cf = local_calibration['Cf']
        kl_global = local_calibration['kl_global']

    D = depth_from_width(W)
    Omega = local_calibration.get('Omega', -1.0)
    Gamma = local_calibration.get('Gamma', 2.5)

    # Interpolate kl onto current centerline using DTW alignment
    if 'x_reference' in local_calibration:
        kl_local = _map_to_reference_s(
            local_calibration['x_reference'],
            local_calibration['y_reference'],
            s_kl, kl_s, x, y, s)
    else:
        # calibrate_pair_local doesn't have reference coords — use
        # the pair's own centerline coords if available, else normalised
        kl_local = np.interp(
            (s - s[0]) / (s[-1] - s[0]) if s[-1] > s[0]
            else np.zeros_like(s),
            (s_kl - s_kl[0]) / (s_kl[-1] - s_kl[0]),
            kl_s)
    # Fill NaN gaps with global kl
    nan_mask = np.isnan(kl_local)
    if np.any(nan_mask):
        kl_local[nan_mask] = kl_global

    # Predicted MR with local kl
    R0_local = nominal_migration_rate(curvature, W, kl_local)
    R1_local = predicted_migration_rate(R0_local, s, D, Cf, Omega, Gamma)

    # Global kl for comparison
    R0_global = nominal_migration_rate(curvature, W, kl_global)
    R1_global = predicted_migration_rate(R0_global, s, D, Cf, Omega, Gamma)

    result = {
        'x': x, 'y': y, 's': s,
        'curvature': curvature,
        'width': width,
        'kl_local': kl_local,
        'kl_global': kl_global,
        'predicted_mr': R1_local,
        'predicted_mr_global': R1_global,
        'D': D, 'Cf': Cf,
    }

    if prediction_years is not None:
        px, py = _migrate_forward_local(
            x.copy(), y.copy(), W, D, Cf, kl_local, s,
            Omega, Gamma, prediction_years, delta_s)
        result['predicted_x'] = px
        result['predicted_y'] = py

        # Uncertainty envelope from segment calibration if available
        if segment_calibration is not None:
            residual_std = segment_calibration.get('residual_std', 0.0)
            (result['predicted_x_low'], result['predicted_y_low'],
             result['predicted_x_high'], result['predicted_y_high']) = \
                _envelope_offsets(px, py, residual_std * prediction_years)

    return result

predict_forward(river, calibration, path=None, delta_s=50, smoothing_factor=1000000.0, pixel_size=None, prediction_years=None)

Predict migration rate from the most recent centerline using calibrated HK model parameters.

Parameters:

Name Type Description Default
river River

Most recent River object.

 required
calibration dict

Output of calibrate_segment. Must have at least 'Cf_median', 'kl_median', 'kl_25', 'kl_75', 'D', 'Omega', 'Gamma'.

 required
path list of tuples

Sub-path to use. If None, uses river.main_path.

 None
delta_s float

Resampling interval (default 50 m).

 50
smoothing_factor float

Spline smoothing factor (default 1e6).

 1000000.0
pixel_size float

Pixel size for width conversion.

 None
prediction_years float

Time horizon in years for computing predicted channel position. If None, only migration rate is returned (no position prediction).

 None

Returns:

Name Type Description
prediction dict
  • 'x', 'y': smoothed centerline coordinates
  • 's': along-channel distance
  • 'curvature': curvature values
  • 'width': channel width
  • 'predicted_mr': predicted migration rate (median kl)
  • 'predicted_mr_low': predicted MR minus residual std
  • 'predicted_mr_high': predicted MR plus residual std
  • 'residual_std': median residual std from calibration pairs
  • 'D', 'Cf', 'kl_median', 'kl_25', 'kl_75': parameters used If prediction_years is given, also includes:
  • 'predicted_x', 'predicted_y': predicted channel position
  • 'predicted_x_low', 'predicted_y_low': lower bound position
  • 'predicted_x_high', 'predicted_y_high': upper bound position
Source code in rivabar/prediction.py 
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
def predict_forward(river, calibration, path=None, delta_s=50,
                    smoothing_factor=1e6, pixel_size=None,
                    prediction_years=None):
    """
    Predict migration rate from the most recent centerline using calibrated
    HK model parameters.

    Parameters
    ----------
    river : River
        Most recent River object.
    calibration : dict
        Output of ``calibrate_segment``.  Must have at least
        ``'Cf_median'``, ``'kl_median'``, ``'kl_25'``, ``'kl_75'``,
        ``'D'``, ``'Omega'``, ``'Gamma'``.
    path : list of tuples, optional
        Sub-path to use.  If *None*, uses ``river.main_path``.
    delta_s : float, optional
        Resampling interval (default 50 m).
    smoothing_factor : float, optional
        Spline smoothing factor (default 1e6).
    pixel_size : float, optional
        Pixel size for width conversion.
    prediction_years : float, optional
        Time horizon in years for computing predicted channel position.
        If *None*, only migration rate is returned (no position prediction).

    Returns
    -------
    prediction : dict
        - 'x', 'y': smoothed centerline coordinates
        - 's': along-channel distance
        - 'curvature': curvature values
        - 'width': channel width
        - 'predicted_mr': predicted migration rate (median kl)
        - 'predicted_mr_low': predicted MR minus residual std
        - 'predicted_mr_high': predicted MR plus residual std
        - 'residual_std': median residual std from calibration pairs
        - 'D', 'Cf', 'kl_median', 'kl_25', 'kl_75': parameters used
        If *prediction_years* is given, also includes:
        - 'predicted_x', 'predicted_y': predicted channel position
        - 'predicted_x_low', 'predicted_y_low': lower bound position
        - 'predicted_x_high', 'predicted_y_high': upper bound position
    """
    from .utils import get_width_and_curvature

    x, y, s, width, curvature, _ = get_width_and_curvature(
        river, delta_s=delta_s, smoothing_factor=smoothing_factor,
        path=path, pixel_size=pixel_size)

    W = np.mean(width)
    D = calibration['D']
    Cf = calibration['Cf_median']
    Omega = calibration.get('Omega', -1.0)
    Gamma = calibration.get('Gamma', 2.5)

    def _predict_with_kl(kl):
        R0 = nominal_migration_rate(curvature, W, kl)
        return predicted_migration_rate(R0, s, D, Cf, Omega, Gamma)

    mr_median = _predict_with_kl(calibration['kl_median'])

    # Uncertainty envelope based on residual std from calibration pairs
    residual_std = calibration.get('residual_std', 0.0)
    mr_low = mr_median - residual_std
    mr_high = mr_median + residual_std

    result = {
        'x': x,
        'y': y,
        's': s,
        'curvature': curvature,
        'width': width,
        'predicted_mr': mr_median,
        'predicted_mr_low': mr_low,
        'predicted_mr_high': mr_high,
        'residual_std': residual_std,
        'D': D,
        'Cf': Cf,
        'kl_median': calibration['kl_median'],
        'kl_25': calibration['kl_25'],
        'kl_75': calibration['kl_75'],
    }

    if prediction_years is not None:
        # Run forward model with median parameters
        px, py = _migrate_forward(
            x.copy(), y.copy(), W, D, Cf, calibration['kl_median'],
            Omega, Gamma, prediction_years, delta_s)
        result['predicted_x'] = px
        result['predicted_y'] = py

        # Build envelope by offsetting the predicted centerline by
        # +/- residual_std * prediction_years along the local normal
        (result['predicted_x_low'], result['predicted_y_low'],
         result['predicted_x_high'], result['predicted_y_high']) = \
            _envelope_offsets(px, py, residual_std * prediction_years)

    return result

temporal_cross_validate(results, df, pair_class_filter='good', n_holdout=1, holdout_date=None, Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5, kl_percentile=75, use_stable_segments=True, use_forward_model=False, delta_s=50, min_time_gap_years=0, max_time_gap_years=None, min_prediction_years=0, use_local_kl=False, window_length=None)

Calibrate on older pairs and evaluate on held-out recent pairs.

Test pairs are those where date1 is the last pre-cutoff scene and date2 is post-cutoff, mimicking the forward-prediction scenario: the model sees only the training-era curvature and must predict migration into unseen time.

Parameters:

Name Type Description Default
results dict

Results dictionary from analyze_river_pairs_filtered or analyze_segment_group.

 required
df DataFrame

DataFrame from create_dataframe_from_results, used to filter pairs by classification.

 required
pair_class_filter str or list of str

Only use pairs with these classifications (default 'good').

 'good'
n_holdout int

Number of the most recent scenes to hold out (default 1). Ignored if holdout_date is given.

 1
holdout_date datetime - like

Explicit cutoff date. Pairs with both dates before this are training; pairs with date1 before and date2 on or after are test. If None, the cutoff is inferred from n_holdout.

 None
Cf_range tuple

Bounds for Cf optimization (default (0.001, 0.02)).

 (0.001, 0.02)
Omega float

Local weight (default -1.0).

 -1.0
Gamma float

Upstream weight (default 2.5).

 2.5
kl_percentile float

Percentile for amplitude matching (default 75).

 75
use_stable_segments bool

Whether to exclude high-variance segments during calibration (default True).

 True
use_forward_model bool

If True, test evaluation uses the iterative forward model (_migrate_forward) to predict the time-2 centerline from time-1 geometry, then compares against the observed time-2 centerline. This accounts for curvature evolution during migration. Requires results['centerline_coords']. If False (default), uses the static single-step HK prediction.

 False
delta_s float

Node spacing (m) for the forward model resampling (default 50). Only used when use_forward_model is True.

 50
min_time_gap_years float

Minimum time gap (in years) for a pair to be included in the training set (default 0). Short-gap pairs tend to have noisier migration estimates, which can inflate kl. This filter only applies to training pairs; test pairs are always evaluated regardless of their time gap.

 0
max_time_gap_years float or None

Maximum time gap (in years) for a pair to be included in the training set (default None, no upper limit). Long-gap pairs may have unreliable migration estimates due to cutoffs and cumulative geometry changes. Only applies to training pairs; test pairs are always evaluated to allow assessing the predictability horizon.

 None
min_prediction_years float

Minimum time (in years) between the cutoff date and a test pair's date2 (default 0). Excludes test pairs whose date2 is too close to the cutoff, where minimal channel change makes prediction trivially easy.

 0
use_local_kl bool

If True, calibrate spatially-varying kl on training pairs and use the aggregated kl(s) profile for test evaluation (both static and forward model). Default False.

 False
window_length float or None

Spatial window for local kl estimation (metres). Auto-estimated if None. Only used when use_local_kl is True.

 None

Returns:

Name Type Description
result dict
  • 'calibration': training-set calibration dict (same format as calibrate_segment output)
  • 'test_pairs': list of dicts, one per test pair, each with:
    • 'pair_idx': index in results
    • 'date1', 'date2': scene dates
    • 'time_gap_years': time between scenes
    • 'r_predicted': Pearson r between HK-predicted and observed MR
    • 'rmse': root mean squared error (m/yr)
    • 'predicted_mr': predicted migration rate array
    • 'observed_mr': observed migration rate array
    • 's': along-channel distance array If use_forward_model is True, also includes:
    • 'predicted_x', 'predicted_y': forward-predicted centerline
    • 'observed_x2', 'observed_y2': observed time-2 centerline
    • 'position_error_mean': mean Euclidean distance between predicted and observed time-2 centerlines (m)
    • 'position_error_median': median distance (m)
    • 'r_forward': correlation between forward-model migration distances and observed migration distances
  • 'cutoff_date': the date used to split train/test
  • 'n_train': number of training pairs
  • 'n_test': number of test pairs
  • 'r_test_median': median correlation across test pairs
  • 'rmse_test_median': median RMSE across test pairs
  • 'train_pair_indices': list of training pair indices
  • 'test_pair_indices': list of test pair indices
  • 'use_forward_model': whether the forward model was used If use_local_kl is True, also includes:
  • 'local_calibration': dict with 'kl_median', 'kl_25', 'kl_75', 's_reference' arrays from the training-set local calibration
Source code in rivabar/prediction.py 
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
def temporal_cross_validate(results, df, pair_class_filter='good',
                            n_holdout=1, holdout_date=None,
                            Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5,
                            kl_percentile=75, use_stable_segments=True,
                            use_forward_model=False, delta_s=50,
                            min_time_gap_years=0,
                            max_time_gap_years=None,
                            min_prediction_years=0,
                            use_local_kl=False, window_length=None):
    """
    Calibrate on older pairs and evaluate on held-out recent pairs.

    Test pairs are those where date1 is the last pre-cutoff scene and date2
    is post-cutoff, mimicking the forward-prediction scenario: the model
    sees only the training-era curvature and must predict migration into
    unseen time.

    Parameters
    ----------
    results : dict
        Results dictionary from ``analyze_river_pairs_filtered`` or
        ``analyze_segment_group``.
    df : pandas.DataFrame
        DataFrame from ``create_dataframe_from_results``, used to filter
        pairs by classification.
    pair_class_filter : str or list of str, optional
        Only use pairs with these classifications (default 'good').
    n_holdout : int, optional
        Number of the most recent *scenes* to hold out (default 1).
        Ignored if *holdout_date* is given.
    holdout_date : datetime-like, optional
        Explicit cutoff date. Pairs with both dates before this are
        training; pairs with date1 before and date2 on or after are test.
        If *None*, the cutoff is inferred from *n_holdout*.
    Cf_range : tuple, optional
        Bounds for Cf optimization (default (0.001, 0.02)).
    Omega : float, optional
        Local weight (default -1.0).
    Gamma : float, optional
        Upstream weight (default 2.5).
    kl_percentile : float, optional
        Percentile for amplitude matching (default 75).
    use_stable_segments : bool, optional
        Whether to exclude high-variance segments during calibration
        (default True).
    use_forward_model : bool, optional
        If True, test evaluation uses the iterative forward model
        (``_migrate_forward``) to predict the time-2 centerline from
        time-1 geometry, then compares against the observed time-2
        centerline. This accounts for curvature evolution during
        migration. Requires ``results['centerline_coords']``.
        If False (default), uses the static single-step HK prediction.
    delta_s : float, optional
        Node spacing (m) for the forward model resampling (default 50).
        Only used when *use_forward_model* is True.
    min_time_gap_years : float, optional
        Minimum time gap (in years) for a pair to be included in the
        training set (default 0). Short-gap pairs tend to have noisier
        migration estimates, which can inflate kl. This filter only
        applies to training pairs; test pairs are always evaluated
        regardless of their time gap.
    max_time_gap_years : float or None, optional
        Maximum time gap (in years) for a pair to be included in the
        training set (default *None*, no upper limit). Long-gap pairs
        may have unreliable migration estimates due to cutoffs and
        cumulative geometry changes. Only applies to training pairs;
        test pairs are always evaluated to allow assessing the
        predictability horizon.
    min_prediction_years : float, optional
        Minimum time (in years) between the cutoff date and a test
        pair's date2 (default 0).  Excludes test pairs whose date2 is
        too close to the cutoff, where minimal channel change makes
        prediction trivially easy.
    use_local_kl : bool, optional
        If True, calibrate spatially-varying kl on training pairs and
        use the aggregated kl(s) profile for test evaluation (both
        static and forward model).  Default False.
    window_length : float or None, optional
        Spatial window for local kl estimation (metres).  Auto-estimated
        if *None*.  Only used when *use_local_kl* is True.

    Returns
    -------
    result : dict
        - 'calibration': training-set calibration dict (same format as
          ``calibrate_segment`` output)
        - 'test_pairs': list of dicts, one per test pair, each with:
            - 'pair_idx': index in results
            - 'date1', 'date2': scene dates
            - 'time_gap_years': time between scenes
            - 'r_predicted': Pearson r between HK-predicted and observed MR
            - 'rmse': root mean squared error (m/yr)
            - 'predicted_mr': predicted migration rate array
            - 'observed_mr': observed migration rate array
            - 's': along-channel distance array
            If *use_forward_model* is True, also includes:
            - 'predicted_x', 'predicted_y': forward-predicted centerline
            - 'observed_x2', 'observed_y2': observed time-2 centerline
            - 'position_error_mean': mean Euclidean distance between
              predicted and observed time-2 centerlines (m)
            - 'position_error_median': median distance (m)
            - 'r_forward': correlation between forward-model migration
              distances and observed migration distances
        - 'cutoff_date': the date used to split train/test
        - 'n_train': number of training pairs
        - 'n_test': number of test pairs
        - 'r_test_median': median correlation across test pairs
        - 'rmse_test_median': median RMSE across test pairs
        - 'train_pair_indices': list of training pair indices
        - 'test_pair_indices': list of test pair indices
        - 'use_forward_model': whether the forward model was used
        If *use_local_kl* is True, also includes:
        - 'local_calibration': dict with 'kl_median', 'kl_25', 'kl_75',
          's_reference' arrays from the training-set local calibration
    """

    if isinstance(pair_class_filter, str):
        pair_class_filter = [pair_class_filter]

    # Find qualifying pairs
    qualifying = df.index[df['pair_class'].isin(pair_class_filter)].tolist()
    if len(qualifying) < 2:
        raise ValueError(
            f"Need at least 2 qualifying pairs for cross-validation, "
            f"got {len(qualifying)}")

    # Collect all unique scene dates from qualifying pairs
    all_dates = set()
    for idx in qualifying:
        pi = results['pair_info'][idx]
        all_dates.add(pi['date1'])
        all_dates.add(pi['date2'])
    sorted_dates = sorted(all_dates)

    # Determine cutoff date
    if holdout_date is not None:
        cutoff = holdout_date
    else:
        if n_holdout >= len(sorted_dates):
            raise ValueError(
                f"n_holdout={n_holdout} but only {len(sorted_dates)} "
                f"unique scene dates")
        cutoff = sorted_dates[-n_holdout]

    # Split into training and test pairs
    train_indices = []
    test_indices = []
    for idx in qualifying:
        pi = results['pair_info'][idx]
        if pi['date1'] < cutoff and pi['date2'] < cutoff:
            train_indices.append(idx)
        elif pi['date1'] < cutoff and pi['date2'] >= cutoff:
            test_indices.append(idx)
        # Pairs where both dates are post-cutoff are excluded (no
        # training-era curvature available)

    if len(train_indices) == 0:
        raise ValueError("No training pairs before cutoff date "
                         f"{cutoff}. Try an earlier cutoff or fewer "
                         "holdout scenes.")
    if len(test_indices) == 0:
        raise ValueError("No test pairs crossing the cutoff date "
                         f"{cutoff}. Need pairs where date1 < cutoff "
                         "and date2 >= cutoff.")

    if min_prediction_years > 0:
        cutoff_ts = pd.Timestamp(cutoff)
        min_pred_delta = pd.Timedelta(days=min_prediction_years * 365.25)
        test_indices = [
            idx for idx in test_indices
            if pd.Timestamp(results['pair_info'][idx]['date2'])
            >= cutoff_ts + min_pred_delta]
        if len(test_indices) == 0:
            raise ValueError(
                f"No test pairs with date2 >= {min_prediction_years} "
                f"years after cutoff {cutoff}.")

    # --- Training: calibrate on older pairs ---
    # Filter training pairs by minimum time gap
    if min_time_gap_years > 0:
        train_indices = [idx for idx in train_indices
                         if results['pair_info'][idx]['time_gap_years']
                         >= min_time_gap_years]
        if len(train_indices) == 0:
            raise ValueError(
                f"No training pairs with time gap >= "
                f"{min_time_gap_years} years. Try a smaller threshold.")

    if max_time_gap_years is not None:
        train_indices = [idx for idx in train_indices
                         if results['pair_info'][idx]['time_gap_years']
                         <= max_time_gap_years]
        if len(train_indices) == 0:
            raise ValueError(
                f"No training pairs with time gap <= "
                f"{max_time_gap_years} years.")

    calibration = calibrate_segment(results, df, Cf_range=Cf_range,
                                    Omega=Omega, Gamma=Gamma,
                                    kl_percentile=kl_percentile,
                                    use_stable_segments=use_stable_segments,
                                    pair_indices=train_indices)
    if calibration['n_pairs_calibrated'] == 0:
        raise ValueError("All training pair calibrations failed")

    # --- Local kl training (optional) ---
    local_cal_data = None
    if use_local_kl:
        # Calibrate local kl on each training pair, then aggregate
        local_profiles = []
        local_s_grids = []
        local_coords = []
        for idx in tqdm(train_indices, desc='Calibrating local kl'):
            try:
                lcal = calibrate_pair_local(
                    results, idx,
                    use_stable_segments=use_stable_segments,
                    Cf=calibration['Cf_median'],
                    Omega=Omega, Gamma=Gamma,
                    kl_percentile=kl_percentile,
                    window_length=window_length)
                kl_prof = lcal.get('kl_local_full', lcal['kl_local'])
                s_prof = lcal.get('s_full', lcal['s'])
                coords = results['centerline_coords'][idx]
                local_profiles.append(kl_prof)
                local_s_grids.append(s_prof)
                local_coords.append(coords)
                if window_length is None:
                    window_length = lcal['window_length']
            except Exception as e:
                print(f"  Training pair {idx}: local kl failed - {e}")

        if local_profiles:
            # Use the first training pair's centerline as reference
            x_ref = local_coords[0]['x1']
            y_ref = local_coords[0]['y1']
            from .utils import compute_s_distance
            s_ref = compute_s_distance(x_ref, y_ref)

            remapped = []
            for kl_prof, s_prof, coords in zip(
                    local_profiles, local_s_grids, local_coords):
                kl_on_ref = _map_to_reference_s(
                    coords['x1'], coords['y1'], s_prof, kl_prof,
                    x_ref, y_ref, s_ref)
                if not np.all(np.isnan(kl_on_ref)):
                    remapped.append(kl_on_ref)

            if remapped:
                stack = np.array(remapped)
                local_cal_data = {
                    'kl_median': np.nanmedian(stack, axis=0),
                    'kl_25': np.nanpercentile(stack, 25, axis=0),
                    'kl_75': np.nanpercentile(stack, 75, axis=0),
                    's_reference': s_ref,
                    'x_reference': x_ref,
                    'y_reference': y_ref,
                    'Cf': calibration['Cf_median'],
                    'D': calibration['D'],
                    'window_length': window_length,
                    'n_profiles': len(remapped),
                }

    # --- Test: evaluate on held-out pairs ---
    Cf_trained = calibration['Cf_median']
    kl_trained = calibration['kl_median']
    D_trained = calibration['D']

    if use_forward_model:
        from .utils import correlate_curves, compute_migration_distances

    test_results = []
    for idx in tqdm(test_indices, desc='Evaluating test pairs'):
        pi = results['pair_info'][idx]
        curvature = results['curvatures'][idx]
        distances = results['migration_distances'][idx]
        s = results['along_channel_distances'][idx]
        time_gap_years = pi['time_gap_years']
        W = np.mean(pi['width1'])

        observed_mr = distances / time_gap_years

        # --- Determine kl for this test pair ---
        if use_local_kl and local_cal_data is not None:
            # Map the aggregated kl(s) profile onto this test pair's
            # centerline using DTW geometric alignment
            coords = results['centerline_coords'][idx]
            kl_for_test = _map_to_reference_s(
                local_cal_data['x_reference'],
                local_cal_data['y_reference'],
                local_cal_data['s_reference'],
                local_cal_data['kl_median'],
                coords['x1'], coords['y1'], s)
            # Fill any NaN gaps (e.g. at endpoints) with global kl
            nan_mask = np.isnan(kl_for_test)
            if np.any(nan_mask):
                kl_for_test[nan_mask] = kl_trained
        else:
            kl_for_test = kl_trained

        # --- Static prediction (always computed) ---
        R0 = nominal_migration_rate(curvature, W, kl_for_test)
        R1 = predicted_migration_rate(R0, s, D_trained, Cf_trained,
                                      Omega, Gamma)

        # Evaluate on stable segments only if available
        mask = _stable_segment_mask(pi, len(curvature)) if use_stable_segments else None
        if mask is not None:
            R1_eval = R1[mask]
            obs_eval = observed_mr[mask]
        else:
            R1_eval = R1
            obs_eval = observed_mr

        if len(R1_eval) > 2:
            r = np.corrcoef(R1_eval, obs_eval)[0, 1]
            rmse = np.sqrt(np.mean((R1_eval - obs_eval) ** 2))
        else:
            r = np.nan
            rmse = np.nan

        test_entry = {
            'pair_idx': idx,
            'date1': pi['date1'],
            'date2': pi['date2'],
            'time_gap_years': time_gap_years,
            'r_predicted': r,
            'rmse': rmse,
            'predicted_mr': R1,
            'observed_mr': observed_mr,
            's': s,
        }

        # --- Iterative forward model ---
        if use_forward_model:
            try:
                coords = results['centerline_coords'][idx]
                x1, y1 = coords['x1'].copy(), coords['y1'].copy()
                x2_obs, y2_obs = coords['x2'].copy(), coords['y2'].copy()

                # Forward-integrate from time-1 geometry
                if use_local_kl and local_cal_data is not None:
                    from .utils import compute_s_distance as _cs
                    s1_fwd = _cs(x1, y1)
                    kl_fwd = _map_to_reference_s(
                        local_cal_data['x_reference'],
                        local_cal_data['y_reference'],
                        local_cal_data['s_reference'],
                        local_cal_data['kl_median'],
                        x1, y1, s1_fwd)
                    nan_mask = np.isnan(kl_fwd)
                    if np.any(nan_mask):
                        kl_fwd[nan_mask] = kl_trained
                    px, py = _migrate_forward_local(
                        x1.copy(), y1.copy(), W, D_trained, Cf_trained,
                        kl_fwd, s1_fwd, Omega, Gamma,
                        time_gap_years, delta_s)
                else:
                    px, py = _migrate_forward(
                        x1.copy(), y1.copy(), W, D_trained, Cf_trained,
                        kl_trained, Omega, Gamma, time_gap_years, delta_s)

                # DTW-align predicted time-2 against observed time-2
                p_fwd, q_fwd, cost_fwd = correlate_curves(
                    px, x2_obs, py, y2_obs)
                del cost_fwd

                # Position error: Euclidean distance between DTW-matched
                # points on predicted and observed time-2 centerlines
                pos_errors = np.sqrt(
                    (px[p_fwd] - x2_obs[q_fwd]) ** 2 +
                    (py[p_fwd] - y2_obs[q_fwd]) ** 2)

                # DTW cost: time-1 vs forward-predicted time-2
                # Compare against the observed DTW cost to gauge
                # whether the predicted shape is as realistic as the
                # observed one
                p_t1_pred, q_t1_pred, dtw_cost_matrix = \
                    correlate_curves(x1, px, y1, py)
                dtw_cost_predicted = float(dtw_cost_matrix[-1, -1])
                del dtw_cost_matrix
                dtw_cost_observed = pi.get('dtw_cost', np.nan)

                # Migration distances: forward-predicted vs observed,
                # both measured from time-1 centerline
                fwd_distances, fwd_valid = compute_migration_distances(
                    x1, y1, px, py, p_t1_pred, q_t1_pred)

                # Correlation of forward-model migration vs observed
                # migration (both from time-1), restricted to stable
                # segments when requested
                both_valid = fwd_valid & ~np.isnan(distances)
                stable_mask = _stable_segment_mask(pi, len(curvature)) if use_stable_segments else None
                if stable_mask is not None:
                    both_valid = both_valid & stable_mask
                if np.sum(both_valid) > 2:
                    fwd_mr = fwd_distances[both_valid] / time_gap_years
                    obs_mr_valid = observed_mr[both_valid]
                    r_forward = np.corrcoef(fwd_mr, obs_mr_valid)[0, 1]
                    rmse_forward = np.sqrt(
                        np.mean((fwd_mr - obs_mr_valid) ** 2))
                else:
                    r_forward = np.nan
                    rmse_forward = np.nan

                test_entry.update({
                    'predicted_x': px,
                    'predicted_y': py,
                    'observed_x2': x2_obs,
                    'observed_y2': y2_obs,
                    'position_error_mean': np.mean(pos_errors),
                    'position_error_median': np.median(pos_errors),
                    'r_forward': r_forward,
                    'rmse_forward': rmse_forward,
                    'observed_distances': distances,
                    'forward_distances': fwd_distances,
                    'dtw_cost_predicted': dtw_cost_predicted,
                    'dtw_cost_observed': dtw_cost_observed,
                })
            except Exception as e:
                print(f"  Test pair {idx}: forward model failed - {e}")
                test_entry.update({
                    'r_forward': np.nan,
                    'rmse_forward': np.nan,
                    'position_error_mean': np.nan,
                    'position_error_median': np.nan,
                })

        test_results.append(test_entry)

    r_test = [t['r_predicted'] for t in test_results
              if not np.isnan(t['r_predicted'])]
    rmse_test = [t['rmse'] for t in test_results
                 if not np.isnan(t['rmse'])]

    output = {
        'calibration': calibration,
        'test_pairs': test_results,
        'cutoff_date': cutoff,
        'n_train': len(train_indices),
        'n_test': len(test_indices),
        'r_test_median': np.median(r_test) if r_test else np.nan,
        'rmse_test_median': np.median(rmse_test) if rmse_test else np.nan,
        'train_pair_indices': train_indices,
        'test_pair_indices': test_indices,
        'use_forward_model': use_forward_model,
    }

    if use_forward_model:
        r_fwd = [t['r_forward'] for t in test_results
                 if not np.isnan(t.get('r_forward', np.nan))]
        rmse_fwd = [t['rmse_forward'] for t in test_results
                    if not np.isnan(t.get('rmse_forward', np.nan))]
        pos_err = [t['position_error_mean'] for t in test_results
                   if not np.isnan(t.get('position_error_mean', np.nan))]
        output['r_forward_median'] = np.median(r_fwd) if r_fwd else np.nan
        output['rmse_forward_median'] = (
            np.median(rmse_fwd) if rmse_fwd else np.nan)
        output['position_error_mean_median'] = (
            np.median(pos_err) if pos_err else np.nan)

    if use_local_kl and local_cal_data is not None:
        output['local_calibration'] = local_cal_data

    return output

detect_cutoff_risk(river, path=None, neck_width_threshold=2.0, delta_s=50, smoothing_factor=1000000.0, pixel_size=None, search_distance_wavelengths=0.5)

Identify bends approaching meander cutoff by measuring neck width.

For each bend apex (curvature extremum), measures the minimum distance from the apex to other parts of the centerline beyond the adjacent inflection points. When this distance drops below neck_width_threshold * channel_width, the bend is flagged.

Parameters:

Name Type Description Default
river River

River object.

 required
path list of tuples

Sub-path to use.

 None
neck_width_threshold float

Cutoff warning when neck_width / channel_width < this value (default 2.0).

 2.0
delta_s float

Resampling interval (default 50 m).

 50
smoothing_factor float

Spline smoothing factor (default 1e6).

 1000000.0
pixel_size float

Pixel size for width conversion.

 None
search_distance_wavelengths float

Minimum along-channel distance (in wavelengths) to skip when searching for the nearest point on the centerline (default 0.5).

 0.5

Returns:

Name Type Description
risks list of dict

One entry per bend, each with: - 'bend_index': index of the curvature extremum - 'apex_x', 'apex_y': coordinates of the bend apex - 'along_channel_distance': s-distance of the apex - 'neck_width_m': minimum distance to another part of the channel - 'channel_width_m': local channel width - 'neck_ratio': neck_width / channel_width - 'risk_level': 'high', 'moderate', or 'low'
Source code in rivabar/prediction.py 
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
def detect_cutoff_risk(river, path=None, neck_width_threshold=2.0,
                       delta_s=50, smoothing_factor=1e6, pixel_size=None,
                       search_distance_wavelengths=0.5):
    """
    Identify bends approaching meander cutoff by measuring neck width.

    For each bend apex (curvature extremum), measures the minimum distance
    from the apex to other parts of the centerline beyond the adjacent
    inflection points. When this distance drops below
    ``neck_width_threshold * channel_width``, the bend is flagged.

    Parameters
    ----------
    river : River
        River object.
    path : list of tuples, optional
        Sub-path to use.
    neck_width_threshold : float, optional
        Cutoff warning when neck_width / channel_width < this value
        (default 2.0).
    delta_s : float, optional
        Resampling interval (default 50 m).
    smoothing_factor : float, optional
        Spline smoothing factor (default 1e6).
    pixel_size : float, optional
        Pixel size for width conversion.
    search_distance_wavelengths : float, optional
        Minimum along-channel distance (in wavelengths) to skip when
        searching for the nearest point on the centerline (default 0.5).

    Returns
    -------
    risks : list of dict
        One entry per bend, each with:
        - 'bend_index': index of the curvature extremum
        - 'apex_x', 'apex_y': coordinates of the bend apex
        - 'along_channel_distance': s-distance of the apex
        - 'neck_width_m': minimum distance to another part of the channel
        - 'channel_width_m': local channel width
        - 'neck_ratio': neck_width / channel_width
        - 'risk_level': 'high', 'moderate', or 'low'
    """
    from .utils import get_width_and_curvature

    x, y, s, width, curvature, _ = get_width_and_curvature(
        river, delta_s=delta_s, smoothing_factor=smoothing_factor,
        path=path, pixel_size=pixel_size)

    W_mean = np.mean(width)

    # Estimate wavelength as 2 * pi * W_mean * some factor
    # Use a simpler approach: skip distance = search_distance_wavelengths
    # times mean meander wavelength. Approximate wavelength ~ 12 * W_mean.
    skip_distance = search_distance_wavelengths * 12.0 * W_mean

    # Find curvature extrema (bend apexes)
    from scipy.signal import argrelextrema
    maxima = argrelextrema(np.abs(curvature), np.greater, order=5)[0]

    risks = []
    for apex_idx in maxima:
        apex_s = s[apex_idx]
        local_width = width[apex_idx]

        # Find points far enough away along the channel
        s_diff = np.abs(s - apex_s)
        far_mask = s_diff > skip_distance
        far_indices = np.where(far_mask)[0]

        if len(far_indices) == 0:
            continue

        # Find minimum distance from apex to far-away points
        dists = np.sqrt((x[far_indices] - x[apex_idx])**2 +
                        (y[far_indices] - y[apex_idx])**2)
        neck_width = np.min(dists)
        neck_ratio = neck_width / local_width if local_width > 0 else np.inf

        if neck_ratio > 5.0:
            risk_level = 'low'
        elif neck_ratio > neck_width_threshold:
            risk_level = 'moderate'
        else:
            risk_level = 'high'

        risks.append({
            'bend_index': int(apex_idx),
            'apex_x': x[apex_idx],
            'apex_y': y[apex_idx],
            'along_channel_distance': apex_s,
            'neck_width_m': neck_width,
            'channel_width_m': local_width,
            'neck_ratio': neck_ratio,
            'risk_level': risk_level,
        })

    return risks

track_parameter_stability(pair_calibrations)

Analyze how calibrated parameters vary across time intervals.

Parameters:

Name Type Description Default
pair_calibrations list of dict

Per-pair calibration results from calibrate_segment.

 required

Returns:

Name Type Description
stability dict
  • 'temporal_df': DataFrame with mid_date, D, Cf, kl, r_predicted, time_gap_years
  • 'kl_trend': dict with 'rho' (Spearman), 'p_value'
  • 'Cf_trend': dict with 'rho', 'p_value'
Source code in rivabar/prediction.py 
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
def track_parameter_stability(pair_calibrations):
    """
    Analyze how calibrated parameters vary across time intervals.

    Parameters
    ----------
    pair_calibrations : list of dict
        Per-pair calibration results from ``calibrate_segment``.

    Returns
    -------
    stability : dict
        - 'temporal_df': DataFrame with mid_date, D, Cf, kl, r_predicted,
          time_gap_years
        - 'kl_trend': dict with 'rho' (Spearman), 'p_value'
        - 'Cf_trend': dict with 'rho', 'p_value'
    """

    records = []
    for cal in pair_calibrations:
        mid_date = cal['date1'] + (cal['date2'] - cal['date1']) / 2
        records.append({
            'mid_date': mid_date,
            'date1': cal['date1'],
            'date2': cal['date2'],
            'D': cal['D'],
            'Cf': cal['Cf'],
            'kl': cal['kl'],
            'r_predicted': cal['r_predicted'],
            'time_gap_years': cal['time_gap_years'],
        })

    temporal_df = pd.DataFrame(records).sort_values('mid_date')

    result = {'temporal_df': temporal_df}

    if len(temporal_df) >= 3:
        mid_ordinal = temporal_df['mid_date'].apply(
            lambda d: d.toordinal() if hasattr(d, 'toordinal') else d)

        for param in ['kl', 'Cf']:
            rho, p = stats.spearmanr(mid_ordinal, temporal_df[param])
            result[f'{param}_trend'] = {'rho': rho, 'p_value': p}
    else:
        for param in ['kl', 'Cf']:
            result[f'{param}_trend'] = {'rho': np.nan, 'p_value': np.nan}

    return result

calibrate_local_kl(curvature, observed_mr, s, W, Cf=None, Cf_range=(0.001, 0.02), Omega=-1.0, Gamma=2.5, kl_percentile=75, window_length=None, min_window_points=20, mask=None)

Estimate spatially-varying kl along the channel using a moving window.

Uses a global Cf (optimised once over the full river) and slides a spatial window to estimate local kl by amplitude-matching the HK predicted migration rate to the observed migration rate.

Parameters:

Name Type Description Default
curvature array_like

Curvature along the channel (1/m).

 required
observed_mr array_like

Observed migration rate (m/yr).

 required
s array_like

Along-channel distance (m).

 required
W float

Mean channel width (m).

 required
Cf float or None

Friction factor. If None, optimised globally via calibrate_Cf.

 None
Cf_range tuple

Bounds for Cf optimisation (default (0.001, 0.02)).

 (0.001, 0.02)
Omega float

Local weight (default -1.0).

 -1.0
Gamma float

Upstream weight (default 2.5).

 2.5
kl_percentile float

Percentile for amplitude matching (default 75).

 75
window_length float or None

Spatial window length in metres. If None, auto-estimated from curvature zero-crossings (~2 meander wavelengths).

 None
min_window_points int

Minimum number of points in a window for a valid kl estimate (default 20).

 20
mask array_like of bool

Points to use for the calibration statistics (Cf correlation and window percentiles), e.g. a stable-segment mask. The prediction is always computed on the full contiguous arrays, since the convolution assumes uniform spacing.

 None

Returns:

Name Type Description
result dict
  • 'kl_local': array of kl values (same length as curvature)
  • 'kl_window_centers': along-channel positions of window centres
  • 'kl_window_values': kl at each window centre (before interp)
  • 'Cf': global Cf used
  • 'D': depth estimate
  • 'W': mean channel width
  • 'R1_local': predicted migration rate using kl_local
  • 'R1_global': predicted migration rate using global kl
  • 'kl_global': global kl for comparison
  • 'window_length': window length used (metres)
  • 'observed_mr': observed migration rate (pass-through)
  • 's': along-channel distance (pass-through)
Source code in rivabar/prediction.py 
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
def calibrate_local_kl(curvature, observed_mr, s, W,
                       Cf=None, Cf_range=(0.001, 0.02),
                       Omega=-1.0, Gamma=2.5,
                       kl_percentile=75,
                       window_length=None,
                       min_window_points=20,
                       mask=None):
    """
    Estimate spatially-varying kl along the channel using a moving window.

    Uses a global Cf (optimised once over the full river) and slides a
    spatial window to estimate local kl by amplitude-matching the HK
    predicted migration rate to the observed migration rate.

    Parameters
    ----------
    curvature : array_like
        Curvature along the channel (1/m).
    observed_mr : array_like
        Observed migration rate (m/yr).
    s : array_like
        Along-channel distance (m).
    W : float
        Mean channel width (m).
    Cf : float or None
        Friction factor.  If *None*, optimised globally via
        ``calibrate_Cf``.
    Cf_range : tuple, optional
        Bounds for Cf optimisation (default (0.001, 0.02)).
    Omega : float, optional
        Local weight (default -1.0).
    Gamma : float, optional
        Upstream weight (default 2.5).
    kl_percentile : float, optional
        Percentile for amplitude matching (default 75).
    window_length : float or None
        Spatial window length in metres.  If *None*, auto-estimated from
        curvature zero-crossings (~2 meander wavelengths).
    min_window_points : int, optional
        Minimum number of points in a window for a valid kl estimate
        (default 20).
    mask : array_like of bool, optional
        Points to use for the calibration statistics (Cf correlation and
        window percentiles), e.g. a stable-segment mask. The prediction is
        always computed on the full contiguous arrays, since the convolution
        assumes uniform spacing.

    Returns
    -------
    result : dict
        - 'kl_local': array of kl values (same length as curvature)
        - 'kl_window_centers': along-channel positions of window centres
        - 'kl_window_values': kl at each window centre (before interp)
        - 'Cf': global Cf used
        - 'D': depth estimate
        - 'W': mean channel width
        - 'R1_local': predicted migration rate using kl_local
        - 'R1_global': predicted migration rate using global kl
        - 'kl_global': global kl for comparison
        - 'window_length': window length used (metres)
        - 'observed_mr': observed migration rate (pass-through)
        - 's': along-channel distance (pass-through)
    """
    curvature = np.asarray(curvature, dtype=float)
    observed_mr = np.asarray(observed_mr, dtype=float)
    s = np.asarray(s, dtype=float)
    if mask is None:
        mask = np.ones(len(curvature), dtype=bool)

    # --- Step 1: depth ---
    D = depth_from_width(W)

    # --- Step 2: global Cf ---
    # Use kl=1 for initial R0 so that kl cancels out of the Cf optimisation
    R0_unit = nominal_migration_rate(curvature, W, 1.0)
    if Cf is None:
        Cf, _ = calibrate_Cf(R0_unit, observed_mr, s, D, Cf_range,
                             Omega, Gamma, mask=mask)

    # --- Step 3: global kl for comparison ---
    R1_unit = predicted_migration_rate(R0_unit, s, D, Cf, Omega, Gamma)
    p_obs = np.percentile(np.abs(observed_mr[mask]), kl_percentile)
    p_r1 = np.percentile(np.abs(R1_unit[mask]), kl_percentile)
    kl_global = p_obs / p_r1 if p_r1 > 0 else 1.0

    R0_global = nominal_migration_rate(curvature, W, kl_global)
    R1_global = predicted_migration_rate(R0_global, s, D, Cf, Omega, Gamma)

    # --- Step 4: auto-estimate window length ---
    if window_length is None:
        window_length = _estimate_window_length(curvature, s)

    # --- Step 5: sliding window local kl ---
    half_win = window_length / 2.0
    step = window_length / 2.0  # 50% overlap
    s_min, s_max = s[0], s[-1]

    window_centers = []
    kl_values = []
    sc = s_min + half_win
    while sc <= s_max - half_win:
        in_window = (s >= sc - half_win) & (s <= sc + half_win) & mask
        n_pts = np.sum(in_window)
        if n_pts >= min_window_points:
            p_obs_local = np.percentile(np.abs(observed_mr[in_window]),
                                        kl_percentile)
            p_r1_local = np.percentile(np.abs(R1_unit[in_window]),
                                       kl_percentile)
            kl_local = p_obs_local / p_r1_local if p_r1_local > 0 else kl_global
            window_centers.append(sc)
            kl_values.append(kl_local)
        sc += step

    window_centers = np.array(window_centers)
    kl_values = np.array(kl_values)

    # --- Step 6: interpolate to every point ---
    if len(window_centers) >= 2:
        kl_local = np.interp(s, window_centers, kl_values)
    elif len(window_centers) == 1:
        kl_local = np.full_like(s, kl_values[0])
    else:
        # Not enough data for windowing — fall back to global kl
        kl_local = np.full_like(s, kl_global)

    # --- Step 7: compute R1 with local kl ---
    # R1_unit was computed with kl=1, so R1_local = kl_local * R1_unit
    R1_local = kl_local * R1_unit

    return {
        'kl_local': kl_local,
        'kl_window_centers': window_centers,
        'kl_window_values': kl_values,
        'Cf': Cf,
        'D': D,
        'W': W,
        'R1_local': R1_local,
        'R1_global': R1_global,
        'kl_global': kl_global,
        'window_length': window_length,
        'observed_mr': observed_mr,
        's': s,
    }

calibrate_pair_local(results, pair_idx, use_stable_segments=True, **kwargs)

Calibrate local kl for a single pair from analysis results.

Thin wrapper around calibrate_local_kl that extracts arrays from the results dictionary (analogous to calibrate_pair).

Parameters:

Name Type Description Default
results dict

Results from analyze_river_pairs_filtered or analyze_segment_group.

 required
pair_idx int

Index of the pair.

 required
use_stable_segments bool

If True, mask out high-variance regions before calibration (default True).

 True
**kwargs

Forwarded to calibrate_local_kl.

 {}

Returns:

Name Type Description
result dict

Output of calibrate_local_kl, plus 'date1', 'date2', 'time_gap_years', 'pair_idx'.
Source code in rivabar/prediction.py 
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
def calibrate_pair_local(results, pair_idx, use_stable_segments=True,
                         **kwargs):
    """
    Calibrate local kl for a single pair from analysis results.

    Thin wrapper around ``calibrate_local_kl`` that extracts arrays from
    the ``results`` dictionary (analogous to ``calibrate_pair``).

    Parameters
    ----------
    results : dict
        Results from ``analyze_river_pairs_filtered`` or
        ``analyze_segment_group``.
    pair_idx : int
        Index of the pair.
    use_stable_segments : bool, optional
        If True, mask out high-variance regions before calibration
        (default True).
    **kwargs
        Forwarded to ``calibrate_local_kl``.

    Returns
    -------
    result : dict
        Output of ``calibrate_local_kl``, plus ``'date1'``, ``'date2'``,
        ``'time_gap_years'``, ``'pair_idx'``.
    """
    pair_info = results['pair_info'][pair_idx]
    curvature = results['curvatures'][pair_idx]
    distances = results['migration_distances'][pair_idx]
    s = results['along_channel_distances'][pair_idx]
    time_gap_years = pair_info['time_gap_years']
    W = np.mean(pair_info['width1'])
    observed_mr = distances / time_gap_years

    # Restrict calibration statistics to stable segments if requested.
    # The arrays passed to calibrate_local_kl are always the full contiguous
    # ones: slicing out high-variance regions before the convolution would
    # concatenate spatially disjoint points and distort the kernel.
    mask = _stable_segment_mask(pair_info, len(curvature)) if use_stable_segments else None

    cal = calibrate_local_kl(curvature, observed_mr, s, W,
                             mask=mask, **kwargs)

    # Everything is already on the full extent; keep the *_full keys as
    # aliases for backward compatibility with existing consumers
    cal['kl_local_full'] = cal['kl_local']
    cal['R1_local_full'] = cal['R1_local']
    cal['R1_global_full'] = cal['R1_global']
    cal['observed_mr_full'] = observed_mr
    cal['s_full'] = s

    cal['date1'] = pair_info['date1']
    cal['date2'] = pair_info['date2']
    cal['time_gap_years'] = time_gap_years
    cal['pair_idx'] = pair_idx
    return cal

calibrate_segment_local(results, df, pair_class_filter='good', window_length=None, time_window_years=5.0, time_step_years=None, use_stable_segments=True, min_time_gap_years=0, reference_river=None, delta_s=100, **kwargs)

Aggregate spatially-varying kl profiles across pairs with optional temporal windowing.

For each qualifying pair, calibrate_pair_local produces a kl(s) profile. These profiles are DTW-aligned to a reference centerline and aggregated (median, IQR). With temporal windowing, a sliding time window selects subsets of pairs, yielding kl(s, t).

Parameters:

Name Type Description Default
results dict

Results from analyze_river_pairs_filtered.

 required
df DataFrame

Analysis DataFrame with pair_class column.

 required
pair_class_filter str or list of str

Pair class(es) to include (default 'good').

 'good'
window_length float or None

Spatial window in metres (auto-estimated if None).

 None
time_window_years float or None

Width of the temporal window in years. If None, all pairs are used (no temporal windowing).

 5.0
time_step_years float or None

Step between temporal window centres (default: half of time_window_years).

 None
use_stable_segments bool

Exclude high-variance regions (default True).

 True
min_time_gap_years float

Minimum pair time gap (default 0).

 0
reference_river River or None

River object to use as the spatial reference for the common s-grid. If None, uses the last river in results['filtered_rivers'].

 None
delta_s float

Resampling interval for the reference s-grid (default 100 m).

 100
**kwargs

Forwarded to calibrate_local_kl.

 {}

Returns:

Name Type Description
result dict

Always contains: - 's_reference': common along-channel distance grid - 'x_reference', 'y_reference': reference centerline coords - 'Cf_global': global Cf used - 'window_length': spatial window used

Without temporal window (time_window_years=None): - 'kl_local_median': median kl(s) profile - 'kl_local_25', 'kl_local_75': IQR bounds - 'n_pairs': number of pairs used - 'pair_indices': list of pair indices used

With temporal window: - 'time_centers': list of datetime window centres - 'kl_local_series': list of dicts, one per time centre, each with 'kl_median', 'kl_25', 'kl_75', 'n_pairs', 'pair_indices'
Source code in rivabar/prediction.py 
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
def calibrate_segment_local(results, df, pair_class_filter='good',
                            window_length=None,
                            time_window_years=5.0,
                            time_step_years=None,
                            use_stable_segments=True,
                            min_time_gap_years=0,
                            reference_river=None,
                            delta_s=100,
                            **kwargs):
    """
    Aggregate spatially-varying kl profiles across pairs with optional
    temporal windowing.

    For each qualifying pair, ``calibrate_pair_local`` produces a kl(s)
    profile.  These profiles are DTW-aligned to a reference centerline
    and aggregated (median, IQR).  With temporal windowing, a sliding
    time window selects subsets of pairs, yielding kl(s, t).

    Parameters
    ----------
    results : dict
        Results from ``analyze_river_pairs_filtered``.
    df : DataFrame
        Analysis DataFrame with ``pair_class`` column.
    pair_class_filter : str or list of str
        Pair class(es) to include (default ``'good'``).
    window_length : float or None
        Spatial window in metres (auto-estimated if *None*).
    time_window_years : float or None
        Width of the temporal window in years.  If *None*, all pairs are
        used (no temporal windowing).
    time_step_years : float or None
        Step between temporal window centres (default: half of
        *time_window_years*).
    use_stable_segments : bool, optional
        Exclude high-variance regions (default True).
    min_time_gap_years : float, optional
        Minimum pair time gap (default 0).
    reference_river : River or None
        River object to use as the spatial reference for the common
        s-grid.  If *None*, uses the last river in
        ``results['filtered_rivers']``.
    delta_s : float, optional
        Resampling interval for the reference s-grid (default 100 m).
    **kwargs
        Forwarded to ``calibrate_local_kl``.

    Returns
    -------
    result : dict
        Always contains:
        - 's_reference': common along-channel distance grid
        - 'x_reference', 'y_reference': reference centerline coords
        - 'Cf_global': global Cf used
        - 'window_length': spatial window used

        Without temporal window (``time_window_years=None``):
        - 'kl_local_median': median kl(s) profile
        - 'kl_local_25', 'kl_local_75': IQR bounds
        - 'n_pairs': number of pairs used
        - 'pair_indices': list of pair indices used

        With temporal window:
        - 'time_centers': list of datetime window centres
        - 'kl_local_series': list of dicts, one per time centre, each
          with 'kl_median', 'kl_25', 'kl_75', 'n_pairs', 'pair_indices'
    """
    from .utils import get_width_and_curvature

    if isinstance(pair_class_filter, str):
        pair_class_filter = [pair_class_filter]

    # --- Identify qualifying pairs ---
    qualifying = df.index[df['pair_class'].isin(pair_class_filter)].tolist()
    if min_time_gap_years > 0:
        qualifying = [idx for idx in qualifying
                      if results['pair_info'][idx]['time_gap_years']
                      >= min_time_gap_years]
    if not qualifying:
        raise ValueError("No qualifying pairs found")

    # --- Build reference centerline ---
    if reference_river is None:
        filtered_rivers = results.get('filtered_rivers', [])
        if not filtered_rivers:
            raise ValueError(
                "results does not contain 'filtered_rivers'. "
                "Pass a River object via the reference_river parameter.")
        reference_river = filtered_rivers[-1]
    ref_river = reference_river

    pixel_size = kwargs.pop('pixel_size', None)
    x_ref, y_ref, s_ref, w_ref, _, _ = get_width_and_curvature(
        ref_river, delta_s=delta_s, pixel_size=pixel_size)

    # --- Calibrate each pair and remap to reference ---
    pair_calibrations = {}
    for idx in qualifying:
        try:
            cal = calibrate_pair_local(results, idx,
                                       use_stable_segments=use_stable_segments,
                                       window_length=window_length,
                                       **kwargs)
            coords = results['centerline_coords'][idx]
            kl_profile = cal.get('kl_local_full', cal['kl_local'])
            s_profile = cal.get('s_full', cal['s'])

            kl_on_ref = _map_to_reference_s(
                coords['x1'], coords['y1'], s_profile, kl_profile,
                x_ref, y_ref, s_ref)

            pair_calibrations[idx] = {
                'kl_on_ref': kl_on_ref,
                'cal': cal,
                'date1': cal['date1'],
                'date2': cal['date2'],
            }
            if window_length is None:
                window_length = cal['window_length']
        except Exception as e:
            print(f"  Pair {idx}: local calibration failed - {e}")

    if not pair_calibrations:
        raise ValueError("All pair calibrations failed")

    first_cal = next(iter(pair_calibrations.values()))['cal']
    Cf_global = first_cal['Cf']

    output = {
        's_reference': s_ref,
        'x_reference': x_ref,
        'y_reference': y_ref,
        'Cf_global': Cf_global,
        'window_length': window_length,
    }

    def _aggregate_pairs(indices):
        """Stack kl profiles for given pair indices and return stats."""
        profiles = []
        for idx in indices:
            kl = pair_calibrations[idx]['kl_on_ref']
            if not np.all(np.isnan(kl)):
                profiles.append(kl)
        if not profiles:
            n = len(s_ref)
            return {
                'kl_median': np.full(n, np.nan),
                'kl_25': np.full(n, np.nan),
                'kl_75': np.full(n, np.nan),
                'n_pairs': 0,
                'pair_indices': [],
            }
        stack = np.array(profiles)
        return {
            'kl_median': np.nanmedian(stack, axis=0),
            'kl_25': np.nanpercentile(stack, 25, axis=0),
            'kl_75': np.nanpercentile(stack, 75, axis=0),
            'n_pairs': len(profiles),
            'pair_indices': list(indices),
        }

    if time_window_years is None:
        # --- No temporal windowing: aggregate all ---
        agg = _aggregate_pairs(list(pair_calibrations.keys()))
        output['kl_local_median'] = agg['kl_median']
        output['kl_local_25'] = agg['kl_25']
        output['kl_local_75'] = agg['kl_75']
        output['n_pairs'] = agg['n_pairs']
        output['pair_indices'] = agg['pair_indices']
    else:
        # --- Temporal windowing ---
        if time_step_years is None:
            time_step_years = time_window_years / 2.0

        pair_midpoints = {}
        for idx, info in pair_calibrations.items():
            d1 = info['date1']
            d2 = info['date2']
            if isinstance(d1, str):
                d1 = pd.Timestamp(d1)
            if isinstance(d2, str):
                d2 = pd.Timestamp(d2)
            pair_midpoints[idx] = d1 + (d2 - d1) / 2

        all_midpoints = sorted(pair_midpoints.values())
        t_min = all_midpoints[0]
        t_max = all_midpoints[-1]

        half_win = pd.Timedelta(days=time_window_years * 365.25 / 2)
        step = pd.Timedelta(days=time_step_years * 365.25)

        time_centers = []
        kl_series = []
        tc = t_min + half_win
        while tc <= t_max - half_win:
            in_window = [idx for idx, mp in pair_midpoints.items()
                         if abs((mp - tc).total_seconds())
                         <= half_win.total_seconds()]
            if in_window:
                agg = _aggregate_pairs(in_window)
                time_centers.append(tc)
                kl_series.append(agg)
            tc += step

        if not time_centers:
            agg = _aggregate_pairs(list(pair_calibrations.keys()))
            time_centers = [t_min + (t_max - t_min) / 2]
            kl_series = [agg]

        output['time_centers'] = time_centers
        output['kl_local_series'] = kl_series

    return output