MazamaSpatialPlots: Thematic Plots for Mazama Spatial Datasets

Documented in stateMap

#' @title State level thematic map
#' @description Uses the \pkg{tmap} package to generate a thematic map at the
#' state level. Input consists of a dataframe with \code{stateCode} identifiers.
#'
#' Data to plot is specified with \code{parameter} argument. If \code{parameter}
#' is mult-valued, mutliple plots will be generated and displayed as "facets".
#'
#' The returned object is a \pkg{tmap} ggplot object which can be further
#' modified with tmap or ggplot options.
#'
#' @param data Dataframe containing values to plot. This dataframe
#' must contain a column named \code{stateCode} with the 2-character state code.
#' @param parameter Name of the column in \code{data} to use for coloring the map.
#' @param state_SFDF simple features data frame with US states. It's data
#' \code{@slot} must contain a column named \code{stateCode} with the
#' 2-character state code.
#' @param breaks Numeric vector of break points.
#' @param palette A vector of colors or palette name from the \pkg{cols4all} package
#' (see \code{\link[cols4all:c4a]{cols4all::c4a}}).
#' @param conusOnly Logical specifying Continental US state codes. Ignored when
#' the \code{stateCode} argument is specified.
#' @param stateCode Vector of state codes to include on the map.
#' @param projection Named projection, \emph{e.g.} "EPSG:4326" or "WGS84" or proj4string.
#' @param stateBorderColor Color used for state borders.
#' @param title Vector of text strings to use as individual plot titles.
#' This must be the same length as 'parameter'.
#' @param showLegend Logical specifying whether or not to show the legend.
#' @param legendTitle Text string to use as the legend title.
#' @param legendOrientation Orientation of the legend. Either "portrait" or "landscape".
#' @param legendPosition A \emph{tm_pos} object generated with
#' \code{\link[tmap:tm_pos_in]{tmap::tm_pos_in()}} or
#' \code{\link[tmap:tm_pos_out]{tmap::tm_pos_out()}}.
#'
#' @return A ggplot object.
#'
#' @rdname stateMap
#'
#' @examples
#' library(MazamaSpatialPlots)
#'
#' stateMap(
#'   data = example_US_stateObesity,
#'   parameter = "obesityRate",
#'   stateBorderColor = "white",
#'   title = "2018 Obesity by State"
#' )
#'
#' # Example of customization using tm_layout and breaks parameter
#' stateMap(
#'   data = example_US_stateObesity,
#'   parameter = "obesityRate",
#'   breaks = seq(20,38,3),
#'   stateBorderColor = 'black'
#' ) +
#'   tmap::tm_layout(
#'     frame = TRUE,
#'     frame.double.line = TRUE,
#'     bg.color = "grey85",
#'     inner.margins  = .05
#'   ) +
#'   tmap::tm_title(
#'     text = 'Obesity Rate by State',
#'     size = 1.5,
#'     position = tmap::tm_pos_in("center", "top"),
#'     fontfamily = "serif"
#'   )
#'
#' # Example using stateCode
#' stateMap(
#'   data = example_US_stateObesity,
#'   parameter = "obesityRate",
#'   stateCode = c('ME', 'NH', 'VT', 'MA', 'RI', 'CT'),
#'   stateBorderColor = 'black',
#'   legendPosition = tmap::tm_pos_in("right", "bottom")
#' ) +
#'   tmap::tm_layout(
#'     frame = TRUE,
#'     frame.double.line = TRUE,
#'     bg.color = "grey85",
#'     inner.margins  = .08
#'   ) +
#'   tmap::tm_title(
#'     text = 'Obesity Rates in New England',
#'     size = 1.5,
#'     fontface = 2,
#'     fontfamily = "serif",
#'     position = tmap::tm_pos_in("center", "top")
#'   )
#' @export
#' @importFrom sf st_crs st_bbox
#' @importFrom rlang .data
#' @importFrom tmap tm_shape tm_fill tm_polygons tm_layout
#'

