tests/testthat/test_sw.R

In oce: Analysis of Oceanographic Data

# vim:textwidth=140:expandtab:shiftwidth=4:softtabstop=4

library(oce)
# Table of contents.
#  1. rho and sigma
#  2. potential temperature and temperature scales
#  3. Absolute Salinity and Conservative Temperature
#  4. sound speed
#  5. freezing temperature
#  6. specific heat
#  7. adiabatic lapse rate
#  8. alpha and beta
#  9. swSTrho
# 10. sound absorption (unesco only)
# 11. viscosity (unesco only)
# 12. thermal conductivity (unesco only)
# 13. electrical conductivity
# 14. depth and pressure [FIXME: INCOMPLETE]
# 15. spiciness

# spec vol anom and dens anom
# pressure to depth

# The UNESCO properties are generally tested against UNESCO documents [1].
# The GSW properties are more complicated because its check values are generally
# stated in terms of Absolute Salinity (SA), not Practical Salinity (SP).
# However, we already know that the GSW functions are very well tested
# through tests of the 'gsw' package, so all we really need is a single
# test case for each function, to ensure that eos="gsw" is interpreted
# correctly, etc.
#
# References:
# [1] N. P. Fofonoff and R. C. Millard Jr., 1983.
#     Algorithms for computation of fundamental properties of seawater.
#     UNESCO technical papers in marine science, vol 44.
#     http://www.jodc.go.jp/info/ioc_doc/UNESCO_tech/059832eb.pdf
# [2] Trevor J. McDougall, 1987. Neutral Surfaces, Journal of Physical
#     Oceanography, volume 17, pages 1950-1964.

test_that("1. rho and sigma", {
    # 1. rho and sigma
    # 1.1 UNESCO rho [1 p19]. Note that we must
    # convert to the T68 temperature scale, which was in use at the time
    # that [1] was written.
    sal <- c(0,   0,   0,   0,  35,  35,  35,  35)
    tem <- T90fromT68(c(5,   5,  25,  25,   5,   5,  25,  25))
    pre <- c(0, 1e4,   0, 1e4,   0, 1e4,   0, 1e4)
    rho <- c(999.96675, 1044.12802, 997.04796, 1037.90204, 1027.67547, 1069.48914, 1023.34306, 1062.53817)
    expect_equal(swRho(sal, tem, pre, eos="unesco"), rho)
    # Does swRho() work for a ctd object?
    ctd <- as.ctd(sal, tem, pre)
    expect_equal(swRho(ctd, eos="unesco"), rho)
    # Does accessor work for a ctd object?
    oldEOS <- getOption("oceEOS")
    options(oceEOS="unesco")
    expect_equal(ctd[["density", eos="unesco"]], rho)
    options(oceEOS=oldEOS)
    # check sigma from this
    expect_equal(swRho(sal, tem, pre, eos="unesco")-1000, swSigma(sal, tem, pre, eos="unesco"))
    # 1.2 GSW
    # Since gsw_ functions are tested in the gsw package, we just need a consistency check.
    longitude <- 188
    latitude <- 4
    SP <- 35
    t <- 10
    p <- 1000
    SA <- gsw::gsw_SA_from_SP(SP, p, longitude, latitude)
    CT <- gsw::gsw_CT_from_t(SA, t, p)
    # Test density.
    rhoGSW <- gsw::gsw_rho(SA, CT, p)
    rho <- swRho(SP, t, p, longitude, latitude, eos="gsw")
    expect_equal(rhoGSW, rho)
    # Now use density to test sigma (not provided by gsw).
    sigma <- swSigma(SP, t, p, longitude, latitude, eos="gsw")
    expect_equal(rhoGSW-1000, sigma)
    # The following was hard-coded using values from GSW3.03, and it failed with GSW3.05.
    expect_equal(30.818302, swSigma(35, T90fromT68(13), 1000, eos="unesco"), tolerance=0.000001)
    # The sigmaT tests are not from definitive test values, but only provide a check
    # against future changes.
    expect_equal(swSigmaT(35, T90fromT68(13), 1000, eos="unesco"), 26.393538, tolerance=0.000001)

    # Tests from issue 1904
    data(section)
    stn <- section[["station", 100]] # 4.3km deep
    # Ensure that the sigmaTheta is the same whether called directly or with [[
    expect_equal(swSigmaTheta(stn, eos="unesco"), stn[["sigmaTheta", "unesco"]])
    expect_equal(swSigmaTheta(stn, eos="gsw"), stn[["sigmaTheta", "gsw"]])

    # The oceEOS option must be obeyed.
    options(oceEOS="unesco")
    expect_equal(stn[["sigmaTheta"]], stn[["sigma0", "unesco"]])
    options(oceEOS="gsw")
    #(before issue1933) expect_equal(stn[["sigmaTheta"]], stn[["sigma0", "gsw"]])
    expect_false(identical(stn[["sigmaTheta"]], stn[["sigma0", "gsw"]]))

    # sigmaTheta and sigma0 should match within each EOS.
    expect_equal(stn[["sigmaTheta", "unesco"]], stn[["sigma0", "unesco"]])
    #(before issue1933) expect_equal(stn[["sigmaTheta", "gsw"]], stn[["sigma0", "gsw"]])
    expect_false(identical(stn[["sigmaTheta", "gsw"]], stn[["sigma0", "gsw"]]))

    # sigmaTheta and sigma0 ([[ form) should not match between unesco and gsw
    expect_false(identical(stn[["sigmaTheta", "gsw"]], stn[["sigmaTheta", "unesco"]]))
    expect_false(identical(stn[["sigma0", "gsw"]], stn[["sigma0", "unesco"]]))

    # sigmaTheta and sigma0 (function form) should not match between unesco and gsw
    expect_false(identical(swSigmaTheta(stn, eos="unesco"), swSigmaTheta(stn, eos="gsw")))
    expect_false(identical(swSigma0(stn, eos="unesco"), swSigma0(stn, eos="gsw")))
})

