R/threeDspectro.R

#' 3D spectrogram plots from \code{"Wave"} objects
#'
#' @description
#' Create 3D spectrograms from a single object of class \code{"Wave"}.
#'
#' This function works similarly as \code{\link{spectro}} (\code{\link{seewave}} package), with spectrogram data internally computed by \code{\link{spectro}}. However, the 3D plot is generated by \code{\link{persp3D}} (\code{\link{plot3D}} package).
#'
#' @param wave an \code{R} object of the class \code{"Wave"}
#' @param tlim modifications of the time limits (X-axis). Vector with two values in seconds. By default: \code{tlim = NULL}
#' @param flim modifications of the frequency limits (Y-axis). Vector with two values in kHz. By default: \code{flim = NULL}
#' @param samp.grid a logical. If \code{TRUE}, a sampling grid with dimensions \code{x.length} and \code{y.length} is applied over the three-dimensional surface. An amplitude value is collected at each node in the grid and a new matrix of amplitude coordinates is used in the plot. Recommended for faster plots and as protocol for error verification in point acquisition. See also \code{plot.type}. By default: \code{samp.grid = FALSE}
#' @param plot.type allows the choice between `"lines"`, `"surface"` or `"points"`. Beware that `"points"` only applies when \code{samp.grid = TRUE}. \code{plot.type = "surface"} will produce a simplified 3D spectrogram surface using the amplitude values colected by the sampling grid (same output employed by MacLeod et al., 2013). \code{plot.type = "points"} will produce 3D graphs with semilandmarks as points. \code{plot.type = "lines"} will generate profile lines along X or Y-axis (see argument \code{along}). Number and space between lines are controled by arguments \code{skip.lines} and \code{space.lines}, respectively. By default: \code{plot.type = "surface"}
#' @param x.length only applies when \code{samp.grid = TRUE}. Length of sequence (i.e. number of cells per side on sound window) to be used as sampling grid coordinates on the time (X-axis).  By default: \code{x.length = 100}
#' @param y.length only applies when \code{samp.grid = TRUE}. Length of sequence (i.e. number of cells per side on sound window) to be used as sampling grid coordinates on the frequency (Y-axis). By default: \code{y.length = 70}
#' @param log.scale only applies when \code{samp.grid = TRUE}. A logical. If \code{TRUE}, \code{threeDspectro} will use a logarithmic scale on the time (X-axis), which is recommeded when the analyzed sounds present great variation on this axis (e.g. emphasize short duration sounds). If \code{FALSE}, a linear scale is used instead (same as MacLeod et al., 2013). By default: \code{log.scale = FALSE}
#' @param along only applies when \code{plot.type = "lines"}. Lines along X or Y-axis (\code{"x"} or \code{"y"}). By default: \code{along = "x"}
#' @param skip.lines only applies when \code{plot.type = "lines"}. How many lines to skip between each plotted line. Larger numbers imply on less lines plotted. By default: \code{skip.lines = 5}
#' @param space.lines only applies when \code{plot.type = "lines"}. Controls the space between lines, which makes lines wider or narrower. Has to be a number between 0.1 and 0.9. By default: \code{space.lines = 0.6}
#' @param cex only applies when \code{samp.grid = TRUE} and \code{plot.type = "points"}. Similarly as in \code{\link{par}}, intended size for points. By default: \code{cex = 0.5}
#' @param cex.axis Similarly as in \code{\link{par}}, the magnification to be used for axis annotation. By default: \code{cex.axis = 0.5}
#' @param cex.lab Similarly as in \code{\link{par}}, the magnification to be used for x and y labels. By default: \code{cex.lab = 0.9}
#' @param cex.main Similarly as in \code{\link{par}}, the magnification to be used for main titles. By default: \code{cex.main = 1}
#' @param lwd only applies when \code{samp.grid = TRUE} and \code{plot.type = "surface"}. Similarly as in \code{\link{par}}, intended line width for sampling grid. By default: \code{lwd = 0.1}
#' @param plot.exp a logical. If \code{TRUE}, 3D spectrogram plot from \code{wave} object is exported and stored on the folder indicated by \code{store.at}. By default: \code{plot.exp = FALSE}
#' @param store.at only applies when \code{plot.exp = TRUE}. Filepath to the folder where 3D plot will be stored. Should be presented between quotation marks. By default: \code{store.at = NULL} (i.e. user must specify the filepath where spectrogram plots will be stored)
#' @param plot.as only applies when \code{plot.exp = TRUE}. \code{plot.as = "jpeg"} will generate compressed images for quick inspection; \code{plot.as = "tiff"} or \code{"tif"} will generate uncompressed high resolution images that can be edited and used for publication. By default: \code{plot.as = "jpeg"}
#' @param color Color palette to be used for the amplitude (Z-axis). Same default as \code{\link{spectro}}: \code{color=spectro.colors(80)}. See also \code{Details} section.
#' @param scalelab Similarly as \code{\link{plot3D}}, the label to be written on top of the color key. Should be a character string wrapped by \code{expression()}. By default: \code{scalelab = expression("Amplitude (dB)")}. See also \code{\link{colkey}}
#' @param colkey Similarly as \code{\link{plot3D}}, a list with parameters for the color key (legend). By default: \code{colkey = list(plot = TRUE, cex.clab = 0.8, cex.axis = 1, side = 4, length = 1, width = 1, labels = TRUE, tick = TRUE, lty = 1, lwd = 1, lwd.ticks = 1)}. See also \code{\link{colkey}}
#' @param f sampling frequency of \code{Wave} object (in Hz). By default: \code{f = 44100}
#' @param wl length of the window for spectrogram calculation. By default: \code{wl = 512}
#' @param ovlp overlap between two successive windows (in %) for increased spectrogram resolution. By default: \code{ovlp = 70}
#' @param dBlevel absolute amplitude value to be used as relative background on 3D plot. Same as \code{dBlevel} from \code{\link{eigensound}} and \code{\link{align.wave}}. By default: \code{dBlevel = 30}
#' @param resfac resolution factor, in which an value > 1 will increase the resolution. Can be one value or a vector of two numbers, for the x and y values, respectively. \strong{Note:} Same as in \code{\link{persp3D}} (\code{\link{plot3D}} package). By default: \code{resfac = 1}
#' @param rotate.Xaxis rotation of the X-axis. Same as \code{theta} from \code{\link{persp3D}} (\code{\link{plot3D}} package). By default: \code{rotate.Xaxis = 60}
#' @param rotate.Yaxis rotation of the Y-axis. Same as \code{phi} from \code{\link{persp3D}} (\code{\link{plot3D}} package). By default: \code{rotate.Yaxis = 40}
#' @param main main title of output plot. Should be presented between quotation marks. By default: \code{main = "3D spectrogram"}
#'
#' @details
#' Similarly as \code{\link{spectro}} (\code{\link{seewave}} package), any colour palette can be used to describe the amplitude (Z-axis). Some suggestions: \code{seewave::temp.colors, seewave::spectro.colors, seewave::reverse.heat.colors, seewave::reverse.cm.colors, seewave::reverse.topo.colors, grDevices::cm.colors, grDevices::grey.colors, grDevices::heat.colors, grDevices::topo.colors}.
#'
#' @author
#' Pedro Rocha
#'
#' @references
#' MacLeod, N., Krieger, J. & Jones, K. E. (2013). Geometric morphometric approaches to acoustic signal analysis in mammalian biology. \emph{Hystrix, the Italian Journal of Mammalogy, 24}(1), 110-125.
#'
#' Rocha, P. & Romano, P. (2021) The shape of sound: A new \code{R} package that crosses the bridge between Bioacoustics and Geometric Morphometrics. \emph{Methods in Ecology and Evolution, 12}(6), 1115-1121.
#'
#' @seealso
#' \code{\link{spectro}}, \code{\link{seewave}}, \code{\link{eigensound}}, \code{\link{align.wave}}, \code{\link{persp3D}}, \code{\link{plot3D}}, \code{\link{align.wave}}
#'
#' Useful links:
#' \itemize{
#'   \item{\url{https://github.com/p-rocha/SoundShape}}
#'   \item{Report bugs at \url{https://github.com/p-rocha/SoundShape/issues}}}
#'
#'
#' @examples
#'
#' \donttest{
#' # As simple as this
#' threeDspectro(centralis)
#' threeDspectro(cuvieri)
#' threeDspectro(kroyeri)
#'
#' # Controling some arguments
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4))
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4), dBlevel=50)
#'
#' # As points
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4),
#'               samp.grid=TRUE, plot.type="points")
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4),
#'               samp.grid=TRUE, plot.type="points", x.length = 20, y.length = 50)
#'
#' # As lines
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4), plot.type = "lines")
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4),
#'               plot.type = "lines", along="y")
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4),
#'               plot.type = "lines", along="y", skip.lines=18, space.lines=0.8)
#'
#' # Try different colors
#' threeDspectro(cuvieri, color=seewave::reverse.terrain.colors(80),
#'               samp.grid=FALSE, tlim=c(0, 0.5), flim=c(0, 4))
#' threeDspectro(cuvieri, color=seewave::reverse.cm.colors(80),
#'               samp.grid=FALSE, tlim=c(0, 0.5), flim=c(0, 4))
#' threeDspectro(cuvieri, color=grDevices::heat.colors(80),
#'               samp.grid=FALSE, tlim=c(0, 0.5), flim=c(0, 4))
#'
#' # Rotation
#' threeDspectro(cuvieri, tlim=c(0, 0.5), flim=c(0, 4), rotate.Xaxis=40, rotate.Yaxis=50)
#'
#' # Export your graph
#' threeDspectro(cuvieri, plot.exp=TRUE, store.at=tempdir(), tlim=c(0,0.5), flim=c(0,4))
#' }
#'
#' @export
#'
threeDspectro <-  function (wave, tlim = NULL, flim = NULL, samp.grid = FALSE, plot.type = "surface", along="x", skip.lines=5, space.lines = 0.6, x.length=100, y.length=70, lwd=0.1, plot.exp=FALSE, log.scale = FALSE, cex = 0.5, cex.axis=0.5, cex.lab=0.9, cex.main=1, store.at = NULL, plot.as = "jpeg", color = seewave::spectro.colors(80), f = 44100, wl = 512, ovlp = 70, dBlevel = 30, resfac = 1, rotate.Xaxis = 60, rotate.Yaxis = 40, main = "Spectrogram 3D", scalelab=expression("Amplitude (dB)"), colkey = list(plot = TRUE, cex.clab = 0.8, cex.axis = 1, side = 4, length = 1, width = 1, labels = TRUE, tick = TRUE, lty = 1, lwd = 1, lwd.ticks = 1))
{

  if(plot.exp == TRUE && is.null(store.at)){
    stop("Use 'store.at' to specify folder path where 3D spectrogram plots will be stored")}

  if(plot.type == "points" && samp.grid == FALSE){
    stop(" 'plot.type = TRUE' require 'samp.grid = TRUE'. You can also specify 'x.length' and 'y.length', which will control the density of points/semilandmarks ploted")}


  # Acquire spectrogram data
  s <- seewave::spectro(wave, plot=F, f=f, wl=wl, ovlp=ovlp, tlim=tlim, flim=flim)

  # Set background
  for(i in 1:length(s$amp)){if(s$amp[i] == -Inf |s$amp[i] <= -dBlevel)
  {s$amp[i] <- -dBlevel}}

  if(samp.grid==TRUE){

    # Create new sequence to use as coordinates
    freq.seq <- seq(1, length(s$freq), length = y.length)
    ifelse(isTRUE(log.scale),
           time.seq <- 10^(seq(log10(1), log10(length(s$time)), length.out = x.length)),#log
           time.seq <- seq(1, length(s$time), length.out = x.length)) # linear scale

    # Subset original coordinates by new sequence
    time.sub <- s$time[time.seq]
    freq.sub <- s$freq[freq.seq]

    # Subset matrix of amplitude values using new sequences
    amp.sub <- s$amp[freq.seq, time.seq]

    # Assign time and frequency coordinates as column and row names of amplitude matrix
    colnames(amp.sub) <- time.sub
    rownames(amp.sub) <- freq.sub

    # Transform amplitude matrix into semilandmark 3D coordinates
    ind.3D <- as.matrix(stats::setNames(reshape2::melt(t(amp.sub)), c('time', 'freq', 'amp')))

  } # end samp.grid=TRUE


  # Plot semilandmarks as points or sound surface
  if(plot.exp==TRUE){
    if(plot.as == "jpeg")
    {grDevices::jpeg(width =5000,height = 3500, units = "px", res = 500,
                     filename=paste(store.at,"/", main, ".jpg", sep=""))} # compressed images
    if(plot.as=="tiff"|plot.as=="tif")
    {grDevices::tiff(width=5000, height=3500, units="px", res=500,
                     filename=paste(store.at,"/", main, ".tif", sep=""))} # uncompressed images

  } # end plot.exp == TRUE

  if(!isTRUE(samp.grid)){

    if(plot.type=="surface")
    {
      plot3D::persp3D(x = s$time, y = s$freq, z = t(s$amp), colkey=colkey, main=main,
                      theta = rotate.Xaxis, phi = rotate.Yaxis, resfac = resfac,
                      r = 3, expand = 0.5, scale = T, axes = T, ticktype = "detailed",                     nticks=4, col=color, cex.axis=cex.axis, cex.lab=cex.lab,
                      cex.main=cex.main, clab = scalelab, xlab = "Time (s)",
                      ylab = "Frequency (kHz)", zlab = "Amplitude (dB)")
    } # plot surface

    if(plot.type=="lines")
    {
      if(along=="x"){
        plot3D::ribbon3D(x = s$time,
                 y = s$freq[seq(1, nrow(s$amp), by=skip.lines)],
                 z = t(s$amp[ seq(1, nrow(s$amp), by=skip.lines), ]),
                 r = 3, expand = 0.5, scale = T, axes = T,
                 nticks=4, ticktype = "detailed",
                 col=color, cex.axis=cex.axis, cex.lab=cex.lab,
                 cex.main=cex.main, clab = scalelab, xlab = "Time (s)",
                 ylab = "Frequency (kHz)", zlab = "Amplitude (dB)",
                 space = space.lines, along = along,  colkey = colkey, main=main,
                 theta = rotate.Xaxis, phi = rotate.Yaxis)
      }


      if(along=="y"){
        plot3D::ribbon3D(x = s$time[seq(1, ncol(s$amp), by=skip.lines)],
                 y = s$freq,
                 z = t(s$amp[, seq(1, ncol(s$amp), by=skip.lines)]),
                 r = 3, expand = 0.5, scale = T, axes = T,
                 nticks=4, ticktype = "detailed",
                 col=color, cex.axis=cex.axis, cex.lab=cex.lab,
                 cex.main=cex.main, clab = scalelab, xlab = "Time (s)",
                 ylab = "Frequency (kHz)", zlab = "Amplitude (dB)",
                 space = space.lines, along = along, colkey = colkey, main=main,
                 theta = rotate.Xaxis, phi = rotate.Yaxis)
      }

    }

  } # end samp.grid = TRUE


  if(isTRUE(samp.grid)){
    if(plot.type=="surface")
    {plot3D::persp3D(x = time.sub, y = freq.sub, z = t(amp.sub),
                     border = "black", lwd = lwd, theta = rotate.Xaxis,
                     phi = rotate.Yaxis, resfac = resfac, r = 3, expand = 0.5,
                     cex.axis = cex.axis, cex.lab=cex.lab, cex.main=cex.main,
                     scale = T, axes = T, col = color,
                     ticktype = "detailed", nticks = 4, xlab = "Time (s)",
                     ylab = "Frequency (kHz)", zlab = "Amplitude (dB)",
                     main = main, clab = expression("Amplitude (dB)"),
                     colkey=colkey)
    } # end plot.type=="surface"

    if(plot.type=="points") {
      plot3D::scatter3D(x = ind.3D[, 1], y = ind.3D[, 2],
                       z = ind.3D[, 3], pch = 21, cex = cex, theta = rotate.Xaxis,                          phi = rotate.Yaxis, resfac = resfac, r = 3, expand = 0.5,
                       cex.axis = cex.axis, cex.lab=cex.lab, cex.main=cex.main,
                       scale = T, axes = T, col = color,
                       ticktype = "detailed", nticks = 4, xlab = "Time (s)",
                       ylab = "Frequency (kHz)", zlab = "Amplitude (dB)",
                       main = main, clab = expression("Amplitude (dB)"),
                       colkey=colkey)
    } # end plot.type=="points"


  } # end samp.grid==TRUE

  if(plot.exp==TRUE){grDevices::dev.off()} # export plot only if plot.exp==TRUE

} # end function
p-rocha/SoundShape documentation built on Aug. 29, 2023, 10:57 p.m.