stateMap <- function(
    data = NULL,
    parameter = NULL,
    state_SFDF = "USCensusStates_02",
    breaks = NULL,
    palette = "brewer.blues",
    conusOnly = TRUE,
    stateCode = NULL,
    projection = NULL,
    stateBorderColor = "gray50",
    title = NULL,
    showLegend = TRUE,
    legendTitle = NULL,
    legendOrientation = "portrait",
    legendPosition = NULL
) {

  # ----- Validate parameters --------------------------------------------------

  MazamaCoreUtils::stopIfNull(data)
  MazamaCoreUtils::stopIfNull(parameter)

  # Check if data exists
  if ( !exists("data") ) {
    stop("Parameter 'data' is empty.")
  }

  # * Validate 'data' -----

  # Does 'parameter' exist in the incoming dataframe?
  if ( !all(parameter %in% colnames(data)) ) {
    stop(sprintf("'%s' is not found in the incoming 'data' dataframe.", parameter))
  }

  # Does 'data' have the required columns?
  requiredFields <- c("stateCode")
  missingFields <- setdiff(requiredFields, names(data))
  if ( length(missingFields) > 0 ) {
    if ( "state name" %in% tolower(names(data)) ) {
      columnIndex <- which("state name" == tolower(names(data)))
      data$stateCode <- MazamaSpatialUtils::US_stateNameToCode(data[, columnIndex])
    } else {
      stop(paste0("Missing fields in 'data': ",
                  paste0(missingFields, collapse = ", ")))
    }
  }

  # * Validate 'SFDF' -----

  # Accept state SFDF as character string or as object
  if ( is.character(state_SFDF) ) {
    if ( exists(state_SFDF) ) {
      state_SFDF <- get(state_SFDF)
    } else {
      stop(sprintf("State dataset '%s' is not loaded.
  Please load it with MazamaSpatialtUtils::loadSpatialData()",
                   state_SFDF
      ))
    }
  }

  # Does 'state_SFDF' have the required columns?
  requiredSFDFFields <- c("stateCode")
  missingSFDFFields <- setdiff(requiredSFDFFields, names(state_SFDF))
  if ( length(missingSFDFFields) > 0 ) {
    if ( "state name" %in% tolower(names(state_SFDF)) ) {  # if statecode not found, try to create it
      columnIndex <- which("state name" == tolower(names(state_SFDF)))
      state_SFDF$stateCode <- MazamaSpatialUtils::US_stateNameToCode(state_SFDF[, columnIndex])
    } else {
      stop(paste0("Missing fields in 'state_SFDF': ",
                  paste0(missingSFDFFields, collapse = ", ")))
    }
  }

  # * Validate other -----

  # Validate breaks
  if ( !is.null(breaks) && !is.numeric(breaks) ) {
    stop("Parameter 'breaks' must be a numeric vector.")
  }

  # Guarantee that conusOnly is logical
  if ( !is.logical(conusOnly) ) {
    conusOnly <- TRUE
  }

  # Validate incoming stateCode vector if any
  if ( !is.null(stateCode) ) {
    if ( !all(stateCode %in% MazamaSpatialUtils::US_52) ) {
      invalidStateCodes <- setdiff(stateCode, MazamaSpatialUtils::US_52)
      stop(sprintf("Invalid state codes found:  %s",
                   paste0(invalidStateCodes, collapse = ", ")))
    }
  }

  # Convert projection to a CRS object if necessary
  if ( !is.null(projection) ) {
    if ( is.character(projection) ) {
      projection <- sf::st_crs(projection)
    } else if ( "CRS" %in% class(projection) ) {
      # leave it alone
    } else {
      stop(paste0("Parameter 'projection' must be a sf::st_crs object or a valid ",
                  "projection string."))
    }
  }

  # Validate length of 'parameter' and 'title'
  if ( !is.null(title) && length(parameter) != length(title) )
    stop("The lengths of 'parameter' and 'title' must be equal.")

  # Validate legendPosition
  if ( !is.null(legendPosition) ) {
    if ( !"tm_pos" %in% class(legendPosition) ) {
      stop("Parameter 'legendPosition' must be generated with tmap::tmap_pos_in() or tmap::tmap_pos_out()")
    }
  }

  # ----- Subset the SFDF ------------------------------------------------------

  if ( !is.null(stateCode) ) {

    # NOTE:  Subset doesn't work properly unless we change the name here
    incomingStateCode <- stateCode
    state_SFDF <- subset(state_SFDF, state_SFDF$stateCode %in% incomingStateCode)
    data <- data %>% dplyr::filter(.data$stateCode %in% incomingStateCode)

  } else if ( conusOnly ) {

    state_SFDF <- subset(state_SFDF, state_SFDF$stateCode %in% MazamaSpatialUtils::CONUS)
    data <- data %>% dplyr::filter(.data$stateCode %in% MazamaSpatialUtils::CONUS)

  } else {

    # use existing SFDF

  }

  # ----- Create a projection --------------------------------------------------

  # NOTE:  Setting the est most longitude improves plots that include Alaska by
  # NOTE:  assuming the east most longitude cannot be greater than -66.97626
  # NOTE:  (which is the farthest east longitude in the continental US).
  # NOTE:  This wouldn't work if someone wants to include Guam but this function
  # NOTE:  is designed primarily for CONUS.

  if ( is.null(projection) ) {

    # NOTE:  Specifying stateCode takes precedence over specifying conusOnly
    if ( !is.null(stateCode) ) {

      # 1) Get boundaries from stateSFDF
      bbox <- sf::st_bbox(state_SFDF)

      # 2) Calculate lat lo/mid/hi and lon mid
      lat_1 <- bbox[2]
      lat_2 <- bbox[4]
      lat_0 <- (lat_1 + lat_2)/2
      lon_1 <- bbox[1]
      lon_2 <- min(bbox[3], -66.97626) # this is the east most longitude to use
      lon_0 <- (lon_1 + lon_2)/2

      # 3) Create the proj4string text
      projString <- sprintf("+proj=aea +lat_1=%.1f +lat_2=%.1f +lat_0=%.1f +lon_0=%.1f +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs",
                            lat_1, lat_2, lat_0, lon_0)
      projection <- sf::st_crs(projString)

    } else if ( conusOnly ) {

      projection <- sf::st_crs("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")

    } else {

      # TODO:  The 'conusOnly = FALSE' option needs to be replaced with the
      # TODO:  artifical projection that places Alaska and Hawaii in the lower
      # TODO:  left corner. Jon has seen such a projection in an example but
      # TODO:  can't remember exactly where.

      projection <- sf::st_crs("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")

    }

  } else {

    # TODO:  When do we use the passed in projection?
    # Use projection found in SFDF

  }

  # ----- Merge data with SFDF -------------------------------------------------

  # NOTE:  We can use left_join() because 'stateCode' is guaranteed to be in
  # NOTE:  both dataframes. We use "matrix style" subsetting of 'data' to
  # NOTE:  specify that we want all rows and just the columns with "stateCode"
  # NOTE:  and the parameter of interest.

  # Add the incoming data$parameter to 'state_SFDF'.
  state_SFDF <-
    dplyr::left_join(
      state_SFDF,
      data[, c("stateCode", parameter)],
      by = "stateCode"
    )

  # ----- Create plot ----------------------------------------------------------

  gg <-
    tmap::tm_shape(state_SFDF, crs = projection) +
    tmap::tm_fill(
      fill = parameter,
      fill.scale = tmap::tm_scale_intervals(
        breaks = breaks,
        values = palette
      ),
      fill.legend = tmap::tm_legend(
        title = legendTitle,
        show = showLegend,
        orientation = legendOrientation,
        position = legendPosition
      )
    ) +
    tmap::tm_polygons(
      fill_alpha = 0,
      col = stateBorderColor
    ) +
    tmap::tm_title(
      text = title,
      size = 0.9,
      position = tmap::tm_pos_out("center", "top")
    )

  # ----- Return ---------------------------------------------------------------

  return(gg)

}