test_that("potential_temperature (UNESCO)", {
    # 2.1 UNESCO potential temperature
    #
    # The following is an official test value from [1 p44], first with all args,
    # second with a ctd object as an arg. Note the need to convert to the 1968
    # temperature scale, which was used in the UNESCO formulation.
    expect_equal(swTheta(40, T90fromT68(40), 10000, eos="unesco"),
        T90fromT68(36.89073), tolerance=0.00002)
    expect_equal(swTheta(as.ctd(40, T90fromT68(40), 10000), eos="unesco"),
        T90fromT68(36.89073), tolerance=0.00002)
    # Test self-consistency at the surface (also a test of vector reframing)
    tem <- 10 + rnorm(50)
    expect_equal(0, sum(abs(tem - swTheta(rep(35, 50), tem, 0, eos="unesco"))))
})

test_that("SA and CT, sound speed (GSW)", {
    # 2.2 GSW potential temperature
    #
    # Since gsw_ functions are tested in the gsw package, we just need a consistency check.
    SP <- 35
    t <- 13                                 # notation in gsw_...() functions
    p <- 1000
    lon <- 300
    lat <- 30
    ctd <- as.ctd(SP, t, p, longitude=lon, latitude=lat)
    SA <- gsw::gsw_SA_from_SP(SP, p, longitude=lon, latitude=lat)
    thetaGSW <- gsw::gsw_pt_from_t(SA, t, p, p_ref=0)
    theta <- swTheta(ctd, eos="gsw")
    expect_equal(thetaGSW, theta)
    CT <- gsw::gsw_CT_from_t(SA=SA, t=t, p=p)
    expect_equal(SA, swAbsoluteSalinity(salinity=SP, pressure=p, longitude=lon, latitude=lat))
    expect_equal(SA, swAbsoluteSalinity(ctd))
    expect_equal(CT, swConservativeTemperature(salinity=SP, temperature=t, pressure=p, longitude=lon, latitude=lat))
    expect_equal(CT, swConservativeTemperature(ctd))
    expect_equal(1731.995, swSoundSpeed(40, T90fromT68(40), 1e4, eos="unesco"), tolerance=0.001)
    expect_equal(1731.995, swSoundSpeed(as.ctd(40, T90fromT68(40), 1e4), eos="unesco"), tolerance=0.001)
    SA <- gsw::gsw_SA_from_SP(SP=40, p=1e4, longitude=300, latitude=30)
    CT <- gsw::gsw_CT_from_t(SA, 40, 1e4)
    speedGSW <- gsw::gsw_sound_speed(SA, CT, 1e4)
    speed <- swSoundSpeed(salinity=40, temperature=40, pressure=1e4, longitude=300, latitude=30, eos="gsw")
    expect_equal(speedGSW, speed)
})

