jpgattuso/seacarb-git

Seawater carbonate chemistry with R

seacarb-git's Issues

Recent changes 3.1.2

Dear Jim,

I do not understand the logic of the changes made on function error. Why is the uncertainty of the boron concentration considered as a dissociation constant and include in
I suggest instead to consider eBt, on it own as is done with ePt and eSit

The function currently looks like this:
errors <-
function(flag, var1, var2, S=35, T=25, Patm=1, P=0, Pt=0, Sit=0,
evar1=0, evar2=0, eS=0.01, eT=0.01, ePt=0, eSit=0, epK=c(0.004, 0.015, 0.03, 0.01, 0.01, 0.02, 0.02, 0.01),
method="ga", r=0, runs=10000,
k1k2='x', kf='x', ks="d", pHscale="T", b="u74", gas="potential", warn="y")

Perhaps the following would be better:
errors <-
function(flag, var1, var2, S=35, T=25, Patm=1, P=0, Pt=0, Sit=0,
evar1=0, evar2=0, eS=0.01, eT=0.01, ePt=0, eSit=0, eBt= 0.01, epK=c(0.004, 0.015, 0.03, 0.01, 0.01, 0.02, 0.02),
method="ga", r=0, runs=10000,
k1k2='x', kf='x', ks="d", pHscale="T", b="u74", gas="potential", warn="y")

Am I missing something? Please use the new branch 3.1.2 to make the changes.

Best,
Jean-Pierre

Issue with function at

From: Jean-Pierre Gattuso [email protected]
To: Lennart Bach [email protected]
Subject: Re: Seacarb alkalinity function "at"
Date: Thu, 14 Jan 2021 18:52:16 +0100

Dear Lennart,

I also wish you a healthy and enjoyable year 2021!

I apologise for such a belated reply. I have been quite overwhelmed lately. Many thanks for your input which is much
appreciated. This function was added a long time ago and we have not used it since. My colleagues at the lab use a script which
reads the files saved by the titrator and process them without using seacarb. For point #1: this function is for titrations with
HCl unfortified with NaCl (hence molarity and molinity are essentially the same). It is a regrettable mistake that it is not
mentioned in the help file. In the help, “Total alkalinity is estimated using the non-linear least-square procedure described by
Dickson et al. (2007).” is correct as SOP 3b is not mentioned, just the data processing recommended by Dickson et al. is.

So, where do we go from here? Clearly, the at function cannot stay as is. I propose to pull it off seacarb for now. A revised at
function could be developed and I would be happy to also add an air-buoyancy correction function. I do not have the time to do
that on my own but it would be great if you would be willing to proceed. I will assist you to make sure that the package
compiles OK (CRAN is very finicky about that) and you would be the first author of the new functions.

Best regards,
Jean-Pierre

On 6 Jan 2021, at 0:20, Lennart Bach wrote:

I wish you all the best for the New Year, health and good spirit for you and your family.
 
I recently build my alkalinity infrastructure here in Hobart. For the calculation of TA I wanted to use the seacarb function
“at”. When evaluating the titration data I corresponded with Kai Schulz and Matthew Humphreys who cross-checked my titration
data with their TA calculation function for Matlab (Kai) and Python (Matthew). Based on that I noted some potential issues with
the seacarb function. I will list them below.
 
The value C as in seacarb is indicated to describe the “Normality” of the acid (i.e. mol/L). However, the applied Dickson
routine wants mol/kg for the acid. This leads to a ~50 µmol offset and it took my quite long to figure out why my seacarb
calculated values were always way off those calculated by the others. I think the description may need to be corrected to maybe
“molinity” so that it is clear it is mol/kg?
The “weight” as measured with the balance should bey air-buoyancy corrected (see SOP21 in Dickson et al., 2007) in order to get
to “mass”. Seacarb indicates weight is needed as an input, which is I think not 100% correct. If you want I can send you the
R-code for the air-buoyancy correction. This could perhaps be a helpful function in seacarb as well.
This third point is more something I did not understand. I played around with the Etris and pHTris values and they have a
strong influence on the outcome. In principle, the Dickson et al., 2003 method used here does not need these values to begin
with. (Kai and Matthew do not have that in their scripts). So it is unclear to me why it is in seacarb as this seems to differ
from the Dickson approach which is referenced by Seacarb.
This forth point is perhaps a recommendation. It would be great to have EMF_0 as an output value because monitoring this value
can help to monitor the quality of the electrode.

Issues with errors.R

@jamesorr

Jim,

I am working on the Point B manuscript and would like to include some information about uncertainties. Submission is planned at the end of August will it be possible to use the new error function?

In any case, I found issues working on the data set (link below).

