R/paleo_depth.r

In ericbarefoot/barefootr: A Package of Things I Find Useful

#	functions to calculate paleoflow depth from preserved cross-sets.
#	Eric Barefoot
#	Nov 2017
#
#	largely taken from work by Suzanne Leclair, Chris Paola and methods in a review paper by Bradley and Venditti

#	first step is to relate the cross-set (xset) thickness to dune height (H)

#	require packages

# require(stats4)

#' Dune height calculation from cross-set thicknesses
#'
#' \code{xset2H} takes measured cross-set thicknesses (\code{y}) and converts it to dune heights (H).
#'
#' @param y 	A vector of measured cross-set thicknesses in meters.
#' @param mode	Takes two values: \code{simple} and \code{pdfFit}. \code{simple} uses a simple approximation to determine the beta value in Leclair's equation relating x-sets to dune heights. \code{pdfFit} uses a more complicated procedure, fitting a PDF from Paola  + Borgman to get the beta value. \code{pdfFit} is better suited when there are adequate data to estimate the distribution.
#' @param output 	Takes two values: \code{vanilla} and \code{rich}. \code{vanilla} only returns the estimated value of \code{H}, or dune height. \code{rich} returns \code{H}, as well as \code{beta} and either \code{startVal} which is the inital guess of a for the MLE algorithm in \code{pdfFit}, or a span of \code{H} values giving a reasonable range of one particular assumption.
#' @return Either a single numeric if \code{output = vanilla}, or a list if \code{output = rich}.
#' @export

xset2H = function(y, mode = 'simple', output = 'vanilla') {

	#	function takes two modes: 'simple' and 'pdfFit'
	#	y is a vector of cross-set thicknesses -- MUST BE IN METERS

	y = y * 1000 # 	convert to millimeters

	#	gives output in meters

	n = length(y)

	#	first a function to estimate the average dune height via fitting a PDF of x-set thicknesses to find the right parameters
	pdfFit = function(y) {

		#	define the probability density function from Paola and Borgman
		pfun = function(a,s) {
			p = a * exp(-a*s) * (exp(-a*s) + a * s - 1) / (1 - exp(-a*s))^2
			return(p)
		}

		#	construct a likelihood function which takes data as an external argument and the parameter as the internal argument.
		loglikely = function(data) {
			function(a) {
				R = R = pfun(a, data)
				-sum(log(R))
			}
		}

		#	construct likelihood function for these data.
		LL = loglikely(y)

		#	default starting a value for fitting
		#	and error handling sequence for finding an inital guess value for that works and doesnt crash the MLE.
		worked = FALSE
		startVal = 1
		try(stop(''), silent = T)
		j = 1
		while (!worked) {

			try_MLE = try(coef(mle(minuslogl = LL, start = list(a = startVal))), silent = T)

			if (class(try_MLE) == 'try-error' & j < 1000) {

				options(warn = -1)
				msg = as.character(try_MLE)

				if (grepl('initial value', msg)) {
					startVal = startVal / 2
					j = j + 1
#					message(paste(startVal))
				} else {
					stop('Something else happened, its real bad')
				}

			} else {

				worked = TRUE

			}
		}

		#	estimate parameter using MLE function
		a = coef(mle(minuslogl = LL, start = list(a = startVal)))

		#	get beta a la Leclair et al. 0.9 is for dune height/scour correction
		beta = 0.9/a

		#	get dune height from beta from Leclair et al. two ways.
#		H = 2.22 * beta ^ 1.32
		H = 5.3 * beta + 0.001 * beta ^ 2

		list(H = as.numeric(H) / 1000, beta = as.numeric(beta), startVal = startVal)
	}

	#	Here, a simpler function using only the empirical estimations made by Leclair et al.
	simple = function(y) {
		sm = mean(y, na.rm = T) #	find average set thickness
		hsd_tssd = seq(from = 0.7, to = 1.1, by = 0.1) #	find ratio of scour deviation to height deviation. supposedly an average of 0.9, but ranges randomly from 0.7 to 1.1 this function spits out a range.
		beta = sm / (1.64493 / hsd_tssd) #	use that to calculate beta
#		beta = sm / 1.8

#		H = 2.22 * beta ^ 1.32
		H = 5.3 * beta + 0.001 * beta ^ 2 #	use same equation as in pdfFit to get the average height

		list(H = mean(H) / 1000, beta = as.numeric(beta), vecH = H / 1000)
	}

	#	offer some recommendations to the user.
	if(n >= 25 & mode == 'simple') {

		message('this dataset may benefit from running a more complicated PDF based function.')

		H = simple(y)

	} else if (n < 25 & mode == 'simple' ) {

		H = simple(y)

	} else if (n < 25 & mode == 'pdfFit') {

		message('this dataset is small, and may benefit from running a less complicated empirical.')

		H = pdfFit(y)

	} else if (n >= 25 & mode == 'pdfFit') {

		H = pdfFit(y)

	}

	switch(output ,
		   vanilla = return(H$H),
		   rich = return(H)
		  )
}

#' Paleo-depth calculation from dune height
#'
#' \code{H2h} takes estimated dune heights (H) and calculates a flow depth.
#'
#' @param H 	A vector of measured cross-set thicknesses in meters.
#' @param CL	A pre-specified confidence level. If you pick a low confidence (50\%), the range in estimated flow depth will be smaller.
#' @param output 	Takes two values: \code{vanilla} and \code{rich}. \code{vanilla} only returns the estimated value of \code{h}. \code{rich} returns \code{h}, as well as \code{CI}, the confidence interval around the \code{h} estimate.
#' @return Either a single numeric if \code{output = vanilla}, or a list if \code{output = rich}.
#' @export

#	second step is to relate the dune height to the flow depth, following Bradley and Venditti

H2h = function(H, CL = 50, output = 'vanilla') {

	#	H is the height of a dune in meters.
	#	or if H is a vector, it is a set of dune heights
	#	CL is the confidence level in percent

	#	confidence ranges

	CR = data.frame(CL = c(50, 80, 90, 95), upper = c(10.1, 14.6, 23.9, 29.5), lower = c(4.4, 3.1, 2.8, 2.7))

	Cind = which.min(abs(CR$CL - CL)) # match closest one from the table

	h = 6.7 * H	#	using the median parameter from Bradley and Venditti

	CI = list(upper = H * CR$upper[Cind], lower = H * CR$lower[Cind]) # and based on Confidence Level, assign bounds of confidence

	returnval = switch(output ,
					   vanilla = return(h),
					   rich = return(list(h = h, CI = CI))
					  )
}

#	and now a function that combines both into one function, so you can go straight from xsets to flow depth with no extra frills, or depth

#' Paleo-depth calculation from cross-set thicknesses
#'
#' \code{xset2depthSimple} uses the simplest implementation of both \code{xset2H} and \code{H2h}. It uses the vanilla options of each one and spits out an estimate of paleo-depth based on the defaults. Look in the code for the other functions to see what those defaults are.
#' @param y 	A vector of measured cross-set thicknesses in meters.
#' @return The paleo-depth estimate in meters.
#' @export

xset2depthSimple = function(y) {

	#	y is a vector of thicknesses of xsets in meters

	H = xset2H(y)

	h = H2h(H)

	return(h)

}