tests/testthat/test-wmo-guidelines.r

In RcmdrPlugin.lfstat: Rcmdr Plug-in for Low Flow Analysis

context("Examples from 'Guidance of WMO Software Tool for Low-flow Analysis")
data("ngaruroro")
ng_seventies <- subset(ngaruroro, hyear %in% 1970:1979)
ng_eighties <- subset(ngaruroro, hyear %in% 1980:1989)
ng_nineties <- subset(ngaruroro, hyear %in% 1990:1999)

# chapter 2

test_that("2.5.4 Missing values", {
  x <- lfnacheck(ngaruroro)
  expect_equal(x[["total"]], 214)
  expect_equal(x[["percentage"]], 0.0157145)
  expect_equal(names(x), c("total", "percentage", "hydrologicalyear", "duration"))
})

test_that("2.6 Flow indices", {

  # mean flow
  mf <- meanflow(lfobj = ngaruroro, year = "any", yearly = FALSE)
  expect_equal(round(mf, 2), 17.24)


  # q95
  q95 <- Q95(lfobj = ngaruroro, year = "any", yearly = TRUE)
  q95expected <- data.frame(
    hyear = 1964:2001,
    flow = c(3.72250, 5.68980, 6.12395, 5.35860, 3.64350, 4.94860, 4.34840, 5.11780, 4.57975,
             3.20820, 3.33200, 5.35200, 5.64100, 4.82720, 2.82045, 3.46460, 7.82550, 6.43980,
             4.53500, 3.00260, 5.97355, 4.92660, 4.52320, 5.50500, 5.63350, 4.44860, 4.68220,
             4.36020, 5.11100, 4.94580, 3.92500, 5.69320, 7.50850, 4.65920, 3.82960, 4.85320,
             4.41500, 6.78805)
  )

  expect_equal(round(q95,5), q95expected)

  # mean annual minimum
  mam <- MAM(lfobj = ngaruroro, n = 7, year = "any", yearly = FALSE)
  expect_equal(round(mam, 6), 4.380613)


  # base flow index
  bfi <- BFI(lfobj = ngaruroro, year = "any", yearly = FALSE)
  expect_equal(round(bfi, 7), 0.5499149)


  # master recession curve
  mrc <- recession(lfobj = ngaruroro, method = "MRC", seglength = 7,
                   threshold = 70, peaklevel = 0.95,
                   thresbreaks = "fixed", plotMRC = FALSE)
  expect_equal(round(mrc, 5), 19.80968)


  # seasonality index
  si <- seasindex(lfobj = ngaruroro, Q = 95)
  siexpected <- list(theta = 1.0888, D = 63.2505, r = 0.8552)
  expect_equal(lapply(si, round, 4), siexpected)


  # seasonality ratio
  sr <- seasratio(lfobj = ngaruroro, breakdays = c("01/06","01/10"), Q = 95)
  expect_equal(unname(round(sr, 4)), 2.0614)
})


test_that("2.7 Graphics", {
  # not easy to write test for this chapter, but we can test if the functions
  # does not stop with an error
  expect_silent(hydrograph(lfobj = ngaruroro, amin = TRUE))

  # seasonal barchart
  expect_silent(sbplot(ngaruroro))

  # base flow plot
  expect_silent(bfplot(lfobj = ngaruroro, year = 1970))

  # Flow duration curve
  expect_silent(fdc(lfobj = ngaruroro, year = "any", breakdays = c(),
                    ylog = TRUE, xnorm = FALSE,
                    colors = TRUE, legend = TRUE, separate = FALSE))

  # streamflow deficit plot
  expect_silent(streamdefplot(lfobj = ngaruroro, year = 1970, threslevel = 70,
                              thresbreaks = "fixed",
                              breakdays = c("01/06","01/10")))

  # double mass curve
  sub1 <- subset(ngaruroro, hyear %in% 1985:1989)
  sub2 <- subset(ngaruroro, hyear %in% 1990:1995)
  expect_silent(dmcurve(x = sub1 , y = sub2))


})

test_that("2.8.2 Regional Frequency Analysis", {
  # L-moment ratio diagram
  expect_silent(rfaplot(lflist = list(ng_eighties, ng_nineties, ng_seventies),
                        n = 7))


  rfa <- rfa(lflist = list(ng_eighties=ng_eighties,
                           ng_nineties=ng_nineties,
                           ng_seventies=ng_seventies),
             n =7, dist = "wei")


  # Index values
  expect_equal(round(rfa$index, 6),
               c("ng_eighties" = 4.611629,
                 "ng_nineties" = 4.414414,
                 "ng_seventies" = 3.927957))

  # T-year region
  expect_equal(round(lmomRFA::regquant(0.01, rfa), 7),
               c("0.01" = 0.5998283))


})


test_that("2.8 Extreme value", {

  rp <- tyearsn(lfobj = ngaruroro, event = 100, n = 7,
                dist = c("wei", "gev", "ln3", "gum", "pe3"))

  rpexpected <- structure(c(2.596112, 2.596112, 2.475504,
                            2.635668, 2.461607),
                          .Dim = c(1L, 5L),
                          .Dimnames = structure(list(
                            `return period` = "100",
                            distribution = c("wei", "gevR", "ln3", "gum", "pe3")),
                            .Names = c("return period", "distribution")))

  expect_equal(round(rp$T_Years_Event,6), round(rpexpected,6))

})

