3434# ' Note that the default parameters are tailored for Australia and may not be suitable for other
3535# ' regions.
3636# '
37- # ' @inheritParams ref_density
37+ # ' @inheritParams ref_density
3838# ' @param samples A matrix or data.frame containing x, y, predicted-RS, and observed-RS values
3939# ' (in that specific order) for benchmark samples (also known as reference sites). If the
4040# ' \code{data} argument is a SpatRaster, you can provide only the x and y coordinates of the
4141# ' benchmark samples. In this case, the corresponding values will be extracted from the raster
4242# ' layers. Consider extra time for sample value extraction in this case.
43- # ' @param ref_density A matrix or \strong{reference_density} object of normalised HCAS
44- # ' reference density (see \code{\link{ref_density}} and \code{\link{normalise}}).
45- # ' @param xy_stats A vector, mean and standard deviation of coordinates for centre and
43+ # ' @param ref_density A matrix or \strong{reference_density} object of normalised HCAS
44+ # ' reference density (see \code{\link{ref_density}} and \code{\link{normalise}}).
45+ # ' @param xy_stats A vector, mean and standard deviation of coordinates for centre and
4646# ' scaling the coordinate to use as a penalty. The order should be: mean(x), mean(y), sd(x), sd(y).
4747# ' This argument helps achieving consistent results when running benchmarking over multiple tiles.
4848# ' @param xy_penalty Numeric. The spatial distance penalty value for selecting benchmark points.
5151# ' See details section for more information on distance calculation.
5252# ' @param k_pred Integer. Number of nearest predicted RS samples to take.
5353# ' @param k_obs Integer. Number of nearest observed RS sample to takes.
54- # ' @param bin_width Numeric. Specifies the bin width of the reference density. If \code{ref_density} is
55- # ' a \strong{reference_density} object, this value can be read from its attributes and may be left \code{NULL}.
56- # ' The bin width must be consistent between the reference density creation and the benchmarking step to
57- # ' ensure condition is accurately calculated.
58- # ' @param interpolate Logical. Whether to interpolate the reference density for a smoother result.
59- # ' @param offset Integer. Specifies the number of reference density bins that were ignored during
60- # ' normalisation (see \code{\link{normalise}}). If \code{ref_density} is a \strong{reference_density} object,
61- # ' this value can be read from its attributes and may be set \code{NULL}. Similar to bin-width,
62- # ' the \code{offset} must be consistent between the reference density normalisation and the benchmarking
63- # ' step to ensure condition is accurately calculated.
54+ # ' @param bin_width Numeric. Specifies the bin width of the reference density. If \code{ref_density} is
55+ # ' a \strong{reference_density} object, this value can be read from its attributes and may be left \code{NULL}.
56+ # ' The bin width must be consistent between the reference density creation and the benchmarking step to
57+ # ' ensure condition is accurately calculated.
58+ # ' @param interpolate Logical. Whether to interpolate the reference density for a smoother result.
59+ # ' @param offset Integer. Specifies the number of reference density bins that were ignored during
60+ # ' normalisation (see \code{\link{normalise}}). If \code{ref_density} is a \strong{reference_density} object,
61+ # ' this value can be read from its attributes and may be set \code{NULL}. Similar to bin-width,
62+ # ' the \code{offset} must be consistent between the reference density normalisation and the benchmarking
63+ # ' step to ensure condition is accurately calculated.
6464# ' @param confidence Numeric. The confidence value for LDC methods. See details below..
6565# ' @param lambda Numeric. The lambda param for LDC Cauchy weighting...
6666# ' @param exclude_slef Logical. To exclude a benchmark point from assessing itself.
67- # ' @param drop_features Integer vector. Completely remove RS variables from the benchmarking process.
68- # ' Positions are 1-based within the RS feature set, not the full input column order. For
69- # ' consistency, it is recommended to exclude the same variables used in the reference density step; unless
70- # ' you have a specific reason not to.
67+ # ' @param drop_features Integer vector. Completely remove RS variables from the benchmarking process.
68+ # ' Positions are 1-based within the RS feature set, not the full input column order. For
69+ # ' consistency, it is recommended to exclude the same variables used in the reference density step; unless
70+ # ' you have a specific reason not to.
7171# ' @param make_su Logical. To make the uncertainty map or not.
7272# ' @param ... Additional arguments for writing raster outputs e.g. \code{filename},
7373# ' \code{overwrite}, and \code{wopt} from terra \code{\link[terra]{predict}}.
7474# '
75- # ' @seealso \code{\link{ref_density}}, \code{\link{normalise}}, and \code{\link{calibrate}}
75+ # ' @seealso \code{\link{ref_density}}, \code{\link{normalise}}, and \code{\link{calibrate}}
7676# '
7777# ' @return A matrix or SpatRaster, depending on the inputs.
7878# ' @export
8484# '
8585# '
8686# ' }
87- benchmark <- function (
88- data ,
89- samples ,
90- ref_density ,
87+ benchmark <- function (
88+ data ,
89+ samples ,
90+ ref_density ,
9191 xy_stats = c(0 , 0 , 1 , 1 ),
9292 xy_penalty = 0.0 ,
9393 radius_km = 200 ,
@@ -104,64 +104,64 @@ benchmark <- function(
104104 num_threads = - 1 ,
105105 ... ) {
106106
107- # check k_pred and k_obs
108- if (k_pred < k_obs ) stop(" 'k_obs' must be smaller or equal to 'k_pred'." )
109- # check samples and reference density
110- samples <- if (.is_mat(samples )) .check_mat(samples ) else stop(" 'samples' must be a matrix or convertible to one." )
111- ref_density <- if (.is_mat(ref_density )) .check_mat(ref_density ) else stop(" 'ref_density' must be a matrix or convertible to one." )
112- if (nrow(ref_density ) != ncol(ref_density )) warning(" Reference density dimensions are not equal!\n " )
113-
114- if (methods :: is(ref_density , " reference_density" )) {
115- # check for reference density bin_width consistency
116- if (is.null(bin_width )) {
117- bin_width <- attributes(ref_density )$ bin.width
118- } else {
119- if (bin_width != attributes(ref_density )$ bin.width ) {
120- warning(" Provided 'bin_width' differs from reference density attribute." )
121- }
122- }
123- # check for reference density offset consistency
124- if (is.null(offset )) {
125- offset <- attributes(ref_density )$ offset
126- } else {
127- if (offset != attributes(ref_density )$ offset ) {
128- warning(" Provided 'offset' differs from reference density attribute." )
129- }
130- }
131- }
132-
133- # interpolate reference density
134- if (interpolate ) {
135- ref_density <- terra :: as.matrix(
136- terra :: disagg(terra :: rast(ref_density ), fact = 2 , method = " bilinear" ),
137- wide = TRUE
138- )
139- # update the binwidth and offset
107+ # check k_pred and k_obs
108+ if (k_pred < k_obs ) stop(" 'k_obs' must be smaller or equal to 'k_pred'." )
109+ # check samples and reference density
110+ samples <- if (.is_mat(samples )) .check_mat(samples ) else stop(" 'samples' must be a matrix or convertible to one." )
111+ ref_density <- if (.is_mat(ref_density )) .check_mat(ref_density ) else stop(" 'ref_density' must be a matrix or convertible to one." )
112+ if (nrow(ref_density ) != ncol(ref_density )) warning(" Reference density dimensions are not equal!\n " )
113+
114+ if (methods :: is(ref_density , " reference_density" )) {
115+ # check for reference density bin_width consistency
116+ if (is.null(bin_width )) {
117+ bin_width <- attributes(ref_density )$ bin.width
118+ } else {
119+ if (bin_width != attributes(ref_density )$ bin.width ) {
120+ warning(" Provided 'bin_width' differs from reference density attribute." )
121+ }
122+ }
123+ # check for reference density offset consistency
124+ if (is.null(offset )) {
125+ offset <- attributes(ref_density )$ offset
126+ } else {
127+ if (offset != attributes(ref_density )$ offset ) {
128+ warning(" Provided 'offset' differs from reference density attribute." )
129+ }
130+ }
131+ }
132+
133+ # interpolate reference density
134+ if (interpolate ) {
135+ ref_density <- terra :: as.matrix(
136+ terra :: disagg(terra :: rast(ref_density ), fact = 2 , method = " bilinear" ),
137+ wide = TRUE
138+ )
139+ # update the binwidth and offset
140140 bin_width <- bin_width / 2
141141 offset <- offset * 2
142- }
143- # get the bin number after interpolation
144- bin_num <- min(dim(ref_density ))
142+ }
143+ # get the bin number after interpolation
144+ bin_num <- min(dim(ref_density ))
145+
146+ if (.is_mat(data )) {
147+ # check and convert to matrix
148+ data <- .check_mat(data )
145149
146- if (.is_mat(data )) {
147- # check and convert to matrix
148- data <- .check_mat(data )
149-
150- if (ncol(samples ) != ncol(data )) {
151- stop(" Samples must include all raster values (matching column count with 'data')." )
152- }
153-
154- keep_features <- .keep_rs_features(drop_features , .num_rs_vars_mat(samples , " samples" ))
155- samples <- .subset_rs_mat(samples , keep_features )
156- data <- .subset_rs_mat(data , keep_features )
157-
158- tryCatch(
159- {
160- output <- benchmarking(
161- model = list (),
150+ if (ncol(samples ) != ncol(data )) {
151+ stop(" Samples must include all raster values (matching column count with 'data')." )
152+ }
153+
154+ keep_features <- .keep_rs_features(drop_features , .num_rs_vars_mat(samples , " samples" ))
155+ samples <- .subset_rs_mat(samples , keep_features )
156+ data <- .subset_rs_mat(data , keep_features )
157+
158+ tryCatch(
159+ {
160+ output <- .benchmarking(
161+ model = list (),
162162 newdata = data , # rast_stack arg
163163 sample_vals = samples ,
164- ref_density = ref_density ,
164+ ref_density = ref_density ,
165165 xy_stats = xy_stats ,
166166 xy_penalty = xy_penalty ,
167167 radius_km = radius_km ,
@@ -171,15 +171,15 @@ benchmark <- function(
171171 offset = offset ,
172172 k_env = k_pred ,
173173 k_rs = k_obs ,
174- confidence = confidence ,
175- lambda = lambda ,
176- exclude_slef = exclude_slef ,
177- make_su = make_su ,
178- num_threads = num_threads
179- )
180- },
174+ confidence = confidence ,
175+ lambda = lambda ,
176+ exclude_slef = exclude_slef ,
177+ make_su = make_su ,
178+ num_threads = num_threads
179+ )
180+ },
181181 error = function (cond ) {
182- stop(" HCAS benchmarking failed!\n " , cond )
182+ stop(" HCAS benchmarking C++ function failed!\n " , cond )
183183 }
184184 )
185185 } else if (.is_rast(data )) {
@@ -188,26 +188,26 @@ benchmark <- function(
188188 # check and convert to SpatRaster object
189189 data <- .check_rast(data )
190190
191- # sample extraction if needed
192- if (ncol(samples ) == 2 ) {
193- cat(" Extracting sample values...\n " )
194- samples <- cbind(samples , as.matrix(terra :: extract(data , samples , ID = FALSE )))
195- } else if ((ncol(samples ) - 2 ) != terra :: nlyr(data )) {
196- stop(" Sample feature count does not match number of raster layers." )
197- }
198-
199- keep_features <- .keep_rs_features(drop_features , terra :: nlyr(data ) / 2L )
200- samples <- .subset_rs_mat(samples , keep_features )
201- data <- .subset_rs_rast(data , keep_features )
202-
203- tryCatch(
204- {
205- output <- terra :: interpolate(
206- object = data ,
191+ # sample extraction if needed
192+ if (ncol(samples ) == 2 ) {
193+ cat(" Extracting sample values...\n " )
194+ samples <- cbind(samples , as.matrix(terra :: extract(data , samples , ID = FALSE )))
195+ } else if ((ncol(samples ) - 2 ) != terra :: nlyr(data )) {
196+ stop(" Sample feature count does not match number of raster layers." )
197+ }
198+
199+ keep_features <- .keep_rs_features(drop_features , terra :: nlyr(data ) / 2L )
200+ samples <- .subset_rs_mat(samples , keep_features )
201+ data <- .subset_rs_rast(data , keep_features )
202+
203+ tryCatch(
204+ {
205+ output <- terra :: interpolate(
206+ object = data ,
207207 model = list (),
208- fun = benchmarking ,
208+ fun = . benchmarking ,
209209 sample_vals = samples ,
210- ref_density = ref_density ,
210+ ref_density = ref_density ,
211211 xy_stats = xy_stats ,
212212 xy_penalty = xy_penalty ,
213213 radius_km = radius_km ,
@@ -217,16 +217,16 @@ benchmark <- function(
217217 offset = offset ,
218218 k_env = k_pred ,
219219 k_rs = k_obs ,
220- confidence = confidence ,
221- lambda = lambda ,
222- exclude_slef = exclude_slef ,
223- make_su = make_su ,
224- num_threads = num_threads ,
225- ...
226- )
220+ confidence = confidence ,
221+ lambda = lambda ,
222+ exclude_slef = exclude_slef ,
223+ make_su = make_su ,
224+ num_threads = num_threads ,
225+ ...
226+ )
227227 },
228228 error = function (cond ) {
229- stop(" HCAS benchmarking failed!\n " , cond )
229+ stop(" HCAS benchmarking C++ function failed!\n " , cond )
230230 }
231231 )
232232
@@ -240,15 +240,15 @@ benchmark <- function(
240240}
241241
242242
243- # a function to handling predicting with terra
244- benchmarking <- function (model , newdata , make_su , ... ){
245- nr <- nrow(newdata )
246- nc <- make_su + 1
247- col_names <- c(" condition" , " su" )[1 : nc ]
248-
249- dat <- as.matrix(newdata )
250-
251- tryCatch(
243+ # a function to handling predicting with terra
244+ . benchmarking <- function (model , newdata , make_su , ... ){
245+ nr <- nrow(newdata )
246+ nc <- make_su + 1
247+ col_names <- c(" condition" , " su" )[1 : nc ]
248+
249+ dat <- as.matrix(newdata )
250+
251+ tryCatch(
252252 {
253253 hcas_cond <- bench_cpp(
254254 raster_vals = dat ,
0 commit comments