Segmented Power Law

Segmented Power Law#

Open In Colab Binder

There are several approaches to fitting a stage-discharge rating curve. This notebook demonstrates an simple way to approximate the classic approach, which uses a segmented power law.

# Uncomment below to setup Google Colab. It will take a minute or so.
# %%capture
# !pip install pymc==5.1.2
# %env MKL_THREADING_LAYER=GNU
# !pip install git+https://github.com/thodson-usgs/ratingcurve.git

%load_ext autoreload
%autoreload 2

from ratingcurve.ratings import PowerLawRating

Load Data#

Begin by loading the Green Channel dataset

from ratingcurve import data

df = data.load('green channel')
df.head()
datetime stage q q_sigma
0 2020-05-21 14:13:41 [UTC-07:00] 7.04 12199.342 199.172931
1 2020-04-16 14:55:31 [UTC-07:00] 4.43 4921.953 95.425619
2 2020-03-04 13:54:10 [UTC-07:00] 2.99 2331.665 61.860500
3 2020-03-04 13:16:51 [UTC-07:00] 2.94 2289.220 47.886745
4 2020-01-23 11:04:32 [UTC-07:00] 2.96 2408.210 99.522964

and plotting the observations.

ax = df.plot.scatter(x='q', y='stage', marker='o')
ax.set_xlabel("Discharge (cfs)")
ax.set_ylabel("Stage (ft)")

Text(0, 0.5, 'Stage (ft)')
../_images/89e8bd33d94061f5c6bf8adf88a1cd36d722c50c0d5e6ff9e9808ffea66f7d9a.png

Setup model#

Now, setup the rating model. This make take a minute the first time while the model compiles but will be faster on subsequent runs.

powerrating = PowerLawRating(segments=2,
                             prior={'distribution': 'uniform'})

There are of variety of ways to adjust the optimization. Here we’ll use the defaults, which uses ADVI and runs for 200,000 iterations, though the model should coverge well before that.

trace = powerrating.fit(q=df['q'],
                        h=df['stage'],
                        q_sigma=df['q_sigma'])

/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



Convergence achieved at 78800





Interrupted at 78,799 [39%]: Average Loss = 22.626





Sampling: [a, b, hs_, model_q, sigma]





Sampling: [model_q]





/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')









Once fit, we can plot the rating curve.






powerrating.plot()









Sampling: [model_q]





/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')





../_images/2a684ccfaf75382ad25a7c70da7614d37dfa3950d09cca79336c08a64a5896f2.png





or as a table of stage-discharge values.






table = powerrating.table()
table.head()









Sampling: [model_q]





/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')









  
    

      
      stage
      discharge
      median
      gse
    

  
  
    

      0
      2.21
      1384.174800
      1384.088645
      1.017904
    

    

      1
      2.22
      1396.732699
      1396.265913
      1.017979
    

    

      2
      2.23
      1408.924423
      1408.916273
      1.018026
    

    

      3
      2.24
      1420.806259
      1420.354262
      1.017792
    

    

      4
      2.25
      1433.018784
      1432.809323
      1.017745
    

  











Exercise#


What happens if we choose the wrong number of segments?
Increase the number of segments by one and rerun the model.
In fact, we can use this to select the correct number of segments,
which is demonstrated in the model evaluation notebook.






Simulated Example#


This example uses a simulated rating curve, which allows you to test how changing the number of segments affects the rating curve.


First, load the ‘3-segment simulated’ tutorial dataset.






sim_df = data.load('3-segment simulated')









This rating contains observations of every 0.01 inch. increment in stage, which is much more than we’d have for a natural rating.
Try sampling to n_sample=15 or n_sample=30 and see how that affects the model fit.






# subsample the simulated rating curve
n_sample = 30
df = sim_df.sample(n_sample, random_state=12345)

ax = sim_df.plot(x='q', y='stage', color='grey', ls='-', legend=False)
df.plot.scatter(x='q', y='stage', marker='o', color='blue', ax=ax)
ax.set_xlabel("Discharge (cfs)")
ax.set_ylabel("Stage (ft)");








../_images/fad8188929d4cdb129aad43b5b8e13991335d7ce5c774d95ccaf6404e756a6bd.png





Setup a rating model with 3 segments






powerrating = PowerLawRating(segments=3,
                             prior={'distribution': 'uniform'},
                             # prior={'distribution':'normal', 'mu':[5, 8, 11], 'sigma':[1, 1, 0.2]}
                            )









now fit the model using ADVI






trace = powerrating.fit(q=df['q'],
                        h=df['stage'],
                        q_sigma=None,
                        method='advi')









/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')





Finished [100%]: Average Loss = -48.506





Sampling: [a, b, hs_, model_q, sigma]





Sampling: [model_q]





/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')









and visualize the results.






import matplotlib.pyplot as plt

fig, ax = plt.subplots()
powerrating.plot(ax=ax)

# plot the original data for comparison
sim_df.plot(x='q', y='stage', color='red', ls=':', legend=False, ax=ax)









Sampling: [model_q]





/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')



/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/rich/live.py:260: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')





<Axes: xlabel='q', ylabel='Stage'>




../_images/cdba88417ac48235deab5046faf5681ba8fdfb01c691b3615668d7c2a7d02197.png









powerrating.summary(var_names=["b", "a", "sigma", "hs"])










  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      b[0]
      0.98232
      0.00209
      0.98
      0.99
      10224
      9916
      nan
      2.1e-05
      1.5e-05
    

    

      b[1]
      0.3518
      0.0075
      0.34
      0.36
      9527
      9409
      nan
      7.7e-05
      5.4e-05
    

    

      b[2]
      0.048
      0.0569
      -0.043
      0.14
      9927
      9484
      nan
      0.00057
      0.0004
    

    

      a
      -1.00683
      0.00279
      -1
      -1
      9189
      9511
      nan
      2.9e-05
      2e-05
    

    

      sigma
      0.01561
      0.00208
      0.012
      0.019
      9868
      9445
      nan
      2.1e-05
      1.5e-05
    

    

      hs[0, 0]
      4.92776
      0.00122
      4.9
      4.9
      9044
      9624
      nan
      1.3e-05
      9.1e-06
    

    

      hs[1, 0]
      10.5474
      0.0314
      10
      11
      10148
      9789
      nan
      0.00031
      0.00023
    

    

      hs[2, 0]
      12.315
      0.167
      12
      12
      9973
      9537
      nan
      0.0017
      0.0025
    

  





Fitting this model can be tricky.
The most common issue is a poor initialization of the breakpoints.
A fix is under development, but for now, try


    

  1. reinitializing the model PowerLawRating();
    2. 

  2. increasing the number of iterations for the fitting algorithm fit(n=300_000);
    3. 

  3. a prior on the breakpoints, example, try prior={'distribution':'normal', 'mu':[5, 9.5, 10.5], 'sigma':[1, 1, 0.2]}), which implies we know the true breakpoint within +-0.5 ft; or
    4. 

  4. fitting the model with NUTS fit(method='nuts').
    5. 







%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,xarray









Last updated: Tue, 01 Sep 2026

Python implementation: CPython
Python version       : 3.12.14
IPython version      : 9.17.1

pytensor: 3.3.0
xarray  : 2026.7.0

matplotlib : 3.11.1
ratingcurve: 0.1.dev1+ga5a0fe6ba

Watermark: 2.6.0