Add dex

bumbanian · bumbanian · commit e6d52351033c · 2022-04-06T16:10:54.000-06:00
diff --git a/CRAN-RELEASE b/CRAN-RELEASE
@@ -0,0 +1,2 @@
+This package was submitted to CRAN on 2021-06-09.
+Once it is accepted, delete this file and tag the release (commit c39aa1c).
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -26,3 +26,4 @@ License: GPL-3
 Language: en-US
 Encoding: UTF-8
 RoxygenNote: 7.1.1
+LazyData: true
diff --git a/NAMESPACE b/NAMESPACE
@@ -5,6 +5,7 @@ export(iso)
 export(wiDB_sites)
 export(wiDB_data)
 export(wiDB_values)
+export(dex)
 importFrom(httr, GET, content, user_agent)
 importFrom(jsonlite, fromJSON)
 importFrom(doParallel, registerDoParallel, stopImplicitCluster)
@@ -14,4 +15,4 @@ importFrom(R2WinBUGS, monitor)
 importFrom(stats, sd, var)
 importFrom(graphics, abline, lines, par, points)
 importFrom(abind, abind)
-importFrom(utils, read.csv, unzip)
+importFrom(utils, read.csv, unzip, data)
diff --git a/NEWS.md b/NEWS.md
@@ -1,6 +1,8 @@
 # isoWater news
 
 ## isoWater 1.0.1.9000
+* Add dex function
+* Bug and spelling fixes
 
 ## isoWater 1.0.1
 * Bug fixes