test_that("temperature scales", {
    expect_equal(10, T68fromT90(T90fromT68(10)))
})

test_that("freezing temperature", {
    # 5.1 UNESCO freezing temperature [1 p29]
    Tf <- swTFreeze(40, 500, eos="unesco")
    expect_equal(Tf, T90fromT68(-2.588567), tolerance=1e-6)
    # 5.2 GSW freezing temperature. This is actually just a test that
    # the GSW functions are called correctly -- it is *not* using
    # check values. However, we don't have to worry about check
    # values, because in the gsw package test suite, there are multiple
    # tests that check values against those published on the GSW/TEOS10
    # website.
    SA <- gsw::gsw_SA_from_SP(SP=40, p=500, longitude=300, latitude=30)
    TfGSW <- gsw::gsw_t_freezing(SA=SA, p=500, saturation_fraction=1)
    Tf <- swTFreeze(40, 500, longitude=300, latitude=30, eos="gsw")
    expect_equal(TfGSW, Tf)
})

test_that("specific heat", {
    # [1 p31]
    p <- 1e4
    t <- 40
    SP <- 40
    lon <- 300
    lat <- 30
    C <- swSpecificHeat(salinity=SP, temperature=T90fromT68(t), pressure=p, eos="unesco")
    expect_equal(C, 3849.499, tolerance=1e-3)
    SA <- gsw::gsw_SA_from_SP(SP=SP, p=p, longitude=lon, latitude=lat)
    CGSW <- gsw::gsw_cp_t_exact(SA=SA, t=t, p=p)
    C <- swSpecificHeat(salinity=SP, temperature=t, pressure=p, longitude=lon, latitude=lat, eos="gsw")
    expect_equal(CGSW, C)
})

test_that("lapse rate", {
    # [1 p38]
    SP <- 40
    t <- 40
    p <- 1e4
    l <- swLapseRate(salinity=SP, temperature=T90fromT68(t), pressure=p, eos="unesco")
    expect_equal(l, 3.255976e-4, tolerance=1e-7)
    lon <- 300
    lat <- 30
    SA <- gsw::gsw_SA_from_SP(SP=SP, p=p, longitude=lon, latitude=lat)
    CT <- gsw::gsw_CT_from_t(SA=SA, t=t, p=p)
    lGSW <- 1e4*gsw::gsw_adiabatic_lapse_rate_from_CT(SA=SA, CT=CT, p=p) # convert to deg/m
    l <- swLapseRate(salinity=SP, temperature=t, pressure=p, longitude=lon, latitude=lat, eos="gsw")
    expect_equal(lGSW, l)
})

