Areal Reduction Factors see https://tonyladson.wordpress.com/2016/02/27/areal-reduction-factors/

ARF.R

##################################################################################
# 
# ARFs in the new ARR
#
# tony.ladson@gmail.com
# 25 Feb 2016
#
##################################################################################

library(repmis)
library(stringr)
library(testthat)




# Download parameters
my.url <- 'https://dl.dropboxusercontent.com/u/10963448/ARF_params.csv'
params <- source_data(my.url, header = TRUE)

# set parameters up in columns rather than rows
x <- data.frame( t(params[ ,-1 ]))
names(x) <- params[ ,1]
params <- x

params



ARF_long <- function(area, duration, aep, region, params) {
  # Area in km-squre
  # Duration in minutes
  # aep as a fracion e.g. 0.01 
  # region must be one of the 10 valid regions
  # params is a data frame of parameters a,...,i for all regions
  
  # Check the region is valid
  
  all.regions <- names(params)
  all.regions.txt <- str_c(all.regions, collapse = ', ')

  if(!(region %in% all.regions)) stop (str_c('Invalid region. You input "', region, '". Valid regions are ', all.regions.txt))
  
  # For the select region, assignment the parameters to a,...,i
  
  for(i in 1:9) {
    assign(letters[i], params[ ,region][i])
  }
  
  min(1, (1 - a*(area^b - c*log10(duration)) * duration ^-d +
            e*area^f*duration^g * (0.3 + log10(aep)) + 
            h*10^(i*area*duration/1440) * (0.3 + log10(aep))))
    
  
}


ARF_short <- function(area, duration, aep) {

  a <- 0.287	
  b <- 0.265
  c <- 0.439
  d <- 0.36
  e <- 0.00226	
  f <- 0.226
  g <- 0.125
  h <- 0.0141
  i <- -0.021
  j <- 0.213
  
  min(1, (1 - a*(area^b - c*log10(duration)) * duration^(-d) +
          e*area^f*duration^g * (0.3 + log10(aep)) + 
          h * area^j * 10^(i* (1/1440) * (duration - 180)^2)  * (0.3 + log10(aep))))
}



ARF <- function(area, duration, aep, region = NULL, params = NULL) {
  
  # We only need a region and parameters if duration is greater than 12 hours so
  # these arguements are optional
  
  # Define the functions we may need
  
  # Checking inputs
  
  if(duration < 30) stop('Duration must be greater than 30 mins')
  if(duration > 7*24*60) stop('Duration must be less than 7 days')
  if(area > 30000) stop('Area must be less than 30,000 km-square')
  
  
  
  if(duration < 1440 & area >= 1000) 
    stop(str_c("Generalized equations are not applicable for short durations on large areas. ",
               "If area > 1000, duration must be greater than 24 hours (1440 mins)"))
      
     
  # Following the procedure in Scott Podger's spreadsheet
  
  
  if(duration >= 1440){
    if(area >= 10){
      return(max(ARF_long(area, duration, aep, region, params), ARF_short(area, duration = 720, aep)))
    }
    
    # area < 10
    # interpolate based on area between the long duration ARF for 10 km-square
    # and an ARF of 1
    
    ARF.long.10 <- ARF_long(10, duration, aep, region, params)
    return(ARF.long.10 - (ARF.long.10 - 1)*(10 - area)/10)
    
  }
  if(duration <= 720){
    if(area >= 10){
      return(ARF_short(area, duration, aep))
    }
   
    # area < 10
    # interpolate based on area between the short duration ARF for 10 km-square
    # and an ARF of 1
    
    ARF.short.10 <- ARF_short(10, duration, aep)
    return(ARF.short.10 - (ARF.short.10 - 1)*(10 - area)/10)
  }
  
  # duration between 720 and 1440 i.e. between 12 hours and 24 hours
  
  # 1. Calculate short duration for 720 min (12 hours)
  
  ARF.short.720 <- ARF_short(area, 720, aep)
  
  # 2. Calculate long duration for 1440 min (24 hours)
  
  ARF.long.1440 <- ARF_long(area, 1440, aep, region, params)
  
  # 3.  Use the maximum of 1 and 2
  
  ARF.max <- max(ARF.short.720, ARF.long.1440)
  
  # 4. interpolate between 3 and 1
  
  return(ARF.max - (ARF.max - ARF.short.720)*(1440 - duration)/(1440 - 720))
  
}



# Input checks

# Expect error for invalid region

expect_error(ARF(area = 20, duration = 1500, aep = 0.01, region = 'East Coastal', params = params))

# Expect no error for valid regions

valid.regions <- c('East Coast North', 'Semi-arid Inland QLD', 
                   'Tasmania', 'SW WA', 'Central NSW', 'SE Coast',
                   'Southern Semi-arid', 'Southern Temperate', 'Northern Coastal',
                   'Inland Arid')

# ARF should equal 1 if area is zero (and region is valid)

