oce: Analysis of Oceanographic Data

# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4
library(oce)

test_that("gappyIndex", {
    expect_equal(c(3:6, 103:106), gappyIndex(c(1, 101), 2, 4))
})

test_that("approx3d", {
    # Test values from the .c code, before converting to .cpp
    n <- 5
    x <- seq(0, 1, length.out=n)
    y <- seq(0, 1, length.out=n)
    z <- seq(0, 1, length.out=n)
    f <- array(1:n^3, dim=c(length(x), length(y), length(z)))
    # interpolate along a diagonal line
    m <- 10
    xout <- seq(0, 1, length.out=m)
    yout <- seq(0, 1, length.out=m)
    zout <- seq(0, 1, length.out=m)
    approx <- approx3d(x, y, z, f, xout, yout, zout)
    expect_equal(approx,
        c(1,  14.77777778, 28.55555556, 42.33333333, 56.11111111,
            69.88888889, 83.66666667, 97.44444444, 111.22222222, NA))
})

test_that("Coriolis", {
    f <- coriolis(45)
    expect_equal(f, 1.031261e-4, tolerance=1e-6)
})

test_that("despike", { # issue 1067
    min <- 10000
    max <- 20000
    x1 <- c(3715, 7546, 10903, 13386, 15196, 15371, 55748, 71488)
    x2 <- despike(x1, reference="trim", min=min, max=max, replace="reference")
    x3 <- x1
    x3[1:2] <- 10903 # result from approx() with rule=2
    x3[7:8] <- 15371 # result from approx() with rule=2
    expect_equal(x2, x3)
    x4 <- despike(x1, reference="trim", min=min, max=max, replace="NA")
    x5 <- x1
    x5[x5<min] <- NA
    x5[x5>max] <- NA
    expect_equal(x4, x5)
})

test_that("oceConvolve", {
    expect_equal(oceConvolve(c(rep(-1, 10), rep(1, 10)), c(1/4, 1/2, 1/4)),
        c(-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
            -1.0, -0.5, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
            1.0))
})

test_that("oceEdit", {
    data(ctd)
    # metadata
    a <- oceEdit(ctd, "waterDepth", 100)
    expect_equal(a@metadata$waterDepth, 100)
    b <- oceEdit(a, "waterDepth", 200)
    expect_equal(b@metadata$waterDepth, 200)
    c <- oceEdit(b, "metadata@waterDepth", 300)
    expect_equal(c@metadata$waterDepth, 300)
    # data
    scan <- ctd[["scan"]]
    ctd2 <- oceEdit(ctd, "scan", scan+1)
    expect_true("scan" %in% names(ctd2[["data"]]))
    expect_false("scan" %in% names(ctd2[["metadata"]]))
    expect_equal(ctd[["scan"]]+1, ctd2[["scan"]])
    # check "case 1A" of docs (cannot add non-extant item)
    expect_error(oceEdit(ctd, "noSuchThing", 1), "nothing named 'noSuchThing'")
    # check "case 2" of docs (can add non-extant item, if slot provided)
    ctd <- oceEdit(ctd, "metadata@newItem", 1)
    expect_true("newItem" %in% names(ctd[["metadata"]]))
    p <- ctd[["pressure"]]
    ctd <- oceEdit(ctd, action="x@data$pressure <- 0.1 + x@data$pressure")
    expect_equal(p + 0.1, ctd[["pressure"]])
})

test_that("get_bit (unused in oce)", {
    buf <- 0x3a
    bits <- unlist(lapply(7:0, function(i) do_get_bit(buf, i)))
    # NB. 'i' starts at rightmost bit
    expect_equal(c(0, 0, 1, 1, 1, 0, 1, 0), bits)
})

test_that("grad", {
    g <- grad(volcano)
    expect_equal(mean(g$g), 196.982876)

    expect_equal(mean(g$gx), -6.903335218)
    expect_equal(mean(g$gy), -7.009609949)
    expect_equal(g$g[1:3, 1:3],
        matrix(c(86, 86, 86, 91.08238029, 91.08238029,
                91.08238029, 91.08238029, 91.08238029,
                91.08238029), nrow=3))
    expect_equal(g$gx[1:3, 1:3],
        matrix(rep(86, 9), nrow=3))
    expect_equal(g$gy[1:3, 1:3],
        matrix(c(0, 0, 0, 30, 30, 30, 30, 30, 30), nrow=3))
})

test_that("gravity", {
    g <- gravity(45)
    expect_equal(g, 9.8, tolerance=1e-2)
})