test_that("alpha and beta", {
    # 8 alpha and beta
    # 8.1 UNESCO alpha and beta
    # The formulae used are not actually from UNESCO, but are rather given in [2],
    # and the test values come from [2] also.
    #
    # Since [2] gives formula is in terms of theta=10C, we must compute the
    # corresponding in-situ temperature first. Use S=40 and 4000dbar to match
    # his check value.
    tem <- uniroot(function(x) 10-swTheta(40, x, 4000, eos="unesco"), c(9, 12))$root
    # The beta=7.2088e-4 value is from the last sentence of McDougall's Appendix.
    expect_equal(7.2088e-4, swBeta(40, tem, 4000, eos="unesco"), tolerance=4e-8)
    # The alpha/beta=0.34763 is from the left-hand column of McDougall's p1964.
    expect_equal(0.34763, swAlphaOverBeta(40, tem, 4000, eos="unesco"), tolerance=2e-5)
    expect_equal(0.34763*7.20883e-4, swAlpha(40, tem, 4000, eos="unesco"), tolerance=2e-5)
    # 8.1 GSW alpha and beta
    # Check against gsw_ values, which we know to be correct from the gsw test suite.
    SP <- 40
    t <- 10
    p <- 4000
    lon <- 300
    lat <- 30
    a <- swAlpha(SP, t, p, longitude=lon, latitude=lat, eos="gsw")
    b <- swBeta(SP, t, p, longitude=lon, latitude=lat, eos="gsw")
    SA <- gsw::gsw_SA_from_SP(SP=SP, p=p, longitude=lon, latitude=lat)
    CT <- gsw::gsw_CT_from_t(SA=SA, t=t, p=p)
    aGSW <- gsw::gsw_alpha(SA=SA, CT=CT, p=p)
    expect_equal(a, aGSW)
    bGSW <- gsw::gsw_beta(SA=SA, CT=CT, p=p)
    expect_equal(b, bGSW)
    # 8.2 swAlphaOverBeta
    # Ensure that alpha, beta, and alpha/beta are consistent, in both EOS
    sal <- 34
    tem <- 10
    pre <- 100
    expect_equal(swAlphaOverBeta(sal, tem, pre, longitude=300, latitude=30, eos="gsw"),
        swAlpha(sal, tem, pre, longitude=300, latitude=30, eos="gsw") /
        swBeta(sal, tem, pre, longitude=300, latitude=30, eos="gsw"))
    expect_equal(swAlphaOverBeta(sal, tem, pre, eos="unesco"),
        swAlpha(sal, tem, pre, eos="unesco") /
        swBeta(sal, tem, pre, eos="unesco"))
})

test_that("swSTrho", {
    # 9. swSTrho
    # This is used to draw isopycnals on TS diagrams.
    tem <- 10
    rho <- 1022
    pre <- 0
    # 9.1 UNESCO swSTrho
    Su <- swSTrho(T90fromT68(tem), rho, pre, eos="unesco")
    # Next was 28.65114808083 before issue 2044, but the precision
    # of the calculation was increased then, so a new check value
    # is needed here. The relative difference is 1.7e-7, so
    # certainly not a concern, but we want this test suite to
    # check on changes to the code, in addition to checking
    # on test values reflecting external knowledge.
    expect_equal(Su, 28.6511432379484)
    expect_equal(rho, swRho(Su, T90fromT68(tem), 0, eos="unesco"))
    # 9.2 GSW swSTrho
    CT <- gsw::gsw_CT_from_t(Su, tem, pre)
    Sg <- swSTrho(CT, rho, pre, eos="gsw")
    expect_equal(gsw::gsw_rho(Sg, CT, pre), rho)
})

test_that("misc sw calculations", {
    # The following was hard-coded using values from GSW3.03, and it failed with GSW3.05.
    # expect_equal(Sg, 28.7842812841013,tolerance=1e-8)
    tem <- swTSrho(35, 23, 0, eos="unesco") # 26.11301
    # As above, the check value was changed in addressing issue
    # 2044.  The old value was 26.1130113601685, which has
    # relative difference of 1.1e-6 compared to the new value. Note
    # that the value differs because tem differs, reflecting the
    # change in swTSrho(); there is no change to T68fromT90().
    expect_equal(T68fromT90(tem), 26.1130395531654, tolerance=1e-8)
    expect_equal(swRho(35, tem, 0, eos="unesco"), 1023, tolerance=1e-5)
})

test_that("sound absorption", {
    # Compared with Table IV of Fisher & Simmons 1977.
    alpha <- swSoundAbsorption(100e3, 35, 4, 4990) # at 500 atm (4990 dbar of water)
    expect_equal(alpha, 0.0175, tolerance=0.01) # 1% test
    alpha <- swSoundAbsorption(10e3, 35, 4, 0) # expect 0.00083 at 1 atm (0dbar of water)
    expect_equal(alpha, 0.000829, tolerance=0.01) # 1% test
})

test_that("viscosity", {
    # This is just a test against future changes, for
    # the original reference did not provide a test value.
    expect_equal(1000*swViscosity(30, 10), 1.383779, tolerance=0.000002)
})

test_that("thermal conductivity", {
    # Caldwell 1975 table 1 gives 4 digits, i.e. to 1e-6
    joulePerCalorie <- 4.18400
    cmPerM <- 100
    test <- swThermalConductivity(31.5, 10, 1000) / joulePerCalorie / cmPerM
    expect_equal(1000*test, 1.478, tolerance=0.001)
})

