R/gbm.map.R

In SimonDedman/gbm.auto: Automated Boosted Regression Tree Modelling and Mapping Suite

#' Maps of predicted abundance from Boosted Regression Tree modelling
#'
#' Generates maps from the outputs of gbm.step then Gbm.predict.grids, handled
#' automatically within gbm.auto but can be run alone, and generates
#' representativeness surfaces from the output of gbm.rsb.
#'
#' `r lifecycle::badge("superseded")`
#' Superseded by gbm.mapsf on 2023-08-07, but still works.
#'
#' @param x Vector of longitudes, from make.grid in mapplots; x. Order by this
#' (descending) SECOND.
#' @param y Vector of latitudes, from make.grid in mapplots; grids\[,gridslat\].
#' Order by this (descending) first.
#' @param z Vector of abundances generated by gbm.predict.grids, from make.grid
#' in mapplots; grids\[,predabund\].
#' @param byx Longitudinal width of grid cell, from make.grid in mapplots.
#' Autogenerated if left blank.
#' @param byy Latitudinal height of grid cell, from make.grid in mapplots.
#' Autogenerated if left blank.
#' @param grdfun make.grid operand for >=2 values per cell. Default:mean, other
#' options: sum prod min max sd se var.
#' @param mapmain Plot title, has species value appended. Default "Predicted
#' CPUE (numbers per hour): ".
#' @param species Response variable name, from basemap in mapplots;
#' names(samples\[i\]). Defaults to "Response Variable".
#' @param heatcolours Vector for abundance colour scale, defaults to the heatcol from
#' legend.grid and draw.grid in mapplots which is c("white", "yellow", "orange"
#' , "red", "brown4").
#' @param colournumber Number of colours to spread heatcol over, default:8.
#' @param shape Basemap shape to draw, from draw.shape in mapplots. Defaults to
#' NULL which calls gbm.basemap to generate it for you. First read in a shp file e.g. myshape <-
#' sf::st_read(dsn = paste0(savename, ".shp"), layer = savename, quiet = TRUE), then use shape = myshape.
#' @param landcol Colour for 'null' area of map (for marine plots, this is
#' land), from draw.shape in mapplots. Default "grey80" (light grey).
#' @param mapback Basemap background colour, defaults to lightblue (ocean for
#' marine plots).
#' @param legendloc Location on map of legend box, from legend.grid in mapplots,
#'  default bottomright.
#' @param legendtitle The metric of abundance, e.g. CPUE for fisheries, from
#' legend.grid in mapplots. Default "CPUE".
#' @param lejback Background colour of legend, from legend.grid in mapplots.
#' Default "white".
#' @param zero Force include 0-only bin in breaks.grid and thus legend? Default
#' TRUE.
#' @param quantile Set max quantile of data to include in bins, from breaks.grid
#'  in mapplots; lower to e.g. 0.975 cutoff outliers; default 1.
#' @param byxout Export byx to use elsewhere? Default:FALSE.
#' @param breaks Vector of breakpoints for colour scales; default blank,
#' generated automatically.
#' @param byxport Dummy param for package testing for CRAN, ignore.
#' @param ... Additional arguments for legend.grid's ... which passes to legend.
#'
#' @return Species abundance maps using data provided by gbm.auto, and
#' Representativeness Surface Builder maps using data provided by gbm.rsb, to be
#' run in a png/par/gbm.map/dev.off sequence.
#'
#' @details Errors and their origins:
#'
#' Error in seq.default(xlim\[1\], xlim\[2\], by = byx):wrong sign in 'by' argument
#' Check that your lat & long columns are the right way around. Ensure grids
#' data are gridded, i.e. they are in a regular pattern of same/similar lines of lat/lon, even if
#' they're missing sections.
#'
#' Suggested parameter values:
#' z = rsbdf\[,"Unrepresentativeness"\]
#'
#' mapmain = "Unrepresentativeness: "
#'
#' legendtitle = "UnRep 0-1"
#'
#' @export
#' @importFrom grDevices colorRampPalette dev.off grey.colors png
#' @importFrom graphics legend par
#' @importFrom mapplots basemap breaks.grid draw.grid draw.shape legend.grid make.grid
#' @author Simon Dedman, \email{simondedman@@gmail.com}
#' @author Hans Gerritsen
#' @examples
#' \donttest{
#' # Not run: downloads and saves external data.
#' # Suggested code for outputting to png:
#' data(grids)
#' # set working directory somewhere suitable
#' png(filename = "gbmmap.png", width = 7680, height = 7680, units = "px",
#' pointsize = 192, bg = "white", res = NA, family = "", type = "cairo-png")
#' par(mar = c(3.2,3,1.3,0), las = 1, mgp = c(2.1,0.5,0), xpd = FALSE)
#' gbm.map(x = grids[,"Longitude"], y = grids[,"Latitude"], z = grids[,"Effort"]
#' , species = "Effort")
#' dev.off()
#' }
#'
gbm.map <- function(x,        #vector of longitudes, from make.grid in mapplots; x. ##Order by this (descending) SECOND
                    y,        #vector of latitudes, from make.grid in mapplots; grids[,gridslat]. ##Order by this (descending) first
                    z,        #vector of abundances generated by gbm.predict.grids, from make.grid in mapplots; grids[,predabund]
                    byx = NULL,      #longitudinal width of grid cell, from make.grid in mapplots. Auto-generated if left blank.
                    byy = NULL,      #latitudinal height of grid cell, from make.grid in mapplots. Auto-generated if left blank.
                    grdfun = mean, #make.grid operand for >=2 values per cell. Default:mean, other options: sum prod min max sd se var.
                    mapmain = "Predicted CPUE (numbers per hour): ",  #plot title, has species value appended. Default "Predicted CPUE (numbers per hour): "
                    species = "Response Variable",  #Response variable name, from basemap in mapplots; names(samples[i]). Defaults to "Response Variable"
                    heatcolours = c("white", "yellow", "orange","red", "brown4"), # Vector for abundance colour scale, defaults to the heatcol from legend.grid & draw.grid in mapplots.
                    colournumber = 8,   #number of colours to spread heatcol over, default:8
                    shape = NULL,   #basemap shape to draw, from draw.shape in mapplots. Defaults to NULL which calls gbm.basemap to generate it for you
                    landcol = "grey80", #colour for 'null' area of map (for marine plots, this is land), from draw.shape in mapplots. Default "grey80" (light grey)
                    mapback = "lightblue", #basemap background colour, defaults to lightblue (ocean for marine plots)
                    legendloc = "bottomright", #location on map of legend box, from legend.grid in mapplots, default bottomright
                    legendtitle = "CPUE", #the metric of abundance, e.g. CPUE for fisheries, from legend.grid in mapplots. Default "CPUE"
                    lejback = "white",  #background colour of legend, from legend.grid in mapplots. Default "white"
                    zero = TRUE, # force include 0-only bin in breaks.grid and thus legend? Default TRUE
                    quantile = 1, # set max quantile of data to include in bins, from breaks.grid in mapplots; lower to e.g. 0.975 cut-off outliers; default 1
                    byxout = FALSE, # export byx to use elsewhere? Default:FALSE
                    breaks = NULL, # vector of breakpoints for colour scales; default blank, generated automatically
                    byxport = NULL, # addresses devtools::check's no visible binding for global variable https://cran.r-project.org/web/packages/data.table/vignettes/datatable-importing.html#globals
                    ...) # additional arguments for legend.grid's ... which passes to legend, and savedir and others for gbm.basemap