• the errors function breaks when one uses ePk=NULL as you do by default. Perhaps you should use if(is.null(ePk))... to bring in the default values (which should be mentioned, with the source, in errors.Rd.
• I think the units are not mol/kg as indicated in the errors man file are they mol/kg
• In the attached, you will find a short R script for testing. It should be self-explanatory. You will just need to install the dplyr package if you do not use it.

We can speak on the phone if needed.

Archive.zip

Issue with Ks

I get the error below on glodap data centered on the Med and Red Seas. The code of Ks.R lines 92-98 is weird. I believe that you forgot to distinguish to take into consideration the formulations on Dickson and Khoo. Let me know.

Cheers,
Jean-Pierre

---

Error in if (any(is_w & (T > 45 | S > 45 | T < 0 | S < 5))) warning("S and/or T is outside the range of validity of the formulations available for Ks in seacarb.") :
missing value where TRUE/FALSE needed

buffesm does not work

See, for example, the first example mentioned in the man file:

buffesm(flag=8, var1=8.2, var2=0.00234, S=35, T=25, P=0, Pt=0, Sit=0, pHscale="T", kf="pf", k1k2="l", b="u74")

It returns NaNs

Avoiding wrong inputs

@jamesorr

I thought about your question about how we could make sure that users provide parameters in the right unit. We could make a test and provide a warning if, say:

• DIC is outside the range 140e-6 to 3000e-6
• TA is outside the range 600e-6 to 3000e-6
• Sit is outside the range 0 to 300e-6
• Pt is outside the range 0 to 5e-6

These ranges are estimates based on GLODAP v2.

What do you think?

Jean-Pierre

Error in carb package

The carb package releases an error <argument "var1" is missing, with no default>.

This error appears when I use the example df (seacarb_test_P300, seacarb_test_P0, or the x1 x2 df), or my own dataframe that is formatted correctly according to the guidelines for this package. I tried with flag formatted as a number and factor.

str(dframe)
'data.frame': 225700 obs. of 9 variables:
$ flag: num 24 24 24 24 24 24 24 24 24 24 ...
$ var1: num 500 496 489 491 496 ...
$ var2: num 0.00229 0.00229 0.00229 0.00229 0.00229 ...
$ S : num 33.4 33.4 33.4 33.4 33.4 ...
$ T : num 17.9 17.9 18 18 18.1 ...
$ P : num 0 0 0 0 0 0 0 0 0 0 ...
$ Pt : num 0 0 0 0 0 0 0 0 0 0 ...
$ Sit : num 0 0 0 0 0 0 0 0 0 0 ...
$ k1k2: Factor w/ 1 level "m10": 1 1 1 1 1 1 1 1 1 1 ...

carb(dframe)
Error in carb(dframe) : argument "var1" is missing, with no default

Avoiding inputs of incorrect lengths

In going through the function code, I became curious about the way that seacarb handles input variables of varying lengths. Currently, all functions look for the input variable with the maximum length, then force all shorter variables to that length by repeating the first value of that variable. This does deal with the majority of situations, where input values should be either the same length, or 1. However, if for some reason the user passes two inputs of different lengths (e.g. 4 and 5), the current way of dealing with input vector lengths would not throw an error, nor would it return 4 values. It would instead use the first value in the input of length 4 for all five calculations.

The matlab CO2SYS does things a similar way, but has another measure to stop this from occurring; that is, they make sure that there are no more than two unique vector lengths being inputted. The idea with this is to allow for constant inputs (e.g. all calculations with constant salinity), with all other values being equal. However, it would not throw an error for the first case presented (inputs of length 4 and 5).

The solution here is possibly to implement a line of code such that inputs can only be of length 1, or of the length of the maximum input. One way to do this may be:

"Kf" <-
function(S=35,T=25,P=0,kf='x',pHscale="T",Ks_p0=0,Ks_p=0,warn="y"){
    # create vector of input lengths
    input_lengths = c(length(S), length(T), length(P), length(kf), length(pHscale), length(Ks_p0), length(Ks_p))

    # create nK, or max vector length
    nK = max(input_lengths)

    # stop function if any input length is not 1 or the max value in input_lengths. 
    # this error text is copied from CO2SYS. 
    if (sum(! input_lengths %in% c(1, nK))>0){
    stop('INPUT ERROR: Input vectors must all be of same length, or of length 1. ')    
    }
    
    if(length(S)!=nK){S <- rep(S[1], nK)}
    if(length(T)!=nK){T <- rep(T[1], nK)}
    if(length(P)!=nK){P <- rep(P[1], nK)}
    if(length(kf)!=nK){kf <- rep(kf[1], nK)}
    if(length(pHscale)!=nK){pHscale <- rep(pHscale[1], nK)}
 
    ### continue function as normal. 

This should only allow the rest of the function to run if the correct assumption is met: that all inputs have the same length, unless they are of length one.

If this is of interest, I would be happy to fork and work on it. I thought I would ask before putting the work in.

Issue raised by CRAN

Dear maintainers,

This concerns the CRAN packages

[includes seacarb]

maintained by one of you:

[includes Jean-Pierre Gattuso]

R currently only emits a warning when if/while statement is used with a condition of length greater than one, e.g.

if (c(1,2)>0) TRUE [1] TRUE Warning message: In if (c(1, 2) > 0) TRUE : the condition has length > 1 and only the first element will be used

However, such cases are almost always programming errors, and hence this warning will soon be turned into a runtime error. This can already be done in current R-devel by setting environment variable R_CHECK_LENGTH_1_CONDITION=true to make it easier to spot these cases. And this is already done in certain CRAN checks.

A condition of length greater than one has been detected while running tests of your package (and the if/while statement involved was in the code of your package, so it is unlikely that this would have been caused by a dependent package). Detailed reports can be found at

https://github.com/kalibera/rifcond

I’ve modified R slightly to give more detailed debugging information in case of this error, so hopefully the outputs will be useful in finding the cause of the problem. In either case, one can reproduce the problem in normal R-devel via R_CHECK_LENGTH_1_CONDITION=true or in a released version of R by looking at warnings.

Please check out the reports and fix this in your package(s) as soon as possible. In the unlikely case this is caused by a dependent package, please report to the maintainer of that package.

Thank you, Tomas