R/downscale_HUC12.R

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

#' Downscale HUC 12 to HUCs 2-10
#'
#' Converts HUC 12 time-series data to HUC levels 2-10 based on area weighting.
#' @param HUC12_matrix Time-series matrix of HUC 12 values produced through convert_to_huc() function.
#' @param na.cutoff Threshold value to decide when to return an NA value for a region.  Default is 0.5.  In
#' this case, if more than 50 percent of a region is NA, NA is returned.
#' @param verbose Should progress be written to console?
#' @keywords huc, HUC
#' @export downscale_HUC12
#' @examples
#' downscale_HUC12()

downscale_HUC12 <- function(HUC12_matrix,
                            na.cutoff    = 0.5,
                            verbose      = T){
  require(magrittr)

  ###### Define function for applying weights
  # weighting function
  applyWeights <- function( lowerRes_name,
                            conversionList,
                            higherRes_data,
                            na.cutoff ){

    # subset list to current lower res HUC
    conversionList_subset <- conversionList[[ which( names( conversionList ) == lowerRes_name ) ]]
    # remove subset letters (example at NOI[3203])
    conversionList_subset$HUC12 <- substr( conversionList_subset$HUC12, 1, 12 )
    # if repeats, sum weights
    conversionList_subset <- aggregate( conversionList_subset$weights, list( conversionList_subset[ , 1 ] ), sum )
    colnames( conversionList_subset ) <- c( 'HUC12', 'weights' )

    # subset higher res matrix by rownames from conversionList_subset
    higherRes_data_subset <- higherRes_data[ , colnames( higherRes_data ) %in% conversionList_subset[ , 1 ] ]
    # mmm <- conversionList_subset[(conversionList_subset$HUC12 %in% colnames(higherRes_data)) == F, ]
    # colnames(higherRes_data)[nchar(colnames(higherRes_data)) > 12]

    # apply weights

    if ( !is.null( dim( higherRes_data_subset ) ) ){
      # reorder columns
      higherRes_data_subset <- higherRes_data_subset[ , order( colnames( higherRes_data_subset ) ) ]
      # reorder weights
      conversionList_subset <- conversionList_subset[ order( conversionList_subset[ , 1 ] ), ]
      # weight
      #higherRes_data_subset_weighted <- higherRes_data_subset %*% conversionList_subset$weights

      # if entire dataset is NA, return vector of NAs. else, weight
      if ( all( is.na( higherRes_data_subset ) ) ){
        if ( is.vector( higherRes_data_subset ) ){
          return(
            rep( NA, length( higherRes_data_subset ) )
          )
        }else{
          return(
            rep( NA, nrow( higherRes_data_subset ) )
          )
        }
      }
      # weight via apply
      weightFun <- function(rawData, wts, na.cutoff){
        wts.na.sum <- sum( wts[ is.na( rawData ) ] )
        if ( wts.na.sum > na.cutoff ){
          #higherRes_data_subset_weighted <- rep(NA, length(rawData))
          return(
            #rep(NA, length(rawData))
            NA
          )
        }else{
          wts.new <- wts[ !is.na( rawData ) ] / sum( wts[ !is.na( rawData ) ] )
          rawData.new <- rawData[ !is.na( rawData ) ]
          retValue <- tryCatch({
            sum( wts.new * rawData.new )
          }, warning = function( w ){ -9999 })
          return( retValue )
        }
      }
      higherRes_data_subset_weighted <- apply( X         = higherRes_data_subset,
                                               MARGIN    = 1,
                                               FUN       = weightFun,
                                               wts       = conversionList_subset$weights,
                                               na.cutoff = na.cutoff )

      # return
      return( higherRes_data_subset_weighted )



    }else{
      return( higherRes_data_subset )
    }
  }
  ###### Load data
  if ( verbose ) cat( 'Loading conversion values\n')
  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("HUC_conversion_weights", package = "h2hData", envir = environment())

  ###### convert to hucs
  dsMats <- vector( 'list' , length = 6)
  names(dsMats)[1:5] <- names( HUC_conversion_weights )
  names(dsMats)[6]   <- 'HUC12'

  for ( i in 1:( length( dsMats ) - 1 ) ){
    if ( verbose ) cat(
      paste0( 'Converting to ', names( dsMats )[ i ], '\n')
      )
    conversionWeights <- HUC_conversion_weights[[ i ]]
    mat.pre <- do.call(cbind,
                       pbapply::pblapply(X              = names( conversionWeights ),
                                         FUN            = applyWeights,
                                         conversionList = conversionWeights,
                                         higherRes_data = HUC12_matrix,
                                         na.cutoff      = na.cutoff))
    colnames( mat.pre ) <- names( conversionWeights )
    rownames( mat.pre ) <- rownames( HUC12_matrix )
    dsMats[[ i ]] <- mat.pre
  }
  dsMats$HUC12 <- HUC12_matrix

  rm(HUC_conversion_weights, envir = environment())

  return( dsMats )
}
#
# # convert to feather and name correctly
# length(dsMats)
# shead(dsMats$HUC2)
#
#
# library(feather)
# template <- feather::read_feather(path = "C:/Users/ssaxe/Documents/Projects/Model Evaluation2/FeatherFiles/Runoff_NLDAS-VIC_HUC2_month_1979-01-01_2017-06-01.feather")
#
# shead(template)
#
# ## names
# "Runoff_NLDAS-VIC_HUC2_month_1979-01-01_2017-06-01.feather"
#
#
# savePath <- "C:/Users/ssaxe/Documents/Projects/Model Evaluation2/FeatherFiles/"
# for (i in 1:length(dsMats)){
#   print(i)
#   saveName <- paste0(savePath,
#                      "Runoff_MWBM_",
#                      names(dsMats)[i],
#                      "_month_1949-01-01_2010-12-01.feather")
#   write_feather(x    = as.data.frame(dsMats[[i]]),
#                 path = saveName)
# }
#
#
#
#
#
#