# ===== DEBUGGING ==============================================================

if ( FALSE ) {

  # NOTE:  Good idea to clear out the Global environment before running this.

  library(MazamaSpatialPlots)
  mazama_initialize()

  state_SFDF <- USCensusStates_02

  # Set up required variables so we can walk through the code
  data = example_US_stateObesity
  parameter = "obesityRate"
  breaks = NULL
  palette = "brewer.blues"
  conusOnly = TRUE
  stateCode = NULL
  projection = NULL
  stateBorderColor = "white"
  title = "Obesity Rate by state"
  showLegend = TRUE
  legendTitle = NULL
  legendOrientation = "portrait"
  legendPosition = NULL

  # Run the code above and then start walking through the lines of code in the
  # function.
  #
  # You can also source this file and run the function now

  stateMap(
    data = data,
    parameter = parameter,
    state_SFDF = state_SFDF,
    breaks = breaks,
    palette = palette,
    conusOnly = conusOnly,
    stateCode = stateCode,
    projection = projection,
    stateBorderColor = stateBorderColor,
    title = title,
    legendPosition = legendPosition
  )

  ##############################################################################
  # Here is how I would use this function in real life:

  stateMap(
    data,
    state_SFDF = state_SFDF,
    parameter = "obesityRate"
  )

  # Not bad

  stateMap(
    data,
    state_SFDF = state_SFDF,
    parameter = "obesityRate",
    stateBorderColor = "black"
  )

  # Better but still need a title

  stateMap(
    data,
    state_SFDF = state_SFDF,
    parameter = "obesityRate",
    stateBorderColor = "black",
    legendPosition = tmap::tm_pos_in("left", "bottom")
  ) +
    tmap::tm_layout(
      frame = TRUE
    ) +
    tmap::tm_title(
      text = "Obesity rate by state",
      size = 2,
      fontface = "bold"
    )

  # Very nice!

  ##############################################################################
  # More examples

  # Example using tmap style "classic"
  stateMap(
    data = example_US_stateObesity,
    parameter = "obesityRate",
    breaks = seq(20,40,4),
  ) +
    tmap::tm_style(
      'classic'
    ) +
    tmap::tm_title(
      text = 'Obesity Rate by State',
      position = tmap::tm_pos_in("center", "top"),
      size = 1.2
    ) +
    tmap::tm_compass()

  # Example using tmap style "cobalt"
  stateMap(
    data = example_US_stateObesity,
    parameter = "obesityRate",
    breaks = seq(20,40,4)
  ) +
    tmap::tm_style(
      'cobalt'
    ) +
    tmap::tm_title(
      text = 'Obesity Rate by State',
      position = tmap::tm_pos_in("center", "top"),
      size = 1.2
    ) +
    tmap::tm_compass()

  # the previous two examples show how to leverage a preexisting style to create
  # nice plots with few manual specifications

  # Example with conusOnly=FALSE
  stateMap(
    data = example_US_stateObesity,
    parameter = "obesityRate",
    breaks = seq(20,38,3), #increasing color detail
    conusOnly = FALSE ,
    stateBorderColor = 'black',
    legendPosition = tmap::tm_pos_in("right", "top")
  ) +
    tmap::tm_title(
      text = 'Obesity Rates in U.S. States and Territories',
      position = tmap::tm_pos_out("center", "top"),
      size = 2,
      fontface = 2,
      fontfamily = "serif"
    )

}