R/rasterFlow.R

In ssaxe-usgs/prettyWater: Pretty Plots for Attractive Science

#' Raster Plot of Daily Streamflow.
#'
#' Downloads and plots all daily streamflow for a single USGS NWIS streamgage ID.
#' @param site_id USGS NWIS streamgage ID, 8 digits long in character format.
#' @param startDate Start date in YYYY-MM-DD format.  Leave as '' if all time period required.
#' @param endDate End date in YYYY-MM-DD format.  Leave as '' if all time period required.
#' @param main Main title of plot.  If left as NULL, site_id used.
#' @param pCode Parameter code.  Defaults to streamflow.
#' @param dataRaw Column name of value in downloaded raw datafile.
#' @param yearLines Should years be delineated?
#' @param contours Should contour lines be drawn?
#' @param contourQuantiles At what percentiles should contour lines be drawn?
#' @param verbose Should progress be written to console?
#' @param returnData Should downloaded data be returned in raw df format? If so, plot is suppressed.
#' @return Raster/heatmap of streamflow and downloaded raw data.
#' @export rasterFlow
#' @examples
#' rasterFlow()

rasterFlow <- function(site_id          = '04085427',
                       startDate        = '',
                       endDate          = '',
                       main             = NULL,
                       pCode            = '00060',
                       dataRaw          = 'X_00060_00003',
                       yearLines        = F,
                       contours         = F,
                       contourQuantiles = c(0.60, 0.95),
                       verbose          = F,
                       returnData       = F){
  require('dataRetrieval')
  require('lubridate')
  require('dplyr')
  require('tidyr')
  require('RColorBrewer')
  require('fields')
  
  # Import data
  if ( verbose ) cat( 'Downloading data\n' )
  raw <- dataRetrieval::readNWISdv( site_id[ 1 ], pCode, startDate, endDate )
  if (verbose) cat('Download complete\n')
  if (returnData) return(raw)
  
  # Make sure dataRaw argument is a column in downloaded dataset
  if ( !( dataRaw %in% colnames( raw ) ) ){
    cat( "'dataRaw' not available as a column name.  Downloaded dataset is returned, plotting will not occur.\n" )
    if (returnData){
      return( raw )
    }else{
      return(colnames(raw))
    }
  }
  # rename data column as 'values'
  colnames( raw )[ which( colnames( raw ) == dataRaw ) ] <- 'Values'
  
  # Subset to dates and data
  formal <- raw %>%
    dplyr::select( Date, Values ) %>%
    dplyr::mutate( DoY   = lubridate::yday(  Date ) ) %>%
    dplyr::mutate( Year  = lubridate::year(  Date ) ) %>%
    dplyr::mutate( Month = lubridate::month( Date ) ) %>%
    dplyr::mutate( WaterYear = Year )
  formal$waterYear[ formal$Month > 10 ] <- formal$waterYear[ formal$Month > 10 ] + 1
  formal <- formal %>%
    dplyr::select( WaterYear, DoY, Values )
  
  # convert to wide format (matrix)
  mat <- tidyr::spread( data  = formal,
                        key   = DoY,
                        value = Values )
  rownames( mat ) <- mat$WaterYear
  mat <- mat %>%
    dplyr::select( -WaterYear ) %>%
    as.matrix( . )
  mat[ mat == -999999 ] <- NA
  
  # Remove leapday blanks
  mat <- mat[, -366]
  
  # Reorder matrix
  #mat <- mat[ , c( 244:366, 1:243 ) ]
  mat <- mat[ , c( 244:365, 1:243 ) ]
  
  
  
  ### Plotting ###
  # define color palette
  getPalette <- colorRampPalette( c( 'lightgrey', 'grey', 'orange', '#ef3b2c', 'forestgreen', 'yellow', 'blue' ),
                                  bias = 4.5 )

  # Format plot area
  if ( verbose ) cat( 'Plotting\n' )
  par( mar = c( 0, 0, 0, 2 ),
       oma = c( 3, 3, 3, 2) )
  
  # make plot
  fields::image.plot( (t( mat )), col = getPalette( 1000 ), axes = F, horizontal = F)
  
  # add box
  box( )
  
  # if contours, add
  if ( contours ){
    contour( x      = t( mat ),
             add    = T,
             levels = quantile( x     = mat,
                                probs = contourQuantiles,
                                na.rm = T ),
             drawlabels = T,
             labels = paste0(contourQuantiles * 100, '%')) 
  }
  
  # get plot X coordinates
  xRange <- par()$usr[ 1:2 ]
  xDif   <- xRange[ 2 ] - xRange[ 1 ]
  xStep  <- xDif / ncol( mat )
  xSeq.outer <- seq( from = xRange[ 1 ],
                     to   = xRange[ 2 ],
                     by   = xStep )
  xSeq.inner <- seq( from = xRange[ 1 ] + ( xStep / 2 ),
                     to   = xRange[ 2 ] - ( xStep / 2 ),
                     by   = xStep )
  # add lines for months
  monthLengths <- c( 0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 )
  monthLengths <- c( 0, 31, 30, 31, 31, 28, 31, 30, 31, 30, 31, 31, 30 )
  
  endMonths    <- cumsum( monthLengths )
  abline( v   = xSeq.outer[ endMonths[ -length( endMonths ) ] ],
          lty = 3 )
  # add month names
  midPt <- NULL
  for ( i in 2:length( endMonths ) ){
    midPt <- c( midPt, mean( c( endMonths[ i ], endMonths[ i-1 ] ) ) )
  }
  # add x labels
  mtext( text = month.abb[c(10:12, 1:9)],
         side  = 1,
         line  = 0.3,
         at    = xSeq.inner[ round( midPt ) ],
         las   = 2,
         font  = 2 )
  
  # get plot Y coordinates
  yRange <- par()$usr[ 3:4 ]
  yDif   <- yRange[ 2 ] - yRange[ 1 ]
  yStep  <- yDif / nrow( mat )
  ySeq.outer <- seq( from = yRange[ 1 ],
                     to   = yRange[ 2 ],
                     by   = yStep )
  ySeq.inner <- seq( from = yRange[ 1 ] + ( yStep / 2 ),
                     to   = yRange[ 2 ] - ( yStep / 2 ),
                     by   = yStep )
  if ( yearLines ) abline( h   = ySeq.outer,
                           lwd = 0.05 )
  
  # select y labels with "pretty"
  prettyLabels <- rownames( mat )[ base::pretty( x     = 1:nrow( mat ),
                                                 min.n = 2 ) 
                                   ]
  prettyLabels.locs <- match( x     = prettyLabels,
                              table = rownames( mat ) )
  
  mtext( text = prettyLabels,
         side = 2,
         line = 0.3,
         at   = ySeq.inner[ prettyLabels.locs ],
         las  = 1,
         font = 2)
  
  # Add main title
  if ( is.null( main ) ){
    ttext <- site_id
  }else{
    ttext <- main
  }
  mtext( text = paste( 'Gage No.', ttext ),
         font = 2,
         cex  = 1.2,
         line = 1.5 )
  
  mtext( text = paste0('Discharge (cfs), ', rownames( mat )[ 1 ], ' - ', rownames( mat )[ nrow( mat ) ] ),
         font = 2,
         cex  = 1.1,
         line = 0.5 )
  # done
  if ( verbose ) cat( 'Completed\n' )
  if (returnData){
    return( raw )
  }else{
    return( )
  }
}