R/waveStatsZC.R

Defines functions almostZero waveStatsZC

Documented in almostZero waveStatsZC

# Filename: waveStatsZC.R
# Based on the function zero_crossing.m by Urs Neumeier, v 1.06
# From http://neumeier.perso.ch/matlab/waves.html
# 
# 
# Author: Luke Miller  Mar 22, 2017
###############################################################################

#' Calculate wave statistics using zero-crossing method
#' 
#' Calculate ocean wave summary statistics, including significant
#' wave height and period. 
#' 
#' Based on an upward zero-crossing algorithm originally
#' provided by Urs Neumeier, v1.06. 
#' 
#' @param data A numeric vector of water surface height data. The data do not
#' need to be detrended prior to use. Typical units = meters
#' @param Fs Sampling frequency of the data, in Hz. 
#' @param threshold The minimum height necessary for a zero-crossing event to 
#' be considered a wave. 
#' @param plot Set to TRUE if summary histograms of wave heights and wave
#' periods are desired.
#' @return A list object containing summary statistic values.
#' 
#' \itemize{
#'   \item \code{Hsig} Mean of the highest 1/3 of waves in the data set. Units = 
#'   same as input surface heights.
#' 
#'   \item \code{Hmean}  Overall mean wave height, for all waves (bigger than 
#'   threshold). 
#' 
#'   \item \code{H10} Mean height of the upper 10\% of all waves. 
#' 
#'   \item \code{Hmax} Maximum wave height in the input data.
#' 
#'   \item \code{Tmean} Mean period of all waves (bigger than threshold). 
#'   Units = seconds.
#' 
#'   \item \code{Tsig} Mean period of \code{Hsig} (highest 1/3 of waves).
#' }
#' @references Original MATLAB function by Urs Neumeier:  
#' http://neumeier.perso.ch/matlab/waves.html
#' @seealso \code{\link{waveStatsSP}} for wave statistics determined using 
#' spectral analysis methods.
#' @examples
#' data(wavedata)
#' waveStatsZC(data = wavedata$SurfaceHeight.m, Fs = 4, plot = TRUE)
#' @export

