R/convert_to_nwis.R

In ssaxe-usgs/hru2nwis:

#' Convert HRU Data to NWIS
#'
#' Converts time-series HRU data to overlying NWIS 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 NWIS region.  Default is 0.5.  In
#' this case, if more than 50 percent of a NWIS region is NA, NA is returned.
#' @param dateVector Vector of dates in "YYYY-MM-DD" format to be added.
#' @return NWIS matrix of converted values.
#' @export convert_to_nwis
#' @examples
#' convert_to_nwis()

convert_to_nwis <- 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
  hru_to_nwis <- hru2nwis::hru_to_nwis

  # define function
  calcNWIS <- 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_nwis,
                                        FUN          = calcNWIS,
                                        hru.values   = hru.data,
                                        na.cutoff    = na.cutoff )
  # combine into matrix
  convertedValues <- do.call( cbind, convertedValues )
  # clear excess
  rm( hru_to_nwis, envir = environment( ) )
  gc( )

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

  # order data
  allData <- convertedValues[ , order( colnames( convertedValues ) ) ]

  # return
  return( allData )
}