test_that("3 Case study", {

  # figure 3.3
  expect_silent(sbplot(ngaruroro))

  # text: seasonality ratio
  sr <- seasratio(lfobj = ngaruroro, breakdays = c("01/06","01/10"), Q = 95)
  expect_equal(unname(round(sr, 2)),  2.06)

  # figure 3.4
  expect_silent(bfplot(lfobj = ngaruroro, year = 1973))

  # table 3.1
  mrc <- recession(lfobj = ngaruroro, method = "MRC", seglength = 7,
                   threshold = 70, peaklevel = 0.95, seasonbreakdays = c(),
                   thresbreaks = "fixed",thresbreakdays = c("01/06","01/10"),
                   plotMRC = FALSE, trimIRS = 0.1)

  expect_equal(round(mrc, 1), 19.8)

  irs <- recession(lfobj = ngaruroro, method = "IRS", seglength = 7,
                   threshold = 70, peaklevel = 0.95, seasonbreakdays = c(),
                   thresbreaks = "fixed",thresbreakdays = c("01/06","01/10"),
                   plotMRC = FALSE, trimIRS = 0.1)

  expect_equal(round(irs, 1), 20.7)

  # figure 3.5
  expect_silent(fdc(lfobj = ngaruroro, year = "any", breakdays = c(),
                    ylog = TRUE, xnorm = FALSE, colors = TRUE, legend = TRUE,
                    separate = FALSE))

  # table 3.2
  fdc <- fdc(lfobj = ngaruroro, year = "any", breakdays = c(), ylog = TRUE,
             xnorm = FALSE, colors = TRUE, legend = TRUE, separate = FALSE)

  expect_equal(round(fdc[96], 1), 4.4)
  expect_equal(round(fdc[91], 1), 5.3)
  expect_equal(round(fdc[71], 1), 8.4)

  # table 3.3
  mam1 <- MAM(lfobj = ngaruroro, n = 1, year = "any", yearly = FALSE)
  expect_equal(round(mam1, 1), 4.1)

  mam7 <- MAM(lfobj = ngaruroro, n = 7, year = "any", yearly = FALSE)
  expect_equal(round(mam7, 1), 4.4)

  mam30 <- MAM(lfobj = ngaruroro, n = 30, year = "any", yearly = FALSE)
  expect_equal(round(mam30, 1), 5.4)

  # table 3.4
  expect_warning(deficit_table <- lfstat:::streamdefRcmdr(
    lfobj = ngaruroro, pooling = 'IC', threslevel = 90, thresbreaks = 'fixed',
    breakdays = c("01/06","01/10"), MAdays = , tmin = 5, IClevel = 0.1,
    mindur = '5', minvol = '0.5%', table = 'all', plot = FALSE))

  deficit_table <- subset(deficit_table, event.no  %in% 47:55)

  tableexpected <- data.frame(event.no = c(47,48,49,51,55),
                              start = c("1972-11-13", "1972-11-25", "1972-12-17", "1973-01-26", "1973-04-30"),
                              time = c("1972-11-21", "1972-11-30", "1973-01-09", "1973-04-20", "1973-05-04"),
                              end = c("1972-11-21", "1972-11-30", "1973-01-09", "1973-04-20", "1973-05-04"),
                              volume = c(257618.9, 229374.7, 1547786.9, 10016913.6, 283953.6),
                              duration = c(9, 6, 24, 85, 5),
                              dbt = c(9, 6, 24, 79, 5),
                              qmin = c(4.652, 4.713, 3.822, 2.780, 4.234),
                              tqmin = c(1972-11-19, 1972-11-26, 1973-01-05, 1973-03-03, 1973-05-04)
  )

  expect_equal(deficit_table$duration, tableexpected$duration)
  expect_equal(deficit_table$dbt, tableexpected$dbt)
  expect_equal(deficit_table$qmin, tableexpected$qmin)

  #  figure 3.14 left
  expect_warning(ty <- tyearsn(lfobj = ngaruroro, event = 100, n = 1,
                               dist = c("wei", "gevR", "ln3", "gum", "pe3")))

  tyexpected <- structure(c(2.494827, 2.494827, 2.431114, 2.596344, 2.372669),
                          .Dim = c(1L, 5L),
                          .Dimnames = structure(list(
                            `return period` = "100",
                            distribution = c("wei", "gevR", "ln3", "gum", "pe3")),
                            .Names = c("return period", "distribution")))

  expect_equal(round(ty$T_Years_Event,6), round(tyexpected,6))


  # figure 3.14 right
  ty_weibull <- tyearsn(lfobj = ngaruroro, event = 100, n = 1, dist = c("wei"))

  ty_weibullexpected <- structure(2.494827,
                                  .Dim = c(1L, 1L),
                                  .Dimnames = structure(list(
                                    `return period` = "100",
                                    distribution = c("wei")),
                                    .Names = c("return period", "distribution")))

  expect_equal(round(ty_weibull$T_Years_Event,6), round(ty_weibullexpected,6))

})