for(region in valid.regions) {
  expect_equal(ARF(area = 0, duration = 1500, aep = 0.01, region = region, params = params), 1)
}




# Expect error for area > 30000

expect_error(ARF(area = 30001, duration = 1500, aep = 0.01, region = 'Northern Coastal', params = params))

# Expect error for duration > 7 days

expect_error(ARF(area = 30000, duration = 7*24*60+1, aep = 0.01, region = 'Northern Coastal', params = params))

# Expect error if duration < 1440 and area > 1000

expect_error(ARF(area = 1001, duration = 1439, aep = 0.01, region = 'Northern Coastal', params = params))

# Expect error for duration < 30 mins

expect_error(ARF(area = 1001, duration = 29, aep = 0.01, region = 'Northern Coastal', params = params))

# Small areas
# area = 0 should result in ARF = 1

expect_equal(ARF(area = 0, duration = 719, aep = 0.01, region = 'Northern Coastal', params = params), 1)


expect_equal(ARF(area = 1, duration = 1500, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.996993,
             tolerance = 0.0001,
             scale = 0.996993
)

expect_equal(ARF(area = 0.1, duration = 1500, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.999699,
             tolerance = 0.0001,
             scale = 0.999699
)


#
ARF(area = 1, duration = 0.001, aep = 0.01, region = 'Northern Coastal', params = params)



# Test cases
# Six test cases for numbers
# areas <10 and >= 10
# durations < 720, 720 to 1440 and > 1440
# i.e. < 12 hour, between 12 and 24 hours and greater than 24 hours

# Long duration, duration >= 1440, area > 10




expect_equal(ARF(area = 3000, duration = 1500, aep = 0.01, region = 'SE Coast', params = params), 
             expected = 0.86839335, 
             tolerance = 0.0001,
             scale = 0.86839335
)


expect_equal(ARF(area = 20, duration = 1500, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.961831, 
             tolerance = 0.0001,
             scale = 0.961831
)

# Long duration, duration >= 1440, area < 10

expect_equal(ARF(area = 9, duration = 1500, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.972934, 
             tolerance = 0.0001,
             scale = 0.972934
)

expect_equal(ARF(area = 8, duration = 4000, aep = 0.01, region = 'SE Coast', params = params), 
             expected = 0.991555804, 
             tolerance = 0.0001,
             scale = 0.991555804
)


# Short duration large area

expect_equal(ARF(area = 100, duration = 700, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.917252, 
             tolerance = 0.0001,
             scale = 0.917252
)

expect_equal(ARF(area = 1000, duration = 500, aep = 0.01, region = 'SE Coast', params = params), 
             expected = 0.802047829, 
             tolerance = 0.0001,
             scale = 0.802047829
)




# Short duration small area

expect_equal(ARF(area = 9, duration = 700, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.972344, 
             tolerance = 0.0001,
             scale = 0.972344
)

expect_equal(ARF(area = 6, duration = 500, aep = 0.01, region = 'SE Coast', params = params), 
             expected = 0.978753394, 
             tolerance = 0.0001,
             scale = 0.978753394
)



# middle duration (12 hour - 24 hour)  large area

expect_equal(ARF(area = 1000, duration = 1000, aep = 0.01, region = 'SE Coast', params = params), 
             expected = 0.859737827, 
             tolerance = 0.0001,
             scale = 0.859737827
)


expect_equal(ARF(area = 2000, duration = 1500, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.8640315, 
             tolerance = 0.0001,
             scale = 0.8640315
)


# middle duration (12 hour - 24 hour)  small area 

expect_equal(ARF(area = 9, duration = 1000, aep = 0.01, region = 'Northern Coastal', params = params), 
             expected = 0.971236, 
             tolerance = 0.0001,
             scale = 0.971236
)


expect_equal(ARF(area = 7, duration = 1000, aep = 0.01, region = 'SE Coast', params = params), 
             expected = 0.97991314, 
             tolerance = 0.0001,
             scale = 0.97991314
)







# for durations less than 12 hours we don't need region
# so results should be the same with or without a region

expect_equal(ARF(9, 100, 0.01, region = 'Northern Coastal', params = params), 
             expected = ARF(9, 100, 0.01), tolerance = 0.0001,
             scale =ARF(9, 100, 0.01))


########

my.aep <- c(0.5, 0.2, 0.1, 0.05,0.02, 0.01)
my.duration <- 60 * c(0.5, 1, 2, 3, 6, 12, 24, 48, 72, 96, 120, 144, 168)


x <- expand.grid(aep = my.aep, duration = my.duration)

x$res <- mapply(ARF, aep = x$aep, duration = x$duration, MoreArgs = list(area = 519, region = 'SE Coast', params = params))
 

library(reshape2)

out <- dcast(x, duration ~ aep)  # convert to row x column
out <- out[  , c(1, 7:2)] # sort column by AEP

# write to clipboard
clip <- pipe("pbcopy", "w")                       
write.table(out, file=clip, sep = '\t', row.names = FALSE)                               
close(clip)