@@ -32,33 +32,33 @@ actual_evaporation_waterbody_or_pervious <- function(
3232 fetch_usage <- create_accessor(usage_tuple )
3333 fetch_climate <- create_accessor(climate )
3434 fetch_soil <- create_accessor(soil_properties )
35-
35+
3636 rpot <- fetch_soil(" mean_potential_capillary_rise_rate" )
3737 epot_year <- fetch_climate(" epot_yr" )
38-
38+
3939 # Check input(s)
4040 stopifnot(! anyNA(rpot ))
41-
41+
4242 # Initialise result vector
4343 y <- numeric (length(epot_year ))
44-
44+
4545 # For water bodies, use the potential evaporation
4646 land_types <- fetch_usage(" land_type" )
4747 is_waterbody <- land_type_is_waterbody(land_types )
48-
48+
4949 y [is_waterbody ] <- epot_year [is_waterbody ]
50-
50+
5151 # if all block areas are waterbodies, return
5252 if (all(is_waterbody )) {
5353 return (y )
5454 }
55-
55+
5656 # indices of entries related to any other land_type
5757 i <- which(! is_waterbody )
58-
58+
5959 # otherwise calculate the real evapotranspiration
6060 stopifnot(all(epot_year [i ] > 0 )) # ???
61-
61+
6262 # determine the BAGROV parameter(s) for unsealed surfaces
6363 bagrov_values <- get_bagrov_parameter_unsealed(
6464 g02 = fetch_soil(" g02" )[i ],
@@ -69,41 +69,45 @@ actual_evaporation_waterbody_or_pervious <- function(
6969 epot_summer = fetch_climate(" epot_s" )[i ],
7070 mean_potential_capillary_rise_rate = rpot [i ]
7171 )
72-
72+
7373 if (! is.null(digits )) {
7474 bagrov_values <- cat_and_run(
7575 sprintf(" Rounding BAGROV parameters to %d digits" , digits ),
7676 round(bagrov_values , digits )
7777 )
7878 }
79-
79+
8080 available_water <-
8181 fetch_climate(" prec_yr" )[i ] +
8282 rpot [i ] +
8383 fetch_usage(" irrigation" )[i ]
84-
84+
8585 y [i ] <- real_evapo_transpiration(
8686 potential_evaporation = epot_year [i ],
8787 x_ratio = available_water / epot_year [i ],
8888 bagrov_parameter = bagrov_values
8989 # , use_abimo_algorithm = TRUE
9090 , ...
9191 )
92-
92+
9393 rises <- fetch_soil(" potential_capillary_rise" )
9494 depths <- fetch_soil(" depth_to_water_table" )
95-
95+
9696 # indices of entries related to non-water land_type and capillary rises < 0
9797 j <- which(! is_waterbody & rises < 0 )
98-
98+
9999 y [j ] <- y [j ] + (epot_year [j ] - y [j ]) * exp(depths [j ] / rises [j ])
100-
100+
101101 nas <- rep(NA_real_ , length(y ))
102-
102+
103103 structure(y , bagrovUnsealed = data.frame (
104104 bagrov_eff = `[<-`(nas , i , bagrov_values ),
105- factor_dry = `[<-`(nas , i , get_attribute(bagrov_values , " factor_dry" )),
106- factor_wet = `[<-`(nas , i , get_attribute(bagrov_values , " factor_wet" ))
105+ correction_irrigation_and_unknown_summer = `[<-`(nas , i , get_attribute(
106+ bagrov_values , " correction_irrigation_and_unknown_summer"
107+ )),
108+ correction_known_summer = `[<-`(nas , i , get_attribute(
109+ bagrov_values , " correction_known_summer"
110+ ))
107111 ))
108112}
109113
@@ -120,53 +124,53 @@ get_bagrov_parameter_unsealed <- function(
120124{
121125 # Initialise result vector
122126 y <- numeric (length = length(g02 ))
123-
127+
124128 is_forest <- land_type_is_forest(land_type )
125129 no_forest <- ! is_forest
126-
130+
127131 y [is_forest ] <- lookup_bagrov_forest(g02 [is_forest ])
128-
129- factor_dry <- ifelse(
130- test = irrigation > 0 & is_dry_summer(prec_summer , epot_summer ),
131- yes = irrigation_in_dry_summer_correction_factor(irrigation [no_forest ]),
132+
133+ # It seems that in the original Abimo code, values of zero in the "summer"
134+ # columns were used to indicate "missing"
135+ correction_irrigation_and_unknown_summer <- ifelse(
136+ test = irrigation > 0 & (prec_summer < = 0 & epot_summer < = 0 ),
137+ yes = irrigation_and_unknown_summer_correction(irrigation [no_forest ]),
132138 no = 1
133139 )
134-
140+
135141 y [no_forest ] <- lookup_bagrov_unsealed(g02 [no_forest ], veg_class [no_forest ]) *
136- factor_dry [no_forest ]
137-
138- # in case of a "wet" summer, correct the BAGROV parameter with a factor
139- factor_wet <- ifelse(
140- test = is_wet_summer(prec_summer , epot_summer ),
141- yes = wet_summer_correction_factor(
142- water_availability =
143- prec_summer +
144- irrigation +
145- mean_potential_capillary_rise_rate ,
142+ correction_irrigation_and_unknown_summer [no_forest ]
143+
144+ # in case of known "summer" values for precipitation and evaporation, correct
145+ # the BAGROV parameter with a factor
146+ correction_known_summer <- ifelse(
147+ test = prec_summer > 0 & epot_summer > 0 ,
148+ yes = summer_correction(
149+ water_availability = prec_summer + irrigation + mean_potential_capillary_rise_rate ,
146150 epot_summer = epot_summer
147151 ),
148152 no = 1
149153 )
150-
154+
151155 structure(
152- y * factor_wet ,
153- factor_dry = factor_dry ,
154- factor_wet = factor_wet
156+ y * correction_known_summer ,
157+ correction_irrigation_and_unknown_summer = correction_irrigation_and_unknown_summer ,
158+ correction_known_summer = correction_known_summer
155159 )
156160}
157161
158162# lookup_bagrov_forest ---------------------------------------------------------
159163lookup_bagrov_forest <- function (g02 )
160164{
161165 n <- length(g02 )
162-
166+
163167 if (n == 0L ) {
164168 return (numeric (0 ))
165169 }
166-
170+
167171 breaks <- c(- Inf , 10.0 , 25.0 , Inf )
168172 values <- c(3.0 , 4.0 , 8.0 )
169-
173+
170174 index <- if (n > 1L ) {
171175 findInterval(g02 , breaks , left.open = TRUE )
172176 } else if (g02 < = breaks [2L ]) {
@@ -176,7 +180,7 @@ lookup_bagrov_forest <- function(g02)
176180 } else {
177181 3L
178182 }
179-
183+
180184 values [index ]
181185}
182186
@@ -185,43 +189,43 @@ lookup_bagrov_unsealed <- function(g02, veg_class, do_correction = TRUE)
185189{
186190 # Calculate the k index (integer)
187191 k <- veg_class_to_k_index(veg_class )
188-
192+
189193 # Calculate result based on the k index
190194 y <-
191195 BAGROV_COEFFICIENTS [k ] +
192196 BAGROV_COEFFICIENTS [k + 1L ] * g02 +
193197 BAGROV_COEFFICIENTS [k + 2L ] * g02 ^ 2
194-
198+
195199 # Return y if no correction is required
196200 if (! do_correction ) {
197201 return (y )
198202 }
199-
203+
200204 # Apply correction where needed
201205 i <- which(
202206 (y > = 2.0 & veg_class < 60 ) |
203207 (g02 > = 20.0 & veg_class > = 60 )
204208 )
205-
209+
206210 y [i ] <-
207211 BAGROV_COEFFICIENTS [k [i ] - 2L ] * g02 [i ] +
208212 BAGROV_COEFFICIENTS [k [i ] - 1L ]
209-
213+
210214 y
211215}
212216
213217# veg_class_to_k_index -------------------------------------------------------------
214218veg_class_to_k_index <- function (veg_class )
215219{
216220 k <- as.integer(ifelse(veg_class < 50 , veg_class / 5 , veg_class / 10 + 5 ))
217-
221+
218222 # make sure that k is at least 1
219223 k <- pmax(1L , k )
220-
224+
221225 # if k is at least 4, reduce it by one
222226 selected <- k > = 4L
223227 k [selected ] <- k [selected ] - 1L
224-
228+
225229 5L * pmin(k , 13L ) - 2L
226230}
227231
@@ -245,47 +249,31 @@ BAGROV_COEFFICIENTS <- c(
245249 0.33895 , 3.721 , 6.69999 , - 0.07 , 0.013
246250)
247251
248- # is_dry_summer ----------------------------------------------------------------
249- # TODO: Remove redundancy with is_wet_summer.
250- # Variables are (almost!) one another's opposite!
251- is_dry_summer <- function (prec_summer , epot_summer )
252- {
253- prec_summer < = 0 & epot_summer < = 0
254- }
255-
256- # irrigation_in_dry_summer_correction_factor -----------------------------------
257- irrigation_in_dry_summer_correction_factor <- function (irrigation )
252+ # irrigation_and_unknown_summer_correction -------------------------------------
253+ irrigation_and_unknown_summer_correction <- function (irrigation )
258254{
259255 0.9985 + 0.00284 * irrigation - 0.00000379762 * irrigation ^ 2
260256}
261257
262- # is_wet_summer ----------------------------------------------------------------
263- # TODO: Remove redundancy with is_dry_summer.
264- # Variables are (almost!) one another's opposite!
265- is_wet_summer <- function (prec_summer , epot_summer )
266- {
267- prec_summer > 0 & epot_summer > 0
268- }
269-
270- # wet_summer_correction_factor -------------------------------------------------
258+ # summer_correction ------------------------------------------------------------
271259# ' @importFrom stats approx
272- wet_summer_correction_factor <- function (
260+ summer_correction <- function (
273261 water_availability , epot_summer , use_abimo_approx = TRUE
274262)
275263{
276264 xout <- water_availability / epot_summer
277- x <- WET_SUMMER_CORRECTION_MATRIX [, " water_availability" ]
278- y <- WET_SUMMER_CORRECTION_MATRIX [, " correction_factor" ]
279-
265+ x <- SUMMER_CORRECTION_MATRIX [, " water_availability" ]
266+ y <- SUMMER_CORRECTION_MATRIX [, " correction_factor" ]
267+
280268 if (use_abimo_approx ) {
281269 interpolate(x = x , y = y , xout = xout )
282270 } else {
283271 select_columns(stats :: approx(x = x , y = y , xout = xout , rule = 2L ), " y" )
284272 }
285273}
286274
287- # WET_SUMMER_CORRECTION_MATRIX -------------------------------------------------
288- WET_SUMMER_CORRECTION_MATRIX <- matrix (
275+ # SUMMER_CORRECTION_MATRIX ---- -------------------------------------------------
276+ SUMMER_CORRECTION_MATRIX <- matrix (
289277 ncol = 2L ,
290278 byrow = TRUE ,
291279 dimnames = list (
0 commit comments