Compound rating curve fitting

The compound power-law level-discharge relationship for two segments is given as folows:

\[q = \begin{cases} a_1 \, (h-b_1)^{c_1} & \text{ if } h < k, \\ a_2 \, (h-b_2)^{c_2} & \text{ if } h \geq k. \end{cases}\]

Algorithm

For a given k, the rating curve for large discharges is fitted to the gaugings whose level is greater than or equal to $k$. The parameters $(a_2, b_2, c_2)$ are therefore estimated. For levels less than $k$, the rating curve is fitted to the gaugings whose level is less than $k$ by imposing that curve passes through the point $(k, a_2(k-b_2)^{c_2})$ to ensure continuity between the two segments. The search for the optimal $k$ is performed with the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS).

See also rcfit.

Fit on gaugings

The function crcfit can be called with the gauging levels and discharges:

julia> crc = crcfit(h,q)CompoundRatingCurve

   for 0.0 ≤ h < 27.681999497140765
      RatingCurve
         data: Vector{Gauging}[105]
         Parameters:
            a = 64.08997280259226
            b = 26.450047818104405
            c = 1.9289081890688273

   for 27.681999497140765 ≤ h < Inf
      RatingCurve
         data: Vector{Gauging}[89]
         Parameters:
            a = 110.347859154496
            b = 26.7775249837853
            c = 1.4041696598408295

The function returns a CompoundRatingCurve type.

The quality of the fit can be visually assessed in the log space:

x = log.(data.Level)
y = log.(data.Discharge)

obs = layer(x=x, y=y, Geom.point)
model = layer(x->logdischarge(crc, exp(x)), log(26), log(32), Theme(default_color=colorant"red"))
plot(obs, model, Coord.cartesian(ymin=2, ymax=7))

and in the original space:

obs = layer(data, x=:Level, y=:Discharge, Geom.point)
model = layer(h->discharge(crc, h), 26, 32, Theme(default_color=colorant"red"))
plot(obs, model, Coord.cartesian(ymin=0, ymax=1000))

Discharge estimation

With the fitted compound rating curve, the discharge estimation at the level $h_0 = 29$ can be ontained with discharge:

julia> discharge(crc, 29)338.67391359427245

It is also possible to estimate the level corresponding to the discharge $q = 340$ with level:

julia> level(crc, 340)29.006193882896923

Parameter uncertainty

The compound rating curve parameter uncertainty can be estimated by bootstrap with the function cint. For the Sainte-Anne example, the 95% confidence intervals are estimated using 100 boostrap samples of the original gaugings:

res = cint(crc, nboot=100)
println(string("k ∈ [", res[1,1]," , ", res[2,1]," ]"))
println(string("a₁ ∈ [", res[1,2]," , ", res[2,2]," ]"))
println(string("b₁ ∈ [", res[1,3]," , ", res[2,3]," ]"))
println(string("c₁ ∈ [", res[1,4]," , ", res[2,4]," ]"))
println(string("a₂ ∈ [", res[1,5]," , ", res[2,5]," ]"))
println(string("b₂ ∈ [", res[1,6]," , ", res[2,6]," ]"))
println(string("c₂ ∈ [", res[1,7]," , ", res[2,7]," ]"))

k ∈ [27.3485250090447 , 27.90981702390382 ]
a₁ ∈ [54.77716067375237 , 69.09710047545379 ]
b₁ ∈ [26.37662449418904 , 26.493657997015337 ]
c₁ ∈ [1.8037910622187332 , 2.1241854173333756 ]
a₂ ∈ [95.18997757246255 , 118.19054825760479 ]
b₂ ∈ [26.67461000828817 , 26.826245942656225 ]
c₂ ∈ [1.3608020750913712 , 1.503991265383004 ]

Discharge uncertainty

The 95% confidence interval on the discharge estimation at $h_0 = 29$ can be obtained with pint(::CompoundRatingCurve, ::Real, ::Real, ::Real):

julia> pint(crc, 29)2-element Vector{Float64}:
 314.6338316255117
 364.55081501146617

For more details on how this uncertainty is estimated, see the description of pint(::RatingCurve, ::Real, ::Real, ::Real).

The confidence interval for the whole level range of the rating curve can be plotted as follows:

h₀ = range(minimum(data.Level), stop=32, length=100)

# Discharge estimation for each level
q̂₀ = discharge.(crc, h₀)

# 95% confidence intervals of each discharge estimation
res = pint.(crc, h₀)

# Lower bound of interval
qmin = getindex.(res,1)

# Upper bound of interval
qmax = getindex.(res,2)

# Plotting the interval and the gaugings
obs = layer(data, x=:Level, y=:Discharge, Geom.point)
model = layer(x=h₀, y=q̂₀, Geom.line,
    ymin = qmin, ymax = qmax, Geom.ribbon,
    Theme(default_color=colorant"red"))

plot(obs, model)

Fit quality assessment

The BIC (Bayesian Information Criterion) is an index of the quality of the curve fit to the gaugings. Assuming that the errors are normally distributed in the log space for each segment, the BIC of the compund curve in the log space can be obtained as follows:

$\operatorname{bic} = n_1 \log \hat\sigma_{e1}^2 + n_2 \log \hat\sigma_{e2}^2+ 6 \log n,$

where $n_1$ and $n_2$ correspond to the number of gaugings of the first and second segment respectively, $\hat\sigma_{e1}^2$ and $\hat\sigma_{e2}^2$ correspond to the variance of the errors in the log space for the first and second segment respectively, $n$ corresponds to the total number of gauging and the value 6 stands for the number of parameters.

The BIC of a fitted compound rating curve can be obtained with the function bic:

julia> bic(crc)-1262.6034314266365

Compound rating curve fit with a known breakpoint

If the breakpoint $k$ is known by the investigator, the rating curves of the two segments can be fitted with crcfit by specifying the optional argument k. For example, let suppose that the breakpoint for the Sainte-Anne data is at $k = 28$, then the compoud rating curve can be obtained as follows:

julia> crcfit(G, 28)CompoundRatingCurve

   for 0.0 ≤ h < 28.0
      RatingCurve
         data: Vector{Gauging}[122]
         Parameters:
            a = 68.03370837893993
            b = 26.484792318364956
            c = 1.8170551886172388

   for 28.0 ≤ h < Inf
      RatingCurve
         data: Vector{Gauging}[72]
         Parameters:
            a = 121.48662247776792
            b = 26.86182958889689
            c = 1.3543361612729679