test_that("integration", {
    x <- seq(0, 1, length.out=10)
    dx <- x[2] - x[1]
    y <- 2*x + 3*x^2
    A <- integrateTrapezoid(x, y)
    expect_equal(A, 2, tolerance=dx^2) # test for quadratic accuracy
})
test_that("integrateTrapezoid", {
    x <- seq(0, 1, length.out=10)
    y <- rep(1, length(x))
    expect_equal(1, integrateTrapezoid(x, y))
    expect_equal(4, integrateTrapezoid(x, y, xmin=-2, xmax=2))
    expect_equal(9, integrateTrapezoid(rep(1, 10)))
    x <- seq(0, 1, length.out=10)
    y <- 2*x + 3*x^2
    expect_equal(2, integrateTrapezoid(x, y), tolerance=0.01)
    expect_equal(integrateTrapezoid(x, y), integrateTrapezoid(y) * (x[2]-x[1]))
    expect_equal(c(0.0000000000000, 0.0144032921811, 0.0473251028807,
            0.0884773662551, 0.1378600823045, 0.1954732510288,
            0.2613168724280, 0.3353909465021, 0.4176954732510,
            0.5082304526749),
        integrateTrapezoid(x, y, "dA"))
    expect_equal(c(0.0000000000000, 0.0144032921811, 0.0617283950617,
            0.1502057613169, 0.2880658436214, 0.4835390946502,
            0.7448559670782, 1.0802469135802, 1.4979423868313,
            2.0061728395062),
        integrateTrapezoid(x, y, "cA"))
})

test_that("interpBarnes 1D", {
    # These tests are not in comparison to theory, or
    # known values; they simply ensure that results have not
    # changed since 2018-03-11, when the tests were devised.
    data(ctd)
    p <- ctd[["pressure"]]
    y <- rep(1, length(p)) # fake y data, with arbitrary value
    S <- ctd[["salinity"]]
    pg <- pretty(p, n=100)
    g <- interpBarnes(p, y, S, xg=pg, xr=1)
    expect_equal(g$zd[c(1, 10, 100)], c(29.91878482, 29.94385118, 31.44549220))
})

test_that("interpBarnes 2D", {
    # These tests are not in comparison to theory, or
    # known values; they simply ensure that results have not
    # changed since 2016-11-06, when the tests were devised.
    data(wind)
    u <- interpBarnes(wind$x, wind$y, wind$z)
    expect_equal(u$zg[1, 1], 30.962611975027)
    expect_equal(u$zg[5, 1], 20.93550551)
    expect_equal(u$zg[1, 5], 34.2550759)
    expect_equal(u$zg[10, 10], 27.042654784966)
})

test_that("magneticField() handles both POSIX times and dates", {
    A <- magneticField(-63.562, 44.640, as.POSIXct("2013-01-01", tz="UTC"), version=12)$declination
    B <- magneticField(-63.562, 44.640, as.Date("2013-01-01"), version=12)$declination
    expect_equal(A, B, tolerance=1e-8)
})

test_that("magneticField version 12 (why not perfect?)", {
    # test values from http://www.geomag.bgs.ac.uk/data_service/models_compass/wmm_calc.html
    # UPDATE March 3, 2020: I cannot test these old values because that
    # page now only works for present and future dates (and it's quite
    # hard to figure out, frankly).
    expect_equal(-17.976, magneticField(-63.562, 44.640, 2013, version=12)$declination,
        tolerance=0.001)
    expect_equal(67.562, magneticField(-63.562, 44.640, 2013, version=12)$inclination,
        tolerance=0.006) # Q: why does tol=0.001 fail?
    expect_equal(52096, magneticField(-63.562, 44.640, 2013, version=12)$intensity,
        tolerance=16) # Q: why does tol=1 fail?
})

test_that("magneticField version 13 (why not perfect?)", {
    # REF: http://www.geomag.bgs.ac.uk/data_service/models_compass/wmm_calc.html
    # version 13 by default as of oce "develop" branch date 2020-03-03
    mf <- magneticField(-63.562, 44.640, as.POSIXct("2020-03-03 00:00:00", tz="UTC"), version=13)
    mf2 <- magneticField(-63.562, 44.640, as.Date("2020-03-03", tz="UTC"), version=13)
    cbind(mf, mf2) # Q: why so much difference here?
    expect_equal(-16.972, mf$declination, tolerance=0.005) # Q: why does tol=0.001 fail?
    expect_equal(66.855, mf$inclination, tolerance=0.001)
    expect_equal(51498, mf$intensity, tolerance=3) # Q: why does tol=1 fail?
})

test_that("matchBytes with 2 bytes", {
    buf <- as.raw(c(0xa5, 0x11, 0xaa, 0xa5, 0x11, 0x00))
    expect_equal(c(1, 4), matchBytes(buf, 0xa5, 0x11))
})

