R/convert_to_huc12.R

In ssaxe-usgs/hru2huc: Convert HRU Values to HUC 2-12

#' Convert HRU Data to HUC 12
#'
#' Converts time-series HRU data to overlying HUC 12 regions. Returns a zoo object.
#' @param hru.data HRU time series data in matrix format where rows are time steps and columns are HRUs.
#' 109,951 columns required.
#' @param na.cutoff Threshold value to decide when to return an NA value for a HUC 12 region.  Default is 0.5.  In
#' this case, if more than 50 percent of a HUC 12 region is NA, NA is returned.
#' @param dateVector Vector of dates in "YYYY-MM-DD" format to be added.
#' @return HUC12 matrix of converted values.
#' @export convert_to_huc12
#' @examples
#' convert_to_huc12()

convert_to_huc12 <- function(hru.data,
                             na.cutoff    = 0.5,
                             dateVector   = NULL){
  # require 109951 columns
  if (ncol(hru.data) != 109951) stop("Data must have 109,951 columns.")
  # load weights data
  errText <-"'h2hData' package is not installed.  Please install with devtools::install_github(repo = 'ssaxe-usgs/h2hData')."
  if (("h2hData" %in% rownames(installed.packages())) == F) stop(errText)
  data("hru_to_HUC12", package = "h2hData", envir = environment())

  # define function
  calcHUC <- function(wts, hru.values, na.cutoff){
    # reorder weights
    wts     <- wts[ order( wts$HRU ), ]
    wts$HRU <- as.character( wts$HRU )
    wts$HRU[ wts$HRU == "1e+05" ] <- '100000'
    # if more than 50% of the area is NA, return vector of NAs.
    # otherwise, distribute weights to other HRUs evenly
    if ( anyNA( wts$HRU ) ){
      na.weight <- wts$weight[ is.na( wts$HRU ) ]
      if ( na.weight >= na.cutoff ){
        # return NA
        return( rep( NA, nrow( hru.values ) ) )
      }else{
        # distribute NA weights
        pre.Weight          <- wts$weight[ !is.na( wts$HRU ) ]
        add.weights         <- ( pre.Weight / sum( pre.Weight ) ) * na.weight
        post.Weight         <- pre.Weight + add.weights
        wts.official        <- wts[ !is.na( wts$HRU ), ]
        wts.official$weight <- post.Weight
      }
    }else{
      wts.official          <- wts
    }

    # isolate relevant hrus
    col.select <- as.numeric( wts.official$HRU )
    #rawValues <- hru.values[ , col.select + nPass]
    rawValues <- hru.values[ , col.select]

    if ( !is.vector( rawValues ) ){
      # order correctly
      rawValues <- as.matrix( rawValues[ , order( wts.official$HRU ) ] )
      # apply weights
      wtedValues <- ( rawValues %*% wts.official$weight )[ , 1 ]
      # return
      return( wtedValues )
    }else{
      # return
      return( rawValues )
    }
  }

  # apply in loop
  convertedValues <- pbapply::pblapply(X            = hru_to_HUC12,
                                       FUN          = calcHUC,
                                       hru.values   = hru.data,
                                       na.cutoff    = na.cutoff)
  # combine into matrix
  convertedValues <- do.call(cbind, convertedValues)
  # clear excess
  rm(hru_to_HUC12, envir = environment())
  gc()

  # add dates as rownames
  if (!is.null(dateVector)){
    rownames(convertedValues) <- dateVector
  }

  # sum extension names
  extNames <- colnames(convertedValues)[nchar(colnames(convertedValues)) > 12]
  extData  <- convertedValues[, nchar(colnames(convertedValues)) > 12]
  shortData <- convertedValues[, nchar(colnames(convertedValues)) == 12]
  subNames <- substr(extNames, 1, 12)
  summedExt <- t(aggregate(x = t(extData),
                         by = list(subNames),
                         FUN = sum)[,-1])
  colnames(summedExt) <- sort(unique(subNames))
  allData <- cbind(summedExt, shortData)
  allData <- allData[, order(colnames(allData))]

  # return
  return(allData)
}