test_that("electrical conductivity: definitional check values", {
    expect_equal(swCSTp(35, T90fromT68(15), 0, eos="unesco"),  1)
    expect_equal(swCSTp(35, T90fromT68(15), 0, eos="gsw"), 1)
    expect_equal(swSCTp(1, T90fromT68(15), 0, eos="unesco"), 35)
    expect_equal(swSCTp(1, T90fromT68(15), 0, eos="gsw"), 35)
    expect_equal(swSCTp(0.5, 10, 100, eos="unesco"),
        swSCTp(0.5, 10, 100, eos="gsw"))
    # These test values are not against a known standard; rather, they simply
    # assure that there has been no change since a test done on 2019 Mar 23
    # whilst working on issue https://github.com/dankelley/oce/issues/1514
    expect_equal(swSCTp(0.02, 10, 100, eos="gsw"), 0.601398102117915)
    expect_equal(swSCTp(0.02, 10, 100, eos="unesco"), 0.601172086373874)
})

test_that("electrical conductivity: semi-definitional check values (AUTHOR IS CONFUSED ON THESE)", {
    # the C=1 value can be tested directly in gsw, but others are tested against gsw.
    SP <- swSCTp(1.2, T90fromT68(20), 2000, eos="gsw")
    expect_equal(1.2, gsw::gsw_C_from_SP(SP, T90fromT68(20), 2000) / gsw::gsw_C_from_SP(35, T90fromT68(15), 0))
    SP <- swSCTp(0.65, T90fromT68(5), 1500, eos="gsw")
    expect_equal(0.65, gsw::gsw_C_from_SP(SP, T90fromT68(5), 1500) / gsw::gsw_C_from_SP(35, T90fromT68(15), 0))
})

test_that("electrical conductivity: real-data checks", {
    data(ctd)
    # This does not have conductivity, so add it
    salinity <- ctd[["salinity"]]
    temperature <- ctd[["temperature"]]
    pressure <- ctd[["pressure"]]
    conductivity <- swCSTp(salinity, temperature, pressure, eos="unesco")
    ctd <- oceSetData(ctd, name="conductivity", value=conductivity, unit=list(unit=expression(), scale=""))
    S <- swSCTp(ctd)
    expect_equal(S, salinity, tolerance=1e-3)
    # Test that swCSTp() takes both salinity and CTD [issue 630]
    cond1 <- swCSTp(salinity, temperature, pressure, eos="unesco")
    cond2 <- swCSTp(ctd)
    expect_equal(cond1, cond2)
})

test_that("depth and pressure", {
    # 14. depth and pressure
    # The UNESCO test is basically for consistency with old versions, I think,
    # but the GSW test is against gsw_z_from_p(), which is well-tested in
    # the building of the gsw package.
    depth <- swDepth(10000, 30, eos="unesco")
    expect_equal(depth, 9712.653, tolerance=0.001)
    depth <- swDepth(10000, 30, eos="gsw")
    expect_equal(depth, 9713.735, tolerance=0.001)
    pressure <- swPressure(9712.653, 30, eos="unesco")
    expect_equal(pressure, 10000., tolerance=0.001)
    pressure <- swPressure(9712.653, 30, eos="gsw")
    expect_equal(pressure, gsw::gsw_p_from_z(-9712.653, 30), tolerance=0.001)
})

test_that("spiciness", {
    expect_error(swSpice(35, T90fromT68(10), 100, eos="gsw"), "must supply longitude")
    # Q: is this test value from Flament's paper, or is it just a consistency check?
    sp <- swSpice(35, T90fromT68(10), 100, eos="unesco")
    expect_equal(sp, 1.131195, tolerance=0.0000015)
    # compare against direct gsw:: computation
    data(ctd)
    sal <- ctd[["salinity"]]
    tem <- ctd[["temperature"]]
    pre <- ctd[["pressure"]]
    lon <- rep(ctd[["longitude"]], length(sal))
    lat <- rep(ctd[["latitude"]], length(sal))
    piOce <- swSpice(sal, tem, pre, longitude=lon, latitude=lat, eos="gsw")
    piGsw <- gsw::gsw_spiciness0(ctd[["SA"]], ctd[["CT"]])
    expect_equal(piOce, piGsw)
})

