From b7142611d66ba5e7e5249f105fa04875876a186c Mon Sep 17 00:00:00 2001 From: Beaudette Date: Wed, 16 Oct 2024 11:14:16 -0700 Subject: [PATCH] reasonable accounting for #316 , not a generic fix for dynamic res --- R/addVolumeFraction.R | 256 +++++++++++++++++++++------------------ man/addVolumeFraction.Rd | 8 +- 2 files changed, 145 insertions(+), 119 deletions(-) diff --git a/R/addVolumeFraction.R b/R/addVolumeFraction.R index f1048ebe6..39e5e29f5 100644 --- a/R/addVolumeFraction.R +++ b/R/addVolumeFraction.R @@ -1,82 +1,82 @@ ## TODO: this could use an overhaul: ## * symbol size / resolution should be determined in context -## * area approximation should work with a larger grid, after scaling to avoid frational thickness constraints +## * area approximation should work with a larger grid, after scaling to avoid fractional thickness constraints ## * standardize symbology with rect() vs. points() and cex ## -> init checkerboard, select squares, fill squares ## -> number / size of squares is the only parameter ## doesn't work with fractional depths: https://github.com/ncss-tech/aqp/issues/8 + # convert volume pct [0, 100] into DF of points along a res x res grid +# v should not be NA .volume2df <- function(v, depth, res) { - # test for no data - if(is.na(v)) - return(data.frame()) - + # test for >= 100% if(v >= 100) { warning(sprintf("%s is >= 100, likely a data error", v), call. = FALSE) # truncate at 100% v <- 100 } - - # convert volume pct into fraction - v <- v / 100 - - # init matrix with NA - cells <- (depth * res) - m <- matrix(nrow=depth, ncol=res) - m[] <- NA - - # determine number of cells required to symbolize volume fraction - v.n <- round(v * cells) - v.cells <- sample(1:cells, size=v.n) - # mark cells with 1 - m[v.cells] <- 1 - - # convert matrix into data.frame - d <- expand.grid(x=(1:res), y=(1:depth)) - d$val <- m[1:cells] - # keep only those cells with val = 1 - d <- d[which(d$val == 1), ] - - # scrub extra columns and return - d$val <- NULL - return(d) - } + + # convert volume pct into fraction + v <- v / 100 + + # init matrix with NA + cells <- (depth * res) + m <- matrix(nrow = depth, ncol = res) + # set to 1 for debugging output (full grid) + m[] <- NA + + # determine number of cells required to symbolize volume fraction + # v must be >0 -> select a single cell when v is <1 + v.n <- pmax(round(v * cells), 1) + + # random sample of cells, mark with 1 in matrix + v.cells <- sample(1:cells, size = v.n) + m[v.cells] <- 1 + + # convert matrix into data.frame + d <- expand.grid(x = (1:res), y = (1:depth)) + d$val <- m[1:cells] + # keep only those cells with val = 1 + d <- d[which(d$val == 1), ] + + # scrub extra columns and return + d$val <- NULL + return(d) +} -#' @title Symbolize Volume Fraction on a Soil Profile Collection Plot +#' @title Symbolize Volume Fraction within a Soil Profile Collection Plot #' #' @description Symbolize volume fraction on an existing soil profile collection plot. #' -#' @param x a \code{SoilProfileCollection} object -#' @param colname character vector of length 1, naming the column containing volume fraction data (horizon-level attribute) +#' @param x a `SoilProfileCollection` object +#' @param colname character vector of length 1, naming the column containing volume fraction data (horizon-level attribute). Values should be within 0-100 percent. #' @param res integer, resolution of the grid used to symbolize volume fraction #' @param cex.min minimum symbol size #' @param cex.max maximum symbol size #' @param pch integer, plotting character code #' @param col symbol color, either a single color or as many colors as there are horizons in `x` #' -#' @details This function can only be called after plotting a \code{SoilProfileCollection} object. Details associated with a call to \code{plotSPC} are automatically accounted for within this function: e.g. \code{plot.order}, \code{width}, etc.. +#' @details This function can only be called after plotting a `SoilProfileCollection` object. Details associated with a call to `plotSPC()` are automatically accounted for within this function: e.g. `plot.order`, `width`, etc.. #' #' @author D.E. Beaudette #' -#' @seealso \code{\link{plotSPC}} +#' @seealso [plotSPC()] #' #' @export #' ## TODO: symbol size must be controlled by `res` -## TODO: fails with fractional depths https://github.com/ncss-tech/aqp/issues/8 - -addVolumeFraction <- function(x, colname, res=10, cex.min=0.1, cex.max=0.5, pch=1, col='black') { - +addVolumeFraction <- function(x, colname, res = 10, cex.min = 0.1, cex.max = 0.5, pch = 1, col = 'black') { + # color should be either: # single color name / code # vector of colors with length == nrow(x) - + # simplest case, single color vector if(length(col) == 1) { col <- rep(col, times=nrow(x)) @@ -85,90 +85,115 @@ addVolumeFraction <- function(x, colname, res=10, cex.min=0.1, cex.max=0.5, pch= if(length(col) != nrow(x)) stop('length of `col` should be either 1 or nrow(x)', call. = FALSE) } - + # ensure that `colname` is a horizon-level attribute if(! colname %in% horizonNames(x)) { stop(sprintf("%s is not a horizon-level attribute", colname), call. = FALSE) } - + # test for values < 0.5, could be a fraction [0,1] vs. percent [0,100] if(all( na.omit(horizons(x)[[colname]]) < 0.5) ) { message(sprintf("all %s values are < 0.5, likely a fraction vs. percent", colname)) } - - # get plotting details from aqp environment - lsp <- get('last_spc_plot', envir=aqp.env) - w <- lsp$width - plot.order <- lsp$plot.order - # y.offset is a vector of length(x) - depth.offset <- lsp$y.offset - sf <- lsp$scaling.factor - x0 <- lsp$x0 - - # horizontal shrinkage factor - ## TODO: why is this hard-coded ? - w.offset <- w / 7 - - # get top/bottom colnames - hd <- horizonDepths(x) - - # iterate over profiles - for(p.i in 1:length(x)) { - - # get the current pofile, in plotting order - h <- horizons(x[plot.order[p.i], ]) - - ## determine left/right extent of symbols - # 2019-07-15: using relative position as indexed by current profile - x.center <- x0[p.i] - x.left <- x.center - (w - w.offset) - x.right <- x.center + (w - w.offset) - - # iterate over horizons - for(h.i in 1:nrow(h)) { - this.hz <- h[h.i, ] - hz.thick <- this.hz[[hd[2]]] - this.hz[[hd[1]]] - - ## hack until #8 is resolved: round the thickness down - # https://github.com/ncss-tech/aqp/issues/8 - if( hz.thick %% 1 != 0) { - msg <- paste0(paste(h[h.i, hd], collapse = ' - '), depth_units(x)) - message(sprintf('truncating fractional horizon thickness to integer: %s', msg)) - hz.thick <- floor(hz.thick) - } - - - # convert this horizon's data - v <- .volume2df(v=this.hz[[colname]], depth=hz.thick, res=res) - - ## TODO: still throws errors - # just those with data - if(nrow(v) > 0) { - - # get the current color from vector of colors - # typically the same color repeated, but could be as many colors as hz - v$color <- col[h.i] - + + # get plotting details from aqp environment + lsp <- get('last_spc_plot', envir = aqp.env) + w <- lsp$width + plot.order <- lsp$plot.order + + # y.offset is a vector of length(x) + depth.offset <- lsp$y.offset + sf <- lsp$scaling.factor + x0 <- lsp$x0 + + # horizontal shrinkage factor + # hard-coded based on aesthetic tests + w.offset <- w / 7 + + # get top/bottom colnames + hd <- horizonDepths(x) + + # iterate over profiles + for(p.i in 1:length(x)) { + + # get the current pofile, in plotting order + h <- horizons(x[plot.order[p.i], ]) + + ## determine left/right extent of symbols + # 2019-07-15: using relative position as indexed by current profile + x.center <- x0[p.i] + x.left <- x.center - (w - w.offset) + x.right <- x.center + (w - w.offset) + + # iterate over horizons + for(h.i in 1:nrow(h)) { + this.hz <- h[h.i, ] + hz.thick <- this.hz[[hd[2]]] - this.hz[[hd[1]]] + + ## hack until #8 is resolved: round the thickness down + # https://github.com/ncss-tech/aqp/issues/8 + if( hz.thick %% 1 != 0) { + # truncate at closest integer + # minimum value of 1 to ensure something is plotted within very thin horizons + hz.thick <- pmax(floor(hz.thick), 1) + + msg <- sprintf( + "%s%s -> %s%s thick", + paste(h[h.i, hd], collapse = ' - '), + depth_units(x), + hz.thick, + depth_units(x) + ) + message(sprintf('truncating fractional horizon thickness to integer: %s', msg)) + + } + + # proceed only if the current profile/horizon: + # * value not NA + # * value > 0 + # * hz.thick not NA + # * hz.thick > 0 + + if((hz.thick > 0) & (!is.na(hz.thick)) & (this.hz[[colname]] > 0) & (! is.na(this.hz[[colname]]))) { + + # convert this horizon's data + v <- .volume2df(v = this.hz[[colname]], depth = hz.thick, res = res) + + ## TODO: this can be better solved with dynamic 'res' + # if (thickness * volume * res) is too small, there may be no samples returned + # skip + if(nrow(v) < 1) { + # print('skipping...') + next + } + + # get the current color from vector of colors + # typically the same color repeated, but could be as many colors as hz + # may be NA if using a horizon column + v$color <- col[h.i] + # jitter and rescale x-coordinates - v$x <- .rescaleRange(jitter(v$x), x0 = x.left, x1 = x.right) - - # determine horizon depths in current setting - # depth_prime = (depth * scaling factor) + y.offset[i] - y.top <- (this.hz[[hd[1]]] * sf) + depth.offset[p.i] - y.bottom <- (this.hz[[hd[2]]] * sf) + depth.offset[p.i] - # rescale y-coordinates - v$y <- .rescaleRange(jitter(v$y), x0 = y.top, x1 = y.bottom) - - # generate random symbol size - p.cex <- runif(nrow(v), min=cex.min, max=cex.max) - - # add jittered points - # note that color comes from `v` - points(v$x, v$y, cex=p.cex, col=v$color, pch=pch) - } - } # end looping over profiles - - } + v$x <- .rescaleRange(jitter(v$x), x0 = x.left, x1 = x.right) + + # determine horizon depths in current setting + # depth_prime = (depth * scaling factor) + y.offset[i] + y.top <- (this.hz[[hd[1]]] * sf) + depth.offset[p.i] + y.bottom <- (this.hz[[hd[2]]] * sf) + depth.offset[p.i] + + # rescale y-coordinates + v$y <- .rescaleRange(jitter(v$y), x0 = y.top, x1 = y.bottom) + + # generate random symbol size + p.cex <- runif(nrow(v), min = cex.min, max = cex.max) + + # add jittered points + # note that color comes from `v` + points(v$x, v$y, cex = p.cex, col = v$color, pch = pch) + } + + } # end iteration over horizons + } # end iteration over profiles + } @@ -178,3 +203,4 @@ addVolumeFraction <- function(x, colname, res=10, cex.min=0.1, cex.max=0.5, pch= + diff --git a/man/addVolumeFraction.Rd b/man/addVolumeFraction.Rd index a2041ce8f..5169230bf 100644 --- a/man/addVolumeFraction.Rd +++ b/man/addVolumeFraction.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/addVolumeFraction.R \name{addVolumeFraction} \alias{addVolumeFraction} -\title{Symbolize Volume Fraction on a Soil Profile Collection Plot} +\title{Symbolize Volume Fraction within a Soil Profile Collection Plot} \usage{ addVolumeFraction( x, @@ -17,7 +17,7 @@ addVolumeFraction( \arguments{ \item{x}{a \code{SoilProfileCollection} object} -\item{colname}{character vector of length 1, naming the column containing volume fraction data (horizon-level attribute)} +\item{colname}{character vector of length 1, naming the column containing volume fraction data (horizon-level attribute). Values should be within 0-100 percent.} \item{res}{integer, resolution of the grid used to symbolize volume fraction} @@ -33,10 +33,10 @@ addVolumeFraction( Symbolize volume fraction on an existing soil profile collection plot. } \details{ -This function can only be called after plotting a \code{SoilProfileCollection} object. Details associated with a call to \code{plotSPC} are automatically accounted for within this function: e.g. \code{plot.order}, \code{width}, etc.. +This function can only be called after plotting a \code{SoilProfileCollection} object. Details associated with a call to \code{plotSPC()} are automatically accounted for within this function: e.g. \code{plot.order}, \code{width}, etc.. } \seealso{ -\code{\link{plotSPC}} +\code{\link[=plotSPC]{plotSPC()}} } \author{ D.E. Beaudette