diff --git a/R/constructors.R b/R/constructors.R
@@ -18,7 +18,7 @@ iso = function(H, O, Hsd, Osd, HOc = 0){
 #####
 
 mwl = function(HO, plot = TRUE){
-  if(class(HO)[1] != "data.frame"){
+  if(!inherits(HO, "data.frame")){
     stop("HO must be a data.frame")
   }
   if(!is.numeric(HO[,1]) | !is.numeric(HO[,2])){
@@ -69,3 +69,41 @@ mwl = function(HO, plot = TRUE){
   class(mwl) = "mwl"
   return(mwl)
 }
+
+#####
+#--- deuterium excess ----
+#calculate dex for one or more samples
+#####
+
+dex = function(HO, form = "dex", MWL = NULL){
+  if(!inherits(HO, "data.frame")){
+    stop("HO must be a data.frame")
+  }
+  if(!is.numeric(HO[,1]) | !is.numeric(HO[,2])){
+    stop("Non-numeric values in HO")
+  }
+  
+  if(is.null(MWL)){
+    data("GMWL", envir = environment())
+    GMWL = GMWL
+    MWL = GMWL
+  }
+  
+  if(!is.numeric(MWL) | length(MWL) < 2){
+    stop("MWL must be mwl object or numeric vector with slope, intercept")
+  }
+  
+  if(!(form %in% c("dex", "lcex", "both"))){
+    stop("form must be dex, lcex, or both")
+  }
+
+  if(form %in% c("dex", "both")){
+    HO$dex = HO[,1] - 8 * HO[,2]
+  }
+  
+  if(form %in% c("lcex", "both")){
+    HO$lcex = HO[,1] - MWL[1] * HO[,2] - MWL[2]
+  }
+  
+  return(HO)
+}
diff --git a/R/watercomp.R b/R/watercomp.R
@@ -4,14 +4,19 @@
 
 #takes values of observed water (type 'iso'), MWL (see below), hypothesized EL slope value
 #and number of parameter draws
-mwlSource = function(obs, MWL=c(8.01, 9.57, -8.096, 2564532.2, 5.76, 80672), 
-                     slope, stype = 1, ngens=1e4, ncores = 1){
+mwlSource = function(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1){
 
+  if(is.null(MWL)){
+    data("GMWL", envir = environment())
+    GMWL = GMWL
+    MWL = GMWL
+  }
+  
   if(MWL[6] == 80672 & stype == 2){
     warning("Using stype=2 and GMWL is inappropriate for most applications; see man")
   }
   
-  if(length(class(obs)) < 2 | class(obs)[2] != "iso"){
+  if(!inherits(obs, "iso")){
     warning("Expecting iso object for obs, this argument may be formatted incorrectly")
   }
   if(!(stype %in% c(1, 2))){
@@ -101,7 +106,7 @@ mwlSource = function(obs, MWL=c(8.01, 9.57, -8.096, 2564532.2, 5.76, 80672),
 mixSource = function(obs, sources, slope, prior=rep(1,nrow(sources)), 
                    shp=1, ngens=1e5, ncores = 1){
 
-  if(length(class(obs)) < 2 | class(obs)[2] != "iso"){
+  if(!inherits(obs, "iso")){
     warning("Expecting iso object for obs, this argument may be
             formatted incorrectly")
   }
@@ -118,7 +123,7 @@ mixSource = function(obs, sources, slope, prior=rep(1,nrow(sources)),
     obs.vcov[2 + (i - 1) * 2, 2] = obs[i, 4] ^ 2
   }
   
-  if(length(class(sources)) < 2 | class(sources)[2] != "iso"){
+  if(!inherits(sources, "iso")){
     warning("Expecting iso object for sources, this argument may be
             formatted incorrectly")
   }
diff --git a/R/wiDBfunctions.R b/R/wiDBfunctions.R
@@ -7,45 +7,45 @@ wiDB_validate = function(minLat, maxLat, minLong, maxLong, minElev, maxElev,
   qStr = ""
   
   if(!is.null(minLat)){
-    if(class(minLat) != "numeric"){stop("minLat must be numeric")}
+    if(!inherits(minLat, "numeric")){stop("minLat must be numeric")}
     qStr = paste0(qStr, "&minLat=", minLat)
   }
   if(!is.null(maxLat)){
-    if(class(maxLat) != "numeric"){stop("maxLat must be numeric")}
+    if(!inherits(maxLat, "numeric")){stop("maxLat must be numeric")}
     qStr = paste0(qStr, "&maxLat=", maxLat)
   }
   if(!is.null(minLong)){
-    if(class(minLong) != "numeric"){stop("minLong must be numeric")}
+    if(!inherits(minLong, "numeric")){stop("minLong must be numeric")}
     qStr = paste0(qStr, "&minLong=", minLong)
   }
   if(!is.null(maxLong)){
-    if(class(maxLong) != "numeric"){stop("maxLong must be numeric")}
+    if(!inherits(maxLong, "numeric")){stop("maxLong must be numeric")}
     qStr = paste0(qStr, "&maxLong=", maxLong)
   }
   if(!is.null(minElev)){
-    if(class(minElev) != "numeric"){stop("minElev must be numeric")}
+    if(!inherits(minElev, "numeric")){stop("minElev must be numeric")}
     qStr = paste0(qStr, "&minElev=", minElev)
   } 
   if(!is.null(maxElev)){
-    if(class(maxElev) != "numeric"){stop("maxElev must be numeric")}
+    if(!inherits(maxElev, "numeric")){stop("maxElev must be numeric")}
     qStr = paste0(qStr, "&maxElev=", maxElev)
   }
   if(!is.null(minDate)){
-    if(class(minDate) != "character"){stop("minDate must be string")}
+    if(!inherits(minDate, "character")){stop("minDate must be string")}
     td = c(as.numeric(substr(minDate, 1, 4)), as.numeric(substr(minDate, 6, 7)),
            as.numeric(substr(minDate, 9, 10)))
     if(NA %in% td){stop("minDate format must be YYYY-MM-DD")}
     qStr = paste0(qStr, "&minDate=", minDate)
   }
   if(!is.null(maxDate)){
-    if(class(maxDate) != "character"){stop("maxDate must be string")}
+    if(!inherits(maxDate, "character")){stop("maxDate must be string")}
     td = c(as.numeric(substr(maxDate, 1, 4)), as.numeric(substr(maxDate, 6, 7)),
            as.numeric(substr(maxDate, 9, 10)))
     if(NA %in% td){stop("maxDate format must be YYYY-MM-DD")}
     qStr = paste0(qStr, "&maxDate=", maxDate)
   }
   if(!is.null(countries)){
-    if(class(countries) != "character"){stop("countries must be string")}
+    if(!inherits(countries, "character")){stop("countries must be string")}
     countries = gsub("  ", "", countries)
     countries = gsub(" ", "", countries)
     if(length(countries > 1)){
@@ -54,7 +54,7 @@ wiDB_validate = function(minLat, maxLat, minLong, maxLong, minElev, maxElev,
     qStr = paste0(qStr, "&countries=", countries)
   }
   if(!is.null(states)){
-    if(class(states) != "character"){stop("states must be string")}
+    if(!inherits(states, "character")){stop("states must be string")}
     states = gsub("  ", "", states)
     states = gsub(" ", "", states)
     if(length(states > 1)){
@@ -63,7 +63,7 @@ wiDB_validate = function(minLat, maxLat, minLong, maxLong, minElev, maxElev,
     qStr = paste0(qStr, "&states=", states)
   }
   if(!is.null(types)){
-    if(class(types) != "character"){stop("types must be string")}
+    if(!inherits(types, "character")){stop("types must be string")}
     types = gsub("  ", "", types)
     types = gsub(" ", "", types)
     if(length(types > 1)){
@@ -72,7 +72,7 @@ wiDB_validate = function(minLat, maxLat, minLong, maxLong, minElev, maxElev,
     qStr = paste0(qStr, "&types=", types)
   }
   if(!is.null(projects)){
-    if(class(projects) != "character"){stop("projects must be string")}
+    if(!inherits(projects, "character")){stop("projects must be string")}
     projects = gsub("  ", "", projects)
     projects = gsub(" ", "", projects)
     if(length(projects > 1)){
@@ -148,7 +148,7 @@ wiDB_data = function(minLat = NULL, maxLat = NULL, minLong = NULL, maxLong = NUL
     if(!dir.exists(tmpdir)){stop("Unable to create directory")}
   }
   
-  if(class(clean) != "logical"){stop("clean must be TRUE/FALSE")}
+  if(!inherits(clean, "logical")){stop("clean must be TRUE/FALSE")}
   
   flist = c("Site_ID", "Site_Name", "Latitude", "Longitude", "Elevation", 
             "Sample_ID", "Type", "Start_Date", "Start_Time_Zone", 
@@ -158,7 +158,7 @@ wiDB_data = function(minLat = NULL, maxLat = NULL, minLong = NULL, maxLong = NUL
             "Project_ID")
   
   if(!is.null(fields)){
-    if(class(fields) != "character"){stop("fields must be a string")}
+    if(!inherits(fields, "character")){stop("fields must be a string")}
     fields = gsub("  ", "", fields)
     fields = gsub(" ", "", fields)
     fels = strsplit(fields, ",")
diff --git a/README.md b/README.md
@@ -16,5 +16,5 @@ Putman, A. L., & Bowen, G. J. (2019) A global database of the stable isotopic ra
 
 <!-- badges: start -->
   [![Build status](https://github.com/SPATIAL-Lab/isoWater/actions/workflows/r.yml/badge.svg)](https://github.com/SPATIAL-Lab/isoWater/actions)
-  [![codecov](https://codecov.io/gh/SPATIAL-Lab/isoWater/branch/master/graph/badge.svg)](https://codecov.io/gh/SPATIAL-Lab/isoWater) 
+  [![codecov](https://codecov.io/gh/SPATIAL-Lab/isoWater/branch/master/graph/badge.svg)](https://app.codecov.io/gh/SPATIAL-Lab/isoWater) 
   <!-- badges: end -->
diff --git a/data-raw/prepData.R b/data-raw/prepData.R
@@ -53,3 +53,6 @@ mixModel = "model{
 }"
 
 use_data(mwlModel, mixModel, internal = TRUE, overwrite = TRUE)
+
+GMWL = c(8.01, 9.57, -8.096, 2564532.2, 5.76, 80672)
+use_data(GMWL, overwrite = TRUE)
diff --git a/data/GMWL.rda b/data/GMWL.rda
diff --git a/man/GMWL.Rd b/man/GMWL.Rd
@@ -0,0 +1,21 @@
+\name{GMWL}
+\alias{GMWL}
+\docType{data}
+\title{
+Global Meteoric Water Line
+}
+\description{
+Parameters for the Global Meteoric Water Line fit to a global precipitation compilation in Bowen, et al. (2019).
+}
+\usage{data("GMWL")}
+\format{
+  The format is:
+ num [1:6] slope, intercept, average d18O, sum of squares in d18O, root mean square error, number of samples
+}
+\source{
+Bowen et al. (2019) Isotopes in the water cycle: Regional- to global-Scale patterns and applications. \emph{Annual Review of Earth and Planetary Sciences} \strong{47} 453--479. \doi{10.1146/annurev-earth-053018-060220}.
+}
+\examples{
+data(GMWL)
+}
+\keyword{datasets}
diff --git a/man/dex.Rd b/man/dex.Rd
@@ -0,0 +1,44 @@
+\name{dex}
+
+\alias{dex}
+
+\title{
+Deuterium excess
+}
+
+\description{
+Calculates deuterium excess or line-conditioned excess.
+}
+
+\usage{
+dex(HO, form = "dex", MWL = NULL)
+}
+
+\arguments{
+  \item{HO}{
+  data.frame. Hydrogen (column 1) and oxygen (column 2) isotope values for 1 or more water samples.
+  }
+  \item{form}{character. Calculate deuterium excess (\dQuote{dex}), line-conditioned excess (\dQuote{lcex}), or \dQuote{both}.}
+  \item{MWL}{
+  numeric. Vector the first two elements of which contain the meteoric water line slope and intercept (e.g., as created by \code{\link{mwl}}). The default value (if \code{MWL = NULL}) reflects the Global Meteoric Water Line estimated from a global precipitation compilation in Bowen, et al. (2019). Ignored for \code{form = "dex"}.
+  }
+}
+
+\value{
+Returns a copy of \code{HO} with an added field(s) \dQuote{dex} and/or \dQuote{lcex} containing the calculated values. Deuterium excess is calculated following Dansgaard (1964) as: \emph{dex = \eqn{\delta}2H - 8 * \eqn{\delta}18O}, and lc-excess following Landwehr & Coplen (2006) as \emph{lcex = \eqn{\delta}2H - m * \eqn{\delta}18O - b}, where \emph{m} and \emph{b} are the slope and intercept of \code{MWL}, respectively.
+}
+
+\references{
+Bowen et al. (2019) Isotopes in the water cycle: Regional- to global-Scale patterns and applications. \emph{Annual Review of Earth and Planetary Sciences} \strong{47} 453--479. \doi{10.1146/annurev-earth-053018-060220}.
+
+Dansgaard (1964) Stable isotopes in precipitation. \emph{Tellus} \strong{16} 436--468. \doi{10.1111/j.2153-3490.1964.tb00181.x}.
+
+Landwehr & Coplen (2006) Line-conditioned excess: A new method for characterizing stable hydrogen and oxygen isotope ratios in hydrologic systems. In \emph{Isotopes in Environmental Studies}, International Atomic Energy Agency, 132--135. \url{http://www-pub.iaea.org/MTCD/publications/PDF/CSP_26_web.pdf}.
+}
+
+\examples{
+O = runif(10, -15, -2)
+H = O * 8 + 10 + rnorm(10, 0, 6)
+d = dex(data.frame(H, O), form = "both")
+print(d)
+}
diff --git a/man/iso.Rd b/man/iso.Rd
@@ -22,10 +22,10 @@ iso(H, O, Hsd, Osd, HOc = 0)
   numeric. Oxygen isotope value or vector of oxygen isotope values.
   }
   \item{Hsd}{
-  numeric. 1 standard devation uncertainty of \code{H} (value or vector of values).
+  numeric. 1 standard deviation uncertainty of \code{H} (value or vector of values).
   }
   \item{Osd}{
-  numeric. 1 standard devation uncertainty of \code{O} (value or vector of values).
+  numeric. 1 standard deviation uncertainty of \code{O} (value or vector of values).
   }
   \item{HOc}{
   numeric. Covariance of \code{H} and \code{O} uncertainties.
diff --git a/man/mwl.Rd b/man/mwl.Rd
@@ -24,9 +24,9 @@ mwl(HO, plot = TRUE)
 }
 
 \details{
-\code{mwl} will return an error if fewer than 3 sample values are provided and a warning if fewer than 10 samples are provided or if the correlation coefficient betwen H and O values is less than 0.7. Sample values should span a broad enough range of isotope values to strongly constrain the MWL.
+\code{mwl} will return an error if fewer than 3 sample values are provided and a warning if fewer than 10 samples are provided or if the correlation coefficient between H and O values is less than 0.7. Sample values should span a broad enough range of isotope values to strongly constrain the MWL.
 
-Model II (reduced major axis) regression is used to accomodate errors on both isotope values.
+Model II (reduced major axis) regression is used to accommodate errors on both isotope values.
 }
 
 \value{
diff --git a/man/mwlSource.Rd b/man/mwlSource.Rd
@@ -11,8 +11,7 @@ Given parameters describing a meteoric water line in H-O isotope space, generate
 }
 
 \usage{
-mwlSource(obs, MWL=c(8.01, 9.57, -8.096, 2564532.2, 5.76, 80672), 
-                     slope, stype = 1, ngens=1e4, ncores = 1)
+mwlSource(obs, MWL = NULL, slope, stype = 1, ngens=1e4, ncores = 1)
 }
 
 \arguments{
@@ -37,7 +36,7 @@ mwlSource(obs, MWL=c(8.01, 9.57, -8.096, 2564532.2, 5.76, 80672),
 }
 
 \details{
-The prior distribution of source values is constrained by \code{MWL}, which contains the parameters: slope, intercept, average d18O, sum of squares in d18O, root mean square error, and number of samples for an emperically-determined meteoric water line. This object can be created from a H and O isotope dataset using the fuction \code{\link{mwl}}. The default values reflect the Global Meteoric Water Line estimated from a global precipitation compilation in Bowen, et al. (2019). \code{stype} determines how the source uncertainty about the MWL is calculated; the default (1, confidence interval) is appropriate if the source is best represented as an integrated mixture of the samples defining the MWL, whereas option 2 (prediction interval) is appropriate if the source is best represented as a single sample.
+The prior distribution of source values is constrained by \code{MWL}, which contains the parameters: slope, intercept, average d18O, sum of squares in d18O, root mean square error, and number of samples for an empirically-determined meteoric water line. This object can be created from a H and O isotope dataset using the function \code{\link{mwl}}. The default value (if \code{MWL = NULL}) reflects the Global Meteoric Water Line estimated from a global precipitation compilation in Bowen, et al. (2019). \code{stype} determines how the source uncertainty about the MWL is calculated; the default (1, confidence interval) is appropriate if the source is best represented as an integrated mixture of the samples defining the MWL, whereas option 2 (prediction interval) is appropriate if the source is best represented as a single sample.
 
 If \code{ncores} = 1, three chains will be run on a single core. If \code{ncores} > 1, \code{ncores} chains will be run in parallel on \code{ncores} cores. 
 }
diff --git a/man/wiDB_data.Rd b/man/wiDB_data.Rd
@@ -49,7 +49,7 @@ wiDB_data(minLat = NULL, maxLat = NULL, minLong = NULL,
     character. Vector of one or more two-letter state or province codes for query.
   }
   \item{types}{
-    character. Vector of one or more sample types for query. See vocabulary in the \href{https://wateriso.utah.edu/waterisotopes/pages/spatial_db/contribute.html}{wiDB tempate}.
+    character. Vector of one or more sample types for query. See vocabulary in the \href{https://wateriso.utah.edu/waterisotopes/pages/spatial_db/contribute.html}{wiDB template}.
   }
   \item{projects}{
     character. Vector of one or more project codes for query.
diff --git a/man/wiDB_sites.Rd b/man/wiDB_sites.Rd
@@ -48,7 +48,7 @@ wiDB_sites(minLat = NULL, maxLat = NULL, minLong = NULL, maxLong = NULL,
     character. Vector of one or more two-letter state or province codes for query.
   }
   \item{types}{
-    character. Vector of one or more sample types for query. See vocabulary in the \href{https://wateriso.utah.edu/waterisotopes/pages/spatial_db/contribute.html}{wiDB tempate}.
+    character. Vector of one or more sample types for query. See vocabulary in the \href{https://wateriso.utah.edu/waterisotopes/pages/spatial_db/contribute.html}{wiDB template}.
   }
   \item{projects}{
     character. Vector of one or more project codes for query.
diff --git a/tests/testthat/test-mwlSource.R b/tests/testthat/test-mwlSource.R
@@ -1,12 +1,22 @@
 #Prep MWL
 O = runif(10, -15, -2)
 H = O * 8 + 10 + rnorm(10, 0, 6)
+HO = data.frame(H, O)
 
-MWL = mwl(data.frame(H, O), plot = FALSE)
+MWL = mwl(HO, plot = FALSE)
 test_that("mwl works", {
    expect_length(MWL, 6)
 })
 
+d = dex(HO, form = "both")
+test_that("dex works", {
+  expect_s3_class(d, "data.frame")
+  expect_equal(d[,3], H - 8 * O)
+  expect_error(dex(H))
+  expect_error(dex(HO, "lcex", 12))
+  expect_error(dex(HO, "jim"))
+})
+
 obs = iso(-60, -6, 0.5, 0.1, 0)
 slope = c(5, 0.3)
 
diff --git a/vignettes/isoWater.Rmd b/vignettes/isoWater.Rmd

Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,2 @@`
	`1`	`+This package was submitted to CRAN on 2021-06-09.`
	`2`	`+Once it is accepted, delete this file and tag the release (commit c39aa1c).`