wsim_composite.R 5.52 KB
Newer Older
1
#!/usr/bin/env Rscript
Dan Baston's avatar
Dan Baston committed
2 3 4 5 6 7 8 9 10 11 12 13 14 15

# Copyright (c) 2018 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

16 17
wsim.io::logging_init('wsim_composite')

18
suppressMessages({
Dan Baston's avatar
Dan Baston committed
19
  require(Rcpp)
20 21 22 23 24
})

'
Compute composite indicators

25
Usage: wsim_composite (--surplus=<file>)... (--deficit=<file>)... --both_threshold=<value> [--mask=<file>] [--clamp=<value>] --output=<file>
26 27

Options:
Dan Baston's avatar
Dan Baston committed
28 29
--surplus <file>...      One or more variables containing return periods that represent surpluses
--deficit <file>...      One or more variables containing return periods that represent deficits
30
--both_threshold <value> Threshold value for assigning a pixel to both surplus and deficit
31
--output <file>          Output file containing composite indicators
32
--mask <file>            Optional mask to use for computed indicators
33
--clamp <value>          Optional absolute value at which to clamp inputs
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
'->usage

clamp <- function(vals, minval, maxval) {
  pmax(pmin(vals, maxval), minval)
}

vals_for_depth_index <- function(arr, depth) {
  # Susbset the 3D array by providing a matrix of (row, col, level) triplets
  array(arr[cbind(as.vector(row(depth)),
                  as.vector(col(depth)),
                  as.vector(depth))],
        dim=dim(depth))
}

main <- function(raw_args) {
  args <- wsim.io::parse_args(usage,
                              raw_args,
51
                              types=list(both_threshold= 'integer',
52
                                         clamp= 'integer'))
53 54 55 56 57 58 59 60 61 62 63 64 65 66

  if (is.null(args$deficit)) {
    wsim.io::die_with_message("Must supply at least one deficit indicator.")
  }

  if (is.null(args$surplus)) {
    wsim.io::die_with_message("Must supply at least one surplus indicator.")
  }

  outfile <- args$output
  if (!wsim.io::can_write(outfile)) {
    wsim.io::die_with_message("Cannot open ", outfile, " for writing.")
  }

67
  surpluses <- wsim.io::read_vars_to_cube(args$surplus)
Dan Baston's avatar
Dan Baston committed
68
  wsim.io::info('Read surplus values:', paste(dimnames(surpluses)[[3]], collapse=", "))
69

70
  deficits <- wsim.io::read_vars_to_cube(args$deficit)
Dan Baston's avatar
Dan Baston committed
71
  wsim.io::info('Read deficit values:', paste(dimnames(deficits)[[3]], collapse=", "))
72

73 74 75 76 77 78 79
  if (is.null(args$mask)) {
    mask <- 1
  } else {
    mask_data <- wsim.io::read_vars(args$mask)$data[[1]]
    mask <- ifelse(!is.na(mask_data), 1, NA)
  }

Dan Baston's avatar
Dan Baston committed
80
  max_surplus_indices <- wsim.distributions::stack_which_max(surpluses)
81 82 83 84
  max_surplus_values <- vals_for_depth_index(surpluses, max_surplus_indices)
  if (!is.null(args$clamp)) {
    max_surplus_values <- clamp(max_surplus_values, -args$clamp, args$clamp)
  }
85

86
  wsim.io::info('Computed composite surplus.')
87

Dan Baston's avatar
Dan Baston committed
88
  min_deficit_indices <- wsim.distributions::stack_which_min(deficits)
89 90
  min_deficit_values <- vals_for_depth_index(deficits, min_deficit_indices)
  if (!is.null(args$clamp)) {
Dan Baston's avatar
Dan Baston committed
91
    min_deficit_values <- clamp(min_deficit_values, -args$clamp, args$clamp)
92
  }
93

94
  wsim.io::info('Computed composite deficit.')
95 96

  both_values <- ifelse(max_surplus_values > args$both_threshold & min_deficit_values < -(args$both_threshold),
97 98 99 100 101
                        # When above the threshold, take the largest absolute indicator
                        pmax(max_surplus_values, -min_deficit_values),
                        # When below the threshold, default to zero or NA, depending on the underlying
                        # indicators.
                        0 * max_surplus_values * min_deficit_values)
102 103

  cdf_data <- list(
104 105
    deficit= min_deficit_values*mask,
    deficit_cause= min_deficit_indices*mask,
106

107 108
    surplus= max_surplus_values*mask,
    surplus_cause= max_surplus_indices*mask,
109

110
    both= both_values*mask
111 112 113 114 115 116
  )

  cdf_attrs <- list(
    list(var="deficit", key="long_name", val="Composite Deficit Index"),

    list(var="deficit_cause", key="long_name", val="Cause of Deficit"),
117 118
    list(var="deficit_cause", key="flag_values", val=1:dim(deficits)[3], prec="byte"),
    list(var="deficit_cause", key="flag_meanings", val=paste(dimnames(deficits)[[3]], collapse=" "), prec="text"),
119 120 121 122

    list(var="surplus", key="long_name", val="Composite Surplus Index"),

    list(var="surplus_cause", key="long_name", val="Cause of Surplus"),
123 124
    list(var="surplus_cause", key="flag_values", val=1:dim(surpluses)[3], prec="byte"),
    list(var="surplus_cause", key="flag_meanings", val=paste(dimnames(surpluses)[[3]], collapse=" "), prec="text"),
125 126 127 128 129 130 131 132 133 134 135

    list(var="both", key="long_name", val="Composite Combined Surplus & Deficit Index"),
    list(var="both", key="threshold", val=args$both_threshold)
  )

  wsim.io::write_vars_to_cdf(cdf_data,
                             outfile,
                             extent= attr(deficits, 'extent'),
                             attrs= cdf_attrs,
                             prec=list(
                               deficit= "float",
136
                               deficit_cause= "byte",
137
                               surplus= "float",
138
                               surplus_cause= "byte",
139 140 141
                               both= "float"
                             ))

142
  wsim.io::info("Wrote composite indicators to", outfile)
143 144
}

145
tryCatch(main(commandArgs(trailingOnly=TRUE)), error=wsim.io::die_with_message)