test_that("matchBytes with 3 bytes", {
    buf <- as.raw(c(0xa5, 0x11, 0xaa, 0x12, 0xa5, 0x11, 0xaa, 0x99))
    expect_equal(c(1, 5), matchBytes(buf, 0xa5, 0x11, 0xaa))
})

test_that("matrixSmooth", {
    # Test for same values after rewriting the C code in C++.
    data(volcano)
    v <- matrixSmooth(volcano)
    ve <- matrix(c(100, 101, 102, 100, 101.1666667, 102.1666667, 101,
            101.8333333, 102.8333333), byrow=FALSE, nrow=3)
    expect_equal(ve, v[1:3, 1:3])
    expect_equal(mean(v), 130.1787262)
    expect_equal(mean(v[, 10]), 121.3429119)
    expect_equal(mean(v[10, ]), 127.9808743)
})

test_that("rotateAboutZ adp", {
    data(adp)
    v11 <- adp[["v"]][1, 1, ]
    angle <- 45
    adpRotated <- rotateAboutZ(adp, angle)
    theta <- angle * pi / 180
    m <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), byrow=TRUE, nrow=2)
    v11Rotated <- as.vector(m %*% cbind(v11[1:2]))
    expect_equal(v11Rotated, adpRotated[["v"]][1, 1, 1:2])
})

test_that("rotateAboutZ adv", {
    data(adv)
    v1 <- adv[["v"]][1, ]
    angle <- 45
    advRotated <- rotateAboutZ(adv, angle)
    theta <- angle * pi / 180
    m <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), byrow=TRUE, nrow=2)
    v1Rotated <- as.vector(m %*% cbind(v1[1:2]))
    expect_equal(v1Rotated, advRotated[["v"]][1, 1:2])
})

test_that("rotateAboutZ cm", {
    data(cm)
    u1 <- cm[["u"]][1]
    v1 <- cm[["v"]][1]
    angle <- 45
    cmRotated <- rotateAboutZ(cm, angle)
    theta <- angle * pi / 180
    m <- matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), byrow=TRUE, nrow=2)
    uvRotated <- as.vector(m %*% cbind(c(u1, v1)))
    expect_equal(uvRotated, c(cmRotated[["u"]][1], cmRotated[["v"]][1]))
})

test_that("runlm", {
    # Test for same values after rewriting the C code in C++.
    x <- 1:8
    y <- c(1.208669331, 1.409418342, 1.594642473, 1.757356091,
        1.891470985, 1.992039086, 2.055449730, 2.079573603)
    # default L
    r <- runlm(x, y)
    expect_equal(4, sum(c("dydx", "L", "x", "y") %in% names(r)))
    expect_equal(r$y, c(1.210171739, 1.402871561, 1.581586126,
            1.740768514, 1.872013473, 1.970487349,
            2.040003801, 2.083375549))
    expect_equal(r$dydx, c(0.19524018261, 0.18668884855, 0.17210118250,
            0.14683764457, 0.11611882136, 0.08116937564,
            0.05596691825, 0.03806434013))
    expect_equal(r$L, 6)
    # specified L
    r <- runlm(x, y, L=4)
    expect_equal(r$dydx, c(0.2007490110, 0.1929865710, 0.1739688745,
            0.1484142560, 0.1173414975, 0.0819893725,
            0.0437672585, 0.0241238730))
    expect_equal(r$y, c(1.208669331, 1.405537122, 1.589014845,
            1.750206410, 1.883084287, 1.982749722,
            2.045628037, 2.079573603))
    expect_equal(r$L, 4)
})


test_that("snakeToCamel", {
    expect_equal(snakeToCamel("PARAMETER_DATA_MODE"),   "parameterDataMode")
    expect_equal(snakeToCamel("PARAMETER"),             "parameter")
    expect_equal(snakeToCamel("HISTORY_QCTEST"),        "historyQctest")
    expect_equal(snakeToCamel("HISTORY_QCTEST", "QC"),  "historyQCTest")
    expect_equal(snakeToCamel("PROFILE_DOXY_QC"),       "profileDoxyQc")
    expect_equal(snakeToCamel("PROFILE_DOXY_QC", "QC"), "profileDoxyQC")
})