test_that("CTD object accessors for derived properties", {
    data(ctd)
    sigma0 <- swSigma0(ctd[["salinity"]], ctd[["temperature"]], ctd[["pressure"]],
        longitude=ctd[["longitude"]], latitude=ctd[["latitude"]])
    sigmaTheta <- swSigmaTheta(ctd[["salinity"]], ctd[["temperature"]], ctd[["pressure"]],
        longitude=ctd[["longitude"]], latitude=ctd[["latitude"]])
    spice <- swSpice(ctd[["salinity"]], ctd[["temperature"]], ctd[["pressure"]],
        longitude=ctd[["longitude"]], latitude=ctd[["latitude"]])
    sigmaTheta <- swSigmaTheta(ctd)
    expect_equal(sigma0, ctd[["sigma0"]])
    expect_equal(sigmaTheta, ctd[["sigmaTheta"]])
    expect_equal(spice, ctd[["spice"]])
    expect_equal(sigma0, swSigma0(ctd))
    expect_equal(sigmaTheta, swSigmaTheta(ctd))
    expect_equal(spice, swSpice(ctd))
})

test_that("non-CTD object accessors for derived properties", {
    data(ctd)
    sigma0 <- swSigma0(ctd)
    sigmaTheta <- swSigmaTheta(ctd)
    spice <- swSpice(ctd)
    general <- new("oce")
    general <- oceSetData(general, "temperature", ctd[["temperature"]])
    general <- oceSetData(general, "salinity", ctd[["salinity"]])
    general <- oceSetData(general, "pressure", ctd[["pressure"]])
    general <- oceSetMetadata(general, "longitude", ctd[["longitude"]])
    general <- oceSetMetadata(general, "latitude", ctd[["latitude"]])
    # Test whether [[ works on the three special-case derived quantities (see
    # R/sw.R near line 363 for how '[[' intercepts these three things at the
    # deep level for oce objects, as opposed to CTD objects).
    expect_equal(sigma0, general[["sigma0"]])
    expect_equal(sigma0, swSigma0(general))
    expect_equal(sigmaTheta, general[["sigmaTheta"]])
    expect_equal(sigmaTheta, swSigmaTheta(general))
    expect_equal(spice, general[["spice"]])
    expect_equal(spice, swSpice(general))
})

test_that("swRho handles matrix and array data", {
    data(ctd)
    lon <- ctd[["longitude"]]
    lat <- ctd[["latitude"]]
    sal <- ctd[["salinity"]][1:20]
    tem <- ctd[["temperature"]][1:20]
    pre <- ctd[["pressure"]][1:20]
    Sm <- matrix(sal, nrow=10)
    Tm <- matrix(tem, nrow=10)
    pm <- matrix(pre, nrow=10)
    Sa <- array(sal, dim=c(2, 2, 5))
    Ta <- array(tem, dim=c(2, 2, 5))
    pa <- array(pre, dim=c(2, 2, 5))
    # unesco equation of state
    rho <- swRho(sal, tem, pre, eos="unesco")
    expect_equal(swRho(Sm, Tm, pm, eos="unesco"),
        matrix(rho, nrow=10))
    expect_equal(swRho(Sa, Ta, pa, eos="unesco"),
        array(rho, dim=c(2, 2, 5)))
    # gsw equation of state
    rho <- swRho(sal, tem, pre, longitude=lon, latitude=lat, eos="gsw")
    expect_equal(swRho(Sm, Tm, pm, longitude=lon, latitude=lat, eos="gsw"),
        matrix(rho, nrow=10))
    expect_equal(swRho(Sa, Ta, pa, longitude=lon, latitude=lat, eos="gsw"),
        array(rho, dim=c(2, 2, 5)))
})

test_that("sigma0 works as expected (issue 1933)", {
    # https://github.com/dankelley/oce/issues/1933
    data(section)
    ctd <- section[["station", 100]]
    sigma0 <- ctd[["sigma0", "gsw"]]
    CT <- ctd[["CT"]]
    SA <- ctd[["SA"]]
    sigma0direct <- gsw_sigma0(SA=SA, CT=CT)
    expect_equal(sigma0, sigma0direct)
})