# Generalised Boosting Models, automated map generator. Simon Dedman, 2014-2016, simondedman@gmail.com, https://github.com/SimonDedman/gbm.auto

# Generates maps from the outputs of gbm.step then gbm.predict.grids, handled automatically within gbm.auto but can be run alone, and
# generates representativeness surfaces from the output of gbm.rsb (suggest: z = rsbdf[,"Unrepresentativeness"],
# mapmain = "Unrepresentativeness: ",legendtitle = "UnRep 0-1"). Suggested code for outputting to e.g. png:
# png(...); par(...); gbm.map(...); dev.off()
{
  # utils::globalVariables("byxport") # addresses devtools::check's no visible binding for global variable https://cran.r-project.org/web/packages/data.table/vignettes/datatable-importing.html#globals

  # require(mapplots)
  if (is.null(shape)) { # if no map shape entered, generate bounds and call gbm.basemap
    bounds = c(range(x),range(y))
    #create standard bounds from data, and extra bounds for map aesthetic
    shape <- gbm.basemap(bounds = bounds, extrabounds = TRUE)
  }

  # if user hasn't entered byx or byy values, generate them from the data
  if (is.null(byx)) {
    # work out cell size for uniform square gridded data: Create blank vector for grid length calcs
    bydist <- rep(NA, length(x))
    # and attach it to grids
    cells <- data.frame(LONGITUDE = x, bydist = bydist, stringsAsFactors = FALSE)
    # fill it: if [next longitude minus current longitude] equals [current longitude minus previous longitude], that's a uniform cell.
    # data rounded to prevent tiny fluctuations counting as differences. Need to set that tolerance.
    # Could do 10% of average distance between uniform points, but you don't know that til the end!
    cells[2:(length(x) - 1),"bydist"] <-
      ifelse(round(cells[2:(length(x) - 1),1] - cells[1:(length(x) - 2),1], digits = 5)
             ==
               round(cells[3:length(x),1] - cells[2:(length(x) - 1),1], digits = 5),
             round(cells[2:(length(x) - 1),1] - cells[1:(length(x) - 2),1], digits = 5),
             NA)
    # Take an average of those distances, they should all be identical anyway. Apply it to byx & byy.
    byx <- mean(cells$bydist, na.rm = TRUE)
    # if users grids data are organised in COLUMNS rather than rows, they'll have long blocks of the same longitude
    # the mean distance will be 0, and make.grid will fail when it calls seq(xlim[1], xlim[2], by = byx) and byx == 0
    # in case of this, try latitude instead of longitude
    if (byx == 0) {
      bydist <- rep(NA, length(y))
      cells <- data.frame(LONGITUDE = y, bydist = bydist, stringsAsFactors = FALSE)
      cells[2:(length(y) - 1),"bydist"] <-
        ifelse(round(cells[2:(length(y) - 1),1] - cells[1:(length(y) - 2),1], digits = 5) ==
                 round(cells[3:length(y),1] - cells[2:(length(y) - 1),1], digits = 5),
               round(cells[2:(length(y) - 1),1] - cells[1:(length(y) - 2),1], digits = 5),
               NA)
      byx <- abs(mean(cells$bydist, na.rm = TRUE))
    }
    byy <- byx
    if (byxout) byxport <<- byx
  }
  grd <- mapplots::make.grid(x, y, z, byx, byy, xlim = range(x), ylim = range(y), fun = grdfun)
  # create gridded data. fun defaults to sum which is bad since accidentally
  # overlapping points will look anomalous
  heatcol = colorRampPalette(heatcolours)(colournumber)
  # create heatcol from component parts
  if (is.null(breaks)) breaks <- mapplots::breaks.grid(grd, zero = zero, quantile = quantile, ncol = length(heatcol))
  # if breaks specified, do nothing (used later in draw.grid) else generate them
  if (zero) {heatcol = c("#00000000", colorRampPalette(heatcol)(length(heatcol) - 1))}
  # if zero = TRUE add alpha as 1st colour (1st 2 breakpoints)
  mapplots::basemap(xlim = range(x), ylim = range(y), main = paste0(mapmain, species), bg = mapback, xlab = "Longitude", ylab = "Latitude")
  mapplots::draw.grid(grd, breaks, col = heatcol)
  # plot grd data w/ breaks for colour breakpoints
  mapplots::draw.shape(shape = shape, col = landcol) # add coastline
  mapplots::legend.grid(legendloc, breaks = breaks, type = 2, inset = 0, bg = lejback, title = legendtitle, col = heatcol)
  # breaks=breaks/1000 was causing odd legend? From make.grid help, Hans using to convert kg to t?
  # removed the '...' from legend.grid, preventing optional args from being passed to legend.
  # Legend doesn't have a '...' so ALL args floating in the environment will be passed
  # and expected to be used. Like max.trees and other gbm args. This kills the crab.
  # Could potentially write out the whole of legend.grid (it's small) then call
  # legend with R.utils::doCall("legend" x, y, legend = legend, col = col[ncol:1], pch = pch,
  ## pt.cex = pt.cex, bg = bg, ..., .ignoreUnusedArgs=TRUE))
}