test_that("time-series filtering", {
    # Check against some matlab results.
    b <- rep(1, 5)/5
    a <- 1
    x <- seq(1, 4, by=0.2)
    matlab.res <- c(0.2000, 0.4400, 0.7200, 1.0400, 1.4000, 1.6000, 1.8000, 2.0000, 2.2000,
        2.4000, 2.6000, 2.8000, 3.0000, 3.2000, 3.4000, 3.6000)
    expect_equal(matlab.res, oce.filter(x, a, b))
    # Check against old values.
    b <- rep(1, 5)/5
    a <- 1
    x <- seq(0, 10)
    y <- ifelse(x == 5, 1, 0)
    f <- oce.filter(y, a, b, zero.phase=TRUE)
    expect_equal(f, c(0.00, 0.04, 0.08, 0.12, 0.16, 0.20, 0.16, 0.12,
            0.08, 0.04, 0.00))
})

test_that("gps time", {
    # The GPS test value was calculated as follows:
    # https://www.labsat.co.uk/index.php/en/gps-time-calculator
    # gives week=604 and sec=134336 (for the indicated date), IGNORING
    # leap seconds. However,
    # https://confluence.qps.nl/display/KBE/UTC+to+GPS+Time+Correction#UTCtoGPSTimeCorrection-UTC(CoordinatedUniversalTime)
    # indicates that a 15-second correction was needed for GPS to UTC, so
    # we do that in the test value.
    expect_equal(numberAsPOSIXct(cbind(604, 134336+15), type="gps"),
        as.POSIXct("2011-03-21 13:18:56", tz="UTC"))
    expect_equal(numberAsPOSIXct(cbind(604, 134336+15, 1), type="gps"),
        as.POSIXct("2011-03-21 13:18:56", tz="UTC"))
    # Check internal computation of leap seconds
    ls <- as.POSIXct(.leap.seconds, tz="UTC")
    t <- numberAsPOSIXct(cbind(566, 345615), type="gps")
    expect_equal(15, sum(as.POSIXct("1980-01-01",tz="UTC") < ls & ls < t))
    # For 3-column mode (and this test value) see https://github.com/dankelley/oce/issues/2077
    expect_equal(numberAsPOSIXct(cbind(2250, 502834.494144, 0), type="gps"),
        as.POSIXct("2023-02-24 19:40:16.494", tz="UTC"))
})

test_that("matlab time", {
    # Matlab times; see http://www.mathworks.com/help/matlab/ref/datenum.html
    expect_equal(numberAsPOSIXct(719529, "matlab"), ISOdatetime(1970, 1, 1, 0, 0, 0, tz="UTC"))
    mt <- 7.362007209411687e5
    expect_equal(as.numeric(numberAsPOSIXct(mt, "matlab", tz="UTC")),
        as.numeric(as.POSIXct("2015-08-24 17:18:09", tz="UTC")), tolerance=1)
})


test_that("excel time", {
    # Excel time. I created the test value by entering "Jul 1, 2019" into
    # excel, then copying to a new cell with "paste special" set to
    # "value".
    expect_equal(numberAsPOSIXct(43647.0, "excel"), ISOdatetime(2019, 07, 01, 0, 0, 0, tz="UTC"))
    # Now, check around the erroneous leap-day in 1900, for which oce
    # returns a time of NA, to tell the user that 1900 was not a leap
    # year (i.e. for consistency with what user would get by trying to
    # specify that date in ISOdatetime() or other functions).
    expect_equal(numberAsPOSIXct(59, "excel"), ISOdatetime(1900, 2, 28, 0, 0, 0, tz="UTC"))
    expect_true(is.na(numberAsPOSIXct(60, "excel")))
    expect_equal(numberAsPOSIXct(61, "excel"), ISOdatetime(1900, 03, 01, 0, 0, 0, tz="UTC"))
    expect_equal(numberAsPOSIXct(367, "excel"), ISOdatetime(1901, 01, 01, 0, 0, 0, tz="UTC"))
    expect_equal(numberAsPOSIXct(368, "excel"), ISOdatetime(1901, 01, 02, 0, 0, 0, tz="UTC"))
})

test_that("ncep1 and ncep2 time", {
    # NCEP1 times; test value from
    # http://coastwatch.pfeg.noaa.gov/erddap/convert/time.html?isoTime=2015-09-04T12%3A00%3A00Z&units=hours+since+1800-01-01
    expect_equal(as.numeric(numberAsPOSIXct(1890564, "ncep1")),
        as.numeric(as.POSIXct("2015-09-04 12:00:00", tz="UTC")), tolerance=1)
    # NCEP2 times; see http://www.esrl.noaa.gov/psd/data/gridded/faq.html#3
    # and also https://github.com/dankelley/oce/issues/739, the latter
    # documenting what is essentially a kludge for this to work.
    expect_equal(as.numeric(numberAsPOSIXct(725738, "ncep2")),
        as.numeric(as.POSIXct("1988-01-01 00:00:00", tz="UTC")), tolerance=1)
})