R/EchogramCurtain.R

In dmwarn/EchoNet2FishnO: Estimate Fish Abundance from Acoustic Echoes and Net Catch

#' Create echogram curtain plot from interval-layer data.
#'
#' @param the.dir
#' Character, the directory containing the data that will be used to make the
#' plots.
#' @param d.units
#' Character, either "nph" or "gph".
#' @param spgrp
#' Character vector listing the species groups that should be plotted.
#' Must be entered in the form "X106.L100", for example, for alewife
#' (species code = 106) that are >= 100 mm long. These names are
#' in the form generated by the estimateLake() function.
#' @param QuantileFilter
#' Logical scalar indicating whether or not the user wants to plot all
#' values as they are or if the user wants to assign the values above
#' a particulare percentile the value of the chosen percentile. The intent
#' of this parameter is to allow one to choose to use the plots to find
#' acoustic cells that might have unreasonable density and need to be checked
#' or, instead to plot distribution after confirming that any really high
#' values appear to be real. This allows one to reduce the influence of a
#' few very high values on distribution patterns.
#' @param Percentile
#' Numeric scalar. The data values that are equal to or greater than
#' this percentile will be assigned the value of this percentile.
#' @return
#' This function generates curtain plots for each species group
#' listed and every acoustic transect. The transects are facets,
#' the vertical axis is depth, and the horizontal axis is transect
#' interval number.
#'
#' @details
#' This function is intended to build on the plots in exploreACMT(). The plots
#' provide a view of the data that is not possible when looking at just Sv
#' or TS/sigma, as these two can interact in such a way as to generate fish
#' density that might be higher than expected given the appearance of the Sv
#' data.One key use for this sort of plot is to look for analysis cells that
#' have density that is higher than expected because of low numbers of targets.
#' @export
#'
#' @examples
#'
#' @importFrom magrittr "%>%"
#' @importFrom lubridate today
#' @importFrom scales "pretty_breaks"
#' @import dplyr tidyr ggplot2
#'

EchogramCurtain <-
  function(the.dir,
           d.units,
           spgrp = c('X106.L0',
                     "X106.L100",
                     "X109.L0",
                     "X109.L90",
                     "X204.L0",
                     "X204.L120"), QuantileFilter = FALSE, Quantile = 0.99) {
    # The remaining code between here and ggplot() fills cells where there is a missing density estimate because there
    # were either no targets or Sv was below threshold. These cases really ought to be filled with a zero.
    # What follows is my jankity attempt to do that.
    #
    the.dir <- choose.dir()
    flist <- list.files(the.dir, pattern = paste0("intlaymeans_",d.units), full.names = TRUE)
    dat <- read.csv(flist)

    # Read in the data

    # 1) create a data.frame() that mean bottom depth for each transect and interval
    bott.parms <- dat %>% group_by(Region_name, Interval) %>%
      summarise(bott = mean(depth.botmin))

    # 2) Create a data.frame that has all possible combinations of Interval-Layer-Layer_depth_min
    # This will be used as the "full" combination of cells given the intervals and bottom depths in each transect
    lay.parms.min <-
      expand.grid(
        Interval = unique(dat$Interval),
        Layer = unique(dat$Layer),
        Layer_depth_min = unique(dat$Layer_depth_min)
      )

    lay.parms.min$d <- (lay.parms.min$Layer * 10) - 10

    # 3) We don't want to plot (or create) data for cells that DON'T exist on a given transect. For example,
    # if Layer_depth_min does not exceed 90 meters for a transect-interval, we don't want to plot data
    # below that depth.
    # We create a subset of the possible cells if all depths were present on all transects (lay.parms.min)
    # that represents the combination of Interval-Layer actually present on the transects
    lay.parms.present <- subset(lay.parms.min, d == Layer_depth_min)

    set1 <-
      merge(bott.parms,
            lay.parms.present,
            by = c("Interval"),
            sort = TRUE)

    dat.fill <-
      merge(
        set1,
        dat,
        by = c("Region_name", "Interval", "Layer", "Layer_depth_min"),
        all.x = TRUE
      )

    dat.fill$depth.botmin.new <-
      ifelse(is.na(dat.fill$depth.botmin),
             dat.fill$bott,
             dat.fill$depth.botmin)

    dat.fill2 <- subset(dat.fill, Layer_depth_min < depth.botmin.new)

    #define cols we want to keep
    keep.cols <-c(
        "Region_name",
        "Interval",
        "Layer_depth_min",
        'fish_ha',
         spgrp)

    put.zero.cols <-
      (c('fish_ha', spgrp))

    dat.fill3 <- dat.fill2[c(keep.cols)]
    dat.fill3[, put.zero.cols][is.na(dat.fill3[, put.zero.cols])] <- 0
    dat.fill3$fish_ha[dat.fill3$fish_ha == Inf |
                        is.na(dat.fill3$fish_ha)] <- 0

    df <- dat.fill3
    df.long <-
      gather(df, spgrp, density,-Region_name,-Layer_depth_min,-Interval, fish_ha) %>%
      filter(spgrp != 'fish_ha')
    # Create the list of spgrp to loop through
    spgrp.list <- unique(df.long$spgrp)
    leg.title <- ifelse(d.units == 'nph', "Fish per hectare", 'Grams per hectare')
    for (i in seq_along(spgrp.list)) {
      inp <- filter(df.long, spgrp == spgrp.list[i])
      Distance <- inp$Interval *3
      Depth <- inp$Layer_depth_min
      if(QuantileFilter == TRUE) {
        inp$density[inp$density >= quantile(inp$density, 0.99)] <- quantile(inp$density, 0.99)
      } else {
        inp$density <- inp$density
      }
      nam <- spgrp.list[i]
      p1 <-  ggplot() +
        geom_raster(
          data = inp,
          aes(
            x = as.factor(Distance),
            y = -Depth,
            fill = density),
          stat = "identity",
        interpolate = TRUE) +
        scale_fill_viridis_c(
          paste0(leg.title, "  ", spgrp.list[i]),
          option = "plasma",
          breaks = pretty_breaks(n = 10)
        ) +
        scale_x_discrete("Distance (km)", breaks = seq(0, max(Distance), 6)) +
        facet_wrap( ~ Region_name, scales = "free")
      ggsave(paste0(the.dir, "/", spgrp.list[i],"_", Sys.Date(), "_", d.units, "_echogram_curtain.png"), p1,
        dpi = 300,
        width = 12,
        height = 8)

    }
  }