diff --git a/flow/tool/flow.shdw.corr.R b/flow/tool/flow.shdw.corr.R new file mode 100755 index 00000000..85c80cf9 --- /dev/null +++ b/flow/tool/flow.shdw.corr.R @@ -0,0 +1,84 @@ +##### +# Written by John Frank, USDA Forest Service, 240 W. Prospect Rd., Fort Collins, CO 80526, 970-498-1319, john.frank@usda.gov +# This code reads in a 10/20 Hz time series file and applies both the Wyngaard and Zhang (1985) shadowing correction and the "Self/Cross/Top" shadowing correction. +# The Wyngaard and Zhang correction is applied to the CSAT3 as described by Horst, Semmer, and Maclean (2015) +# The "Self/Cross/Top" shadowing correction was created by John Frank (2025) by reevaluating the data and analysis of Frank, Massman, and Ewers (2016) to account for CSAT3 transducer self-shadowing, transducer cross-shadowing, and top/bottom shadowing (the original formulation also included supporting arm shadowing, but it is omitted here due to its lack of significance). + +#------------------------------------------------------------------------------ +rm(list=ls()) #clear all +library(abind) +options(scipen = 8) + +##### +# Define matrices to convert transducer <-> sonic coordinates for CSAT3 +H_StoT = matrix(c(1/4,0.433012701892219,0.866025403784439,-1/2,0,0.866025403784439,1/4,-0.433012701892219,0.866025403784439),nrow=3,ncol=3) +G_TtoS = matrix(c(2/3,-4/3,2/3,1.154700538379252,0,-1.154700538379252,0.384900179459751,0.384900179459751,0.384900179459751),nrow=3,ncol=3) + +V_top = c(0,0,1) # unit vector for the top of the CSAT3 + +# Model parameters for the self/cross/top shadowing correction +d0 = c(0.00159718301425092,0.0113301152734556,0.0633033802230689) +d1 = c(29.2463049972615,24.8602835257800,18.5970206839229) +d2 = c(0.836978194375597,0.874138513424717,0.809209610043660) +d3 = 1.06606236917833 + +##### +# Read in time series data from each 30-minute data file +fname_in = sprintf('g:\\Studies\\Manuscripts\\2015\\Sonic Shadowing\\Analysis\\Beltsville, MD, OPE3\\Time series data\\TOA5_1580.ts_corn_2014_07_16_0000.dat') +TSData <- list() +TSData[[1]] <- as.data.frame(read.table(fname_in,sep = ",",skip=4,na.strings = c(".","NAN"))) +N_TSData = nrow(TSData[[1]]) # Number of rows of time series data in each file + +u <- abind(lapply(TSData,"[",3),along=1) +v <- abind(lapply(TSData,"[",4),along=1) +w <- abind(lapply(TSData,"[",5),along=1) + +Us = matrix(c(u,v,w),nrow = length(u)) # Matrix of time series data in sonic coordinates +Ut = Us %*% H_StoT # Matrix of time series data in transducer coordinates + +##### +# Calculate the angle (theta) between the 3-D wind and the 3 transducers plus the top/bottom (i.e. support structure) +Us_uv = sqrt(Us[,1]^2 + Us[,2]^2) +Us_uvw = sqrt(Us[,1]^2 + Us[,2]^2 + Us[,3]^2) + +theta_transducer = matrix(0,N_TSData,3) +theta_top = matrix(0,N_TSData,1) +for (i in 1:N_TSData) { + theta_transducer[i,1] = acos(abs(( H_StoT[,1] %*% Us[i,]) / Us_uvw[i]))*180/pi + theta_transducer[i,2] = acos(abs(( H_StoT[,2] %*% Us[i,]) / Us_uvw[i]))*180/pi + theta_transducer[i,3] = acos(abs(( H_StoT[,3] %*% Us[i,]) / Us_uvw[i]))*180/pi + theta_top[i,1] = acos(abs(( V_top %*% Us[i,]) / Us_uvw[i]))*180/pi +} + +##### +# Apply the Wyngaard and Zhang sinusoidal correction +c_Wyngaard = 1./(0.84 + 0.16*sin(theta_transducer*pi/180)) +Ut_Wyngaard = Ut*c_Wyngaard +Us_Wyngaard = Ut_Wyngaard %*% G_TtoS # Matrix of time series data in sonic coordinates with Wyngaard and Zhang sinusoidal correction + +##### +# Apply the self/cross/top shadowing correction + +# define self/cross/top shadowing correction function +f_c_self_cross_top = function(d0,d1,d2,d3){ + err_transducer_self = (1-d2[1])*tanh(d0[1]*(theta_transducer-d1[1]))/tanh(d0[1]*(90-d1[1]))+d2[1] + c_transducer_self = 1/err_transducer_self + + err_transducer_cross = (1-d2[2])*tanh(d0[2]*(theta_transducer-d1[2]))/tanh(d0[2]*(90-d1[2]))+d2[2] + c_transducer_cross = 1/err_transducer_cross + + err_top = (1-d2[3])*tanh(d0[3]*(theta_top-d1[3]))/tanh(d0[3]*(90-d1[3]))+d2[3] + c_top = 1/err_top + + c_self_cross_top = matrix(0,N_TSData,3) + c_self_cross_top[,1] = ((c_transducer_self[,1] - 1)^8 + (c_transducer_cross[,2] - 1)^8 + (c_transducer_cross[,3] - 1)^8 + (c_top - 1)^8)^(1/8) + 1 + c_self_cross_top[,2] = ((c_transducer_cross[,1] - 1)^8 + (c_transducer_self[,2] - 1)^8 + (c_transducer_cross[,3] - 1)^8 + (c_top - 1)^8)^(1/8) + 1 + c_self_cross_top[,3] = ((c_transducer_cross[,1] - 1)^8 + (c_transducer_cross[,2] - 1)^8 + (c_transducer_self[,3] - 1)^8 + (c_top - 1)^8)^(1/8) + 1 + c_self_cross_top = c_self_cross_top/d3 # scale this correction to minimize change in horizontal wind relative to the Wyngaard correction + + return(c_self_cross_top) +} + +c_self_cross_top = f_c_self_cross_top(d0,d1,d2,d3) +Ut_self_cross_top = Ut*c_self_cross_top +Us_self_cross_top = Ut_self_cross_top %*% G_TtoS # Matrix of time series data in sonic coordinates with self/cross/top shadowing correction diff --git a/pack/eddy4R.turb/R/def.corr.sct.csat.R b/pack/eddy4R.turb/R/def.corr.sct.csat.R new file mode 100644 index 00000000..dcccb9f1 --- /dev/null +++ b/pack/eddy4R.turb/R/def.corr.sct.csat.R @@ -0,0 +1,197 @@ + +#------------------------------------------------------------------------------- +# def.corr.sct.csat +# +# Apply Self/Cross/Top (SCT) shadowing correction (Frank, 2025) for CSAT3. +# Converts sonic-frame wind to transducer frame, computes SCT coefficients, +# applies correction element-wise, and maps back to sonic coordinates. +#------------------------------------------------------------------------------- + +#' @title Definition function: Self/Cross/Top shadowing correction (CSAT3 formulation) +#' +#' @author +#' John Frank \cr +#' David Durden \email{ddurden@battelleecology.org} +#' +#' @description +#' Applies the Self/Cross/Top (SCT) shadowing correction (Frank, 2025) using CSAT3 +#' transducer geometry. Input wind is in sonic coordinates (N × 3). Output includes +#' corrected wind (sonic and transducer frames), per-transducer SCT correction +#' coefficients, wind–transducer angles, and wind–top-structure angles. +#' +#' @param inpVeloSoni Matrix or data.frame (N × 3) of wind components in sonic +#' coordinates \[m s^-1\]; columns correspond to (veloXaxs [u], veloYaxs [v], veloZaxs [w]). +#' If a data.frame is provided, all three columns must be numeric and will be coerced +#' to a matrix. If names match \code{c("veloXaxs","veloYaxs","veloZaxs")}, they will be +#' ordered accordingly. +#' @param MtrxSoniTran 3 × 3 matrix mapping sonic → transducer coordinates +#' (columns are unit transducer vectors in sonic coords). CSAT3 by default. +#' @param MtrxTranSoni 3 × 3 matrix mapping transducer → sonic coordinates +#' (inverse of \code{MtrxSoniTran}). Default matches the inverse of the default above. +#' @param VectTop Numeric length-3 unit vector for “top/bottom” (support structure) +#' direction (e.g., \code{c(0, 0, 1)}). +#' @param ParaD0 Numeric vector (length 3): model parameter d0 for self (1), cross (2), top (3). +#' Default: \code{c(0.00159718301425092, 0.0113301152734556, 0.0633033802230689)}. +#' @param ParaD1 Numeric vector (length 3): model parameter d1 for self (1), cross (2), top (3). +#' Default: \code{c(29.2463049972615, 24.8602835257800, 18.5970206839229)}. +#' @param ParaD2 Numeric vector (length 3): model parameter d2 for self (1), cross (2), top (3). +#' Default: \code{c(0.836978194375597, 0.874138513424717, 0.809209610043660)}. +#' @param ParaD3 Positive numeric scalar: global scaling to minimize change in horizontal +#' wind relative to Wyngaard correction. Default: \code{1.06606236917833}. +#' @param ParaP Positive numeric scalar: Minkowski exponent used to combine components +#' (default 8, following your script). +#' +#' @return A list with: +#' \itemize{ +#' \item \code{veloWindSoniCorr}: N × 3 data.frame (sonic coords) after SCT correction; +#' columns \code{veloXaxs, veloYaxs, veloZaxs}. +#' \item \code{veloWindTranCorr}: N × 3 matrix (transducer coords) after SCT correction. +#' \item \code{coefSct}: N × 3 matrix of SCT correction coefficients per transducer. +#' \item \code{angTranDeg}: N × 3 matrix of angles (deg) between wind and transducer axes. +#' \item \code{angTopDeg}: N × 1 numeric vector of angles (deg) between wind and \code{VectTop}. +#' } +#' +#' @details +#' Angles \eqn{\theta} are computed between the wind vector and (i) each transducer axis, +#' and (ii) the top-structure direction. For a given parameter triple (d0, d1, d2), +#' the model error is +#' \deqn{ +#' \mathrm{err}(\theta) = \frac{(1 - d_2)\,\tanh\{d_0 (\theta - d_1)\}} +#' {\tanh\{d_0 (90 - d_1)\}} + d_2 \, . +#' } +#' The correction coefficient is \eqn{c(\theta) = 1 / \mathrm{err}(\theta)}. +#' Self, cross, and top terms are combined via a Minkowski-like synthesis with exponent +#' \code{ParaP}, then shifted (+1) and scaled by \code{ParaD3}. The correction is applied +#' element-wise in the transducer frame and mapped back to sonic coordinates. +#' +#' @examples +#' set.seed(1) +#' V <- matrix(rnorm(300), ncol = 3) +#' out <- def.corr.sct.csat( +#' inpVeloSoni = V, +#' VectTop = c(0, 0, 1) +#' # defaults for ParaD0..ParaD3 and ParaP are used +#' ) +#' str(out) +#' +#' @keywords eddy-covariance sonic-anemometer correction csat3 sct +#' @export +def.corr.sct.csat <- function( + inpVeloSoni, + MtrxSoniTran = matrix( + c( 1/4, 0.433012701892219, 0.866025403784439, + -1/2, 0.000000000000000, 0.866025403784439, + 1/4, -0.433012701892219, 0.866025403784439), + nrow = 3, ncol = 3 + ), + MtrxTranSoni = matrix( + c( 2/3, -4/3, 2/3, + 1.154700538379252, 0.000000000000000, -1.154700538379252, + 0.384900179459751, 0.384900179459751, 0.384900179459751), + nrow = 3, ncol = 3 + ), + VectTop = c(0,0,1), + ParaD0 = c(0.00159718301425092, 0.0113301152734556, 0.0633033802230689), + ParaD1 = c(29.2463049972615, 24.8602835257800, 18.5970206839229), + ParaD2 = c(0.836978194375597, 0.874138513424717, 0.809209610043660), + ParaD3 = 1.06606236917833, + ParaP = 8 +) { + # --- coerce/validate input -------------------------------------------------- + if (is.data.frame(inpVeloSoni)) { + if (ncol(inpVeloSoni) != 3L) + stop("inpVeloSoni data.frame must have exactly 3 numeric columns.") + if (!all(vapply(inpVeloSoni, is.numeric, logical(1)))) + stop("All columns of inpVeloSoni data.frame must be numeric.") + nameVarExpc <- c("veloXaxs", "veloYaxs", "veloZaxs") + if (all(nameVarExpc %in% names(inpVeloSoni))) { + inpVeloSoni <- inpVeloSoni[, nameVarExpc] + } + inpVeloSoni <- as.matrix(inpVeloSoni) + } + stopifnot( + is.matrix(inpVeloSoni), ncol(inpVeloSoni) == 3, + is.matrix(MtrxSoniTran), all(dim(MtrxSoniTran) == c(3, 3)), + is.matrix(MtrxTranSoni), all(dim(MtrxTranSoni) == c(3, 3)), + is.numeric(VectTop), length(VectTop) == 3, + is.numeric(ParaD0), length(ParaD0) == 3, + is.numeric(ParaD1), length(ParaD1) == 3, + is.numeric(ParaD2), length(ParaD2) == 3, + is.numeric(ParaD3), ParaD3 > 0, + is.numeric(ParaP), ParaP > 0 + ) + numRow <- nrow(inpVeloSoni) + if (numRow == 0) stop("inpVeloSoni has zero rows") + + # --- map to transducer frame ----------------------------------------------- + veloWindTran <- inpVeloSoni %*% MtrxSoniTran # N × 3 + + # wind magnitude for angle calc (avoid div-by-zero) + magVeloSoni <- sqrt(rowSums(inpVeloSoni^2)) # N + denm <- pmax(magVeloSoni, .Machine$double.eps) + + # --- angles (degrees) ------------------------------------------------------- + # wind–transducer axes + angCosTran <- abs(veloWindTran) / denm # N × 3 (recycling over columns) + angCosTran[angCosTran > 1] <- 1 + angTranDeg <- acos(angCosTran) * 180 / pi # N × 3 + + # wind–top direction + angCosTop <- abs(drop(inpVeloSoni %*% VectTop)) / denm # N + angCosTop[angCosTop > 1] <- 1 + angTopDeg <- acos(angCosTop) * 180 / pi # N + + # --- helper: error model -> coefficient ------------------------------------ + def.err.sct.tanh <- function(angDeg, ParaD0j, ParaD1j, ParaD2j) { + # angDeg can be N×3 (matrix) or N (vector) + (1 - ParaD2j) * tanh(ParaD0j * (angDeg - ParaD1j)) / tanh(d0j * (90 - ParaD1j)) + ParaD2j + } + + # --- evaluate error & coefficients ----------------------------------------- + errSelf <- def.err.sct.tanh(angTranDeg, ParaD0[1], ParaD1[1], ParaD2[1]) # N × 3 + errCros <- def.err.sct.tanh(angTranDeg, ParaD0[2], ParaD1[2], ParaD2[2]) # N × 3 + errTop <- def.err.sct.tanh(angTopDeg, ParaD0[3], ParaD1[3], ParaD2[3]) # N + + coefSelf <- 1 / errSelf + coefCros <- 1 / errCros + coefTop <- 1 / errTop + + # --- Minkowski-like synthesis per transducer arm --------------------------- + coefSct <- matrix(0, numRow, 3) + # Column 1 combines: self@1, cross@2 & @3, and top + coefSct[, 1] <- (((coefSelf[, 1] - 1)^ParaP + + (coefCros[, 2] - 1)^ParaP + + (coefCros[, 3] - 1)^ParaP + + (coefTop - 1)^ParaP)^(1/ParaP)) + 1 + # Column 2 combines: self@2, cross@1 & @3, and top + coefSct[, 2] <- (((coefCros[, 1] - 1)^ParaP + + (coefSelf[, 2] - 1)^ParaP + + (coefCros[, 3] - 1)^ParaP + + (coefTop - 1)^ParaP)^(1/ParaP)) + 1 + # Column 3 combines: self@3, cross@1 & @2, and top + coefSct[, 3] <- (((coefCros[, 1] - 1)^ParaP + + (coefCros[, 2] - 1)^ParaP + + (coefSelf[, 3] - 1)^ParaP + + (coefTop - 1)^ParaP)^(1/ParaP)) + 1 + + # global scaling + coefSct <- coefSct / ParaD3 + + # --- apply correction in transducer frame ---------------------------------- + veloWindTranCorr <- veloWindTran * coefSct # element-wise + veloWindSoniMat <- veloWindTranCorr %*% MtrxTranSoni # back to sonic frame + + # --- convert sonic output to data.frame ------------------------------------ + veloWindSoniCorr <- as.data.frame(veloWindSoniMat) + names(veloWindSoniCorr) <- c("veloXaxs", "veloYaxs", "veloZaxs") + # rownames(veloWindSoniCorr) <- NULL # optional + + # --- return ---------------------------------------------------------------- + list( + veloWindSoniCorr = veloWindSoniCorr, # data.frame (N × 3) + veloWindTranCorr = veloWindTranCorr, # matrix (N × 3) + coefSct = coefSct, # matrix (N × 3) + angTranDeg = angTranDeg, # matrix (N × 3) + angTopDeg = angTopDeg # vector (N) + ) +} diff --git a/pack/eddy4R.turb/R/def.corr.wyzh.csat.R b/pack/eddy4R.turb/R/def.corr.wyzh.csat.R new file mode 100644 index 00000000..7d2b85ea --- /dev/null +++ b/pack/eddy4R.turb/R/def.corr.wyzh.csat.R @@ -0,0 +1,137 @@ + + +#------------------------------------------------------------------------------- +# def.corr.wyzh.csat +# +# Apply Wyngaard & Zhang (1985) sinusoidal shadowing correction in the +# Horst et al. (2015) formulation for CSAT3 geometry. Converts sonic-frame +# wind to transducer frame, applies element-wise correction, and maps back. +# +#------------------------------------------------------------------------------- + +#' @title Definition function: Wyngaard–Zhang sinusoidal shadowing correction (CSAT3 formulation) +#' +#' @author +#' John Frank \cr +#' David Durden \email{ddurden@battelleecology.org} +#' +#' +#' @description +#' Applies the Wyngaard & Zhang (1985) sinusoidal shadowing correction using the +#' Horst et al. (2015) formulation for CSAT3 transducer geometry. Input wind is +#' in sonic coordinates (N × 3). Output includes corrected wind (sonic and +#' transducer frames), per-transducer correction coefficients, and wind–transducer +#' angles. +#' +#' @param inpVeloSoni Matrix or data.frame (N × 3) of wind components in sonic +#' coordinates \[m s^-1\]; columns correspond to (veloXaxs [u], veloYaxs [v], veloZaxs [w]). If a data.frame is +#' provided, all three columns must be numeric and will be coerced to a matrix. +#' @param MtrxSoniTran 3 × 3 matrix mapping sonic → transducer coordinates +#' (columns are unit transducer vectors in sonic coords). Default is CSAT3. +#' @param MtrxTranSoni 3 × 3 matrix mapping transducer → sonic coordinates +#' (inverse of \code{MtrxSoniTran}). Default matches the inverse of the default above. +#' +#' @return A list with: +#' \itemize{ +#' \item \code{veloWindSoniCorr}: N × 3 data.frame (sonic coords) after correction. +#' \item \code{veloWindTranCorr}: N × 3 matrix (transducer coords) after correction. +#' \item \code{coefWyzh}: N × 3 element-wise correction coefficients per transducer. +#' \item \code{angTranDeg}: N × 3 angles (deg) between wind and transducer axes. +#' } +#' +#' @details +#' For each sample and transducer axis, the angle \eqn{\theta} between the wind +#' vector and the axis is computed, and the Wyngaard sinusoidal coefficient is +#' \deqn{ c(\theta) = 1 / (0.84 + 0.16 \sin \theta) .} +#' The correction is applied element-wise in the transducer frame, and the +#' corrected vectors are mapped back to sonic coordinates. +#' +#' @examples +#' set.seed(1) +#' V <- matrix(rnorm(300), ncol = 3) +#' out <- def.corr.wyzh.csat(V) +#' str(out) +#' +#' # With a data.frame: +#' DF <- as.data.frame(V) +#' out2 <- def.corr.wyzh.csat(DF) +#' all.equal(out, out2) # should be TRUE +#' +#' @keywords eddy-covariance sonic-anemometer correction csat3 wyngaard zhang +#' @export +#' +#' +#' +def.corr.wyzh.csat <- function( + inpVeloSoni, + MtrxSoniTran = matrix( + c( 1/4, 0.433012701892219, 0.866025403784439, + -1/2, 0.000000000000000, 0.866025403784439, + 1/4, -0.433012701892219, 0.866025403784439), + nrow = 3, ncol = 3 + ), + MtrxTranSoni = matrix( + c( 2/3, -4/3, 2/3, + 1.154700538379252, 0.000000000000000, -1.154700538379252, + 0.384900179459751, 0.384900179459751, 0.384900179459751), + nrow = 3, ncol = 3 + ) +) { + # --- coerce/validate input -------------------------------------------------- + if (is.data.frame(inpVeloSoni)) { + if (ncol(inpVeloSoni) != 3L) + stop("inpVeloSoni data.frame must have exactly 3 columns.") + if (!all(vapply(inpVeloSoni, is.numeric, logical(1)))) + stop("All columns of inpVeloSoni data.frame must be numeric.") + + + # Reorder if named as expected + nameVarExpc <- c("veloXaxs", "veloYaxs", "veloZaxs") + if (all(nameVarExpc %in% names(inpVeloSoni))) { + inpVeloSoni <- inpVeloSoni[,nameVarExpc] + } + + inpVeloSoni <- as.matrix(inpVeloSoni) + } + stopifnot( + is.matrix(inpVeloSoni), ncol(inpVeloSoni) == 3, + is.matrix(MtrxSoniTran), all(dim(MtrxSoniTran) == c(3, 3)), + is.matrix(MtrxTranSoni), all(dim(MtrxTranSoni) == c(3, 3)) + ) + numRow <- nrow(inpVeloSoni) + if (numRow == 0) stop("inpVeloSoni has zero rows") + + # --- map to transducer frame ----------------------------------------------- + # veloWindTran = V_sonic · (unit transducer vectors in sonic coords) + veloWindTran <- inpVeloSoni %*% MtrxSoniTran # N × 3 + + # wind magnitude for angle calc (avoid div-by-zero) + magVeloSoni <- sqrt(rowSums(inpVeloSoni^2)) # N + denm <- pmax(magVeloSoni, .Machine$double.eps) # N + + # --- wind–transducer angles (degrees) -------------------------------------- + angCosTran <- abs(veloWindTran) / denm # vectorized over columns + angCosTran[angCosTran > 1] <- 1 # numeric safety + angTranDeg <- acos(angCosTran) * 180 / pi # N × 3 + + # --- Wyngaard–Zhang sinusoidal coefficient --------------------------------- + coefWyzh <- 1 / (0.84 + 0.16 * sin(angTranDeg * pi / 180)) # N × 3 + + # --- apply correction in transducer frame ---------------------------------- + veloWindTranCorr <- veloWindTran * coefWyzh # element-wise correction + veloWindSoniCorr <- veloWindTranCorr %*% MtrxTranSoni # back to sonic frame + + # --- convert output data.frame ---------------------------------- + veloWindSoniCorr <- as.data.frame(veloWindSoniCorr) + names(veloWindSoniCorr) <- c("veloXaxs", "veloYaxs", "veloZaxs") + rownames(veloWindSoniCorr) <- NULL + + + # --- return ---------------------------------------------------------------- + list( + veloWindSoniCorr = veloWindSoniCorr, + veloWindTranCorr = veloWindTranCorr, + coefWyzh = coefWyzh, + angTranDeg = angTranDeg + ) +}