waveStatsZC <- function(data, Fs, threshold = NULL, plot = FALSE){
	# Sanity checks
	if (!is.vector(data)){
		stop('data must be a vector of wave heights')
	}
	if (Fs < 0){
		stop('Frequency Fs must be greater than 0')
	}
	
#% The function was written for zero upward-crossing. 
#% To have zero downward-crossing (recommended) leave the next line uncommented
	data <- -data
	
	# detrendHeight function removes a linear trend and should result in a 
	# zero-mean timeseries
	detrended <- oceanwaves::detrendHeight(data) 
	data = detrended[['pt']][] # extract vector of detrended heights 
	
	# Make a shorter version of data with 0's removed
	d0 <- data[!almostZero(data,0)]	
	# Make a vector of indices the same length as 'data'
	back0 <- 1:length(data)
	# Keep only indices that are non-zero in 'data'
	back0 <- back0[!almostZero(data,0)]
 		
	# The code below multiplies each entry in d0 by the next entry in d0
	# and checks if the resulting value is negative (less than 0).
	# If two values are both positive heights, they will 
	# return a positive value, and if two values are both negative 
	# heights, then they will also return a positive value. But if a 
	# set of two values switches from positive to negative, or negative
	# to positive, the result will be a negative number (less than 0), and 
	# this must represent a crossing of the zero height level. 
	# The which() function then returns the indices where this is TRUE.
	f <- which(d0[1:(length(d0)-1)] * d0[2:length(d0)] < 0)
	# Extract the indices in 'data' where the height transitioned across 0
	crossing <- back0[f]
	
	# If the height record started off positive, then the first crossing must
	# be downward, and we will ignore it.
	if (data[1] > 0) {
		crossing <- crossing[-1] # Remove the first crossing index
	}

	# Subset the crossing vector (containing indices) in steps of 2. The 
	# first entry should be a zero upward crossing, since we just removed the 
	# former first entry if it was a downward crossing. Every 2nd entry after 
	# that should also represent an upward zero crossing event. 
	crossing <- crossing[seq(1, length(crossing), by = 2)]

	# Make a 4-column matrix to hold the height, crest, trough, and period data 
	# for each wave.
	wave <- matrix(NA, nrow = length(crossing)-1, ncol = 4)
	
	for (i in 1:(length(crossing)-1)){
		# Get the highest surface height for the interval between each 
		# upward zero crossing
		wave[i,2] <- max(data[crossing[i]:crossing[(i+1)]])
		# Get the lowest surface height for the interval between each 
		# upward zero crossing. This will be made positive for ease of use later
		wave[i,3] <- -min(data[crossing[i]:crossing[(i+1)]])
	}

	# If no wave was found in data, then do nothing. Otherwise, process each
	# wave
	if ( nrow(wave) > 0 ){
		# Calculate the time between each upward zero crossing by calculating
		# the number of rows between each upward crossing and dividing
		# by the sampling rate (Fs, units of Hz). This will be 
		# the period for each wave.
		wave[,4] <- diff(crossing) / Fs		
  
		if (is.null(threshold)){
			# If no minimum threshold for wave height was supplied, calculate a
			# threshold at 1% of the highest wave height (peak to trough) in the 
			# data set.
			threshold <- 0.01 * max(wave[,2] + wave[,3])
		} else if (threshold < 0){
			# If the supplied threshold was less than 0, throw an error
			stop('threshold must be greater than 0')
		}
		
		# Iterate through the 'wave' matrix and remove waves that are too small
		# and join them to the prior wave in the series
		 i <- 0
		 while (i < nrow(wave)){
			 i <- i + 1
			 # First check if the crest of wave i is less than threshold height
			# If the crest is too small, this will join the wave data with the
			# previous wave's data
			 if (wave[i,2] < threshold){
				 if (i != 1){
					# If we are on wave 2 or greater, rewrite wave i-1
					# with the max crest of wave i-1 or i, and the
					# max trough of wave i-i or i, and the
					# duration of wave i-1 and i.
					 wave[i-1, 2:4] <- c(max(wave[(i-1):i, 2]),
							 max(wave[(i-1):i, 3]),
							 sum(wave[(i-1):i, 4]))
				 }
				 # Remove this row that was less than the threshold
				 wave <- wave[-i, ]
				 
			# Next check if the trough was smaller than the threshold
			# If the trough is small, then this will join the wave data with
			# the next wave in the series
			 } else if (wave[i, 3] < threshold){
				 if (i != nrow(wave)){
					 # If we are on wave 2 or greater, rewrite wave i
					 # with the max crest of wave i or i+1, and the
					 # max trough of wave i or i+1, and the
					 # duration of wave i and i+1.
					 wave[i, 2:4] <- c(max(wave[i:(i+1), 2]),
							 max(wave[i:(i+1), 3]),
							 sum(wave[i:(i+1), 4]))
					 # Remove the data for the next wave in the series
					 wave <- wave[-(i+1), ]
				 } else {
					 wave <- wave[-i, ]
				 }
			 }
		 }
		 # Add the crest and trough heights together for each wave to get
		 # overall wave height
		 wave[,1] <- wave[, 2] + wave[, 3]
		
	} 	
					
# 'wave' has 1: wave height, 2: wave crest (Hcm), 3: wave trough (Htm), 4: period.

# Calculation of the wave statistics
	# Get the number of waves available to analyze
	nb <- nrow(wave) 
	# Sort the wave matrix based on the 1st column (wave height), large to small
	wave <- wave[order(wave[,1], decreasing = TRUE),]
	# Calculate the mean of the highest 1/3 of waves in the data set
	Hsig <- mean(wave[1:round(nb*(1/3)), 1])
	# Overall mean wave height, for all waves (bigger than threshold)
	Hmean <- mean(wave[ , 1])
	# Mean height of the upper 10% of all waves
	H10 <- mean(wave[1:round(nb*0.1), 1])
	# Maximum wave height
	Hmax <- max(wave[, 1])
	# Mean period
	Tmean <- mean(wave[, 4])
	# Mean period of Hsig (highest 1/3 of waves)
	Tsig <- mean(wave[1:round(nb*(1/3)), 4])
	
	# Assemble a list to return
	waveStats <- list(Hsig = Hsig, Hmean = Hmean, H10 = H10, Hmax = Hmax,
			Tmean = Tmean, Tsig = Tsig)
	
	if (plot){
		par(mfrow = c(2,1))
		hist(wave[,1], xlab = 'Wave height, m', main = '')
		hist(wave[,4], xlab = 'Wave Period, s', main = '')
		
	}

	# Return a list object with wave statistics
	waveStats

} # end of wave_stats_zc function



################################################################################ 
#' Test whether vector elements are effectively zero
#' 
#' Test whether elements in vector x are effectively zero,
#' within the square root of machine tolerance for a floating point number.
#' 
#' Returns TRUE for all vector elements that are within rounding error of zero.
#' 
#' @param x Numeric vector of values to compared against 0
#' @param y Value to be compared against. Default is 0. 
#' @param tolerance The maximum difference between x and y necessary to call a 
#' value non-zero
#' @return TRUE or FALSE for each element of the vector x. 

almostZero <- function(x, y = 0, tolerance=sqrt(.Machine$double.eps)) {
  diff <- abs(x - y)
  mag <- pmax(abs(x), abs(y))
  ifelse(mag > tolerance, diff/mag <= tolerance, diff <= tolerance)
}
millerlp/oceanwaves documentation built on Sept. 16, 2017, 1:24 a.m.