forked from ropensci/rangr
-
Notifications
You must be signed in to change notification settings - Fork 0
/
initialise.R
798 lines (662 loc) · 25 KB
/
initialise.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
#' Prepare Data Required To Perform A Simulation
#'
#' This funtion generates a `sim_data` object containing all the necessary
#' information required to run a simulation by the [`sim`] function. The
#' input maps (`n1_map` and `K_map`) can be in the Cartesian or longitude/latitude
#' coordinate system.
#'
#'
#' The most time-consuming part of computations performed by the [`sim`]
#' function is the simulation of dispersal. To speed it up, a list containing
#' indexes of target cells at a specified distance from a focal cell
#' is calculated in advance and stored in a `dlist` slot.
#' The `max_dist` parameter sets
#' the maximum distance at which this pre-calculation is performed. If `max_dist`
#' is `NULL`, it is set to 0.99 quantile from the `kernel_fun`.
#' All distance calculations are always based on metres if the input maps are
#' latitude/longitude. For planar input maps, distances are calculated in map
#' units, which are typically metres, but check the [`crs()`][terra::crs]
#' if in doubt.
#'
#' If the input maps are in the Cartesian coordinate system and the grid cells
#' are squares, then
#' the distances between cells are calculated using the [`distance`][terra::distance]
#' function from the `terra` package. These distances are later divided by the
#' resolution of the input maps.
#'
#' For input maps with grid cells in shapes other than squares (e.g. with
#' rectangular cells or longitude/latitude coordinate system), the distance
#' resolution is calculated by finding the shortest distance between each
#' "queen" type neighbor. All distances calculated by the [`distance`][terra::distance]
#' function are further divided by this distance resolution.
#' To avoid discontinuities in the distances at which the target cells are located,
#' an additional parameter `dist_bin` is calculated as half of the maximum
#' distance between each "queen" type neighbour. It is used to expand the
#' distances at which target cells are located from a single number to a range.
#'
#' NA in the input maps represents cells outside the study area.
#'
#' The [`K_get_interpolation`] function can be used to prepare `K_map` that changes
#' over time. This may be useful, when simulating environmental change or
#' exploring the effects of ecological disturbances.
#'
#'
#'
#' @param n1_map [`SpatRaster`][terra::SpatRaster-class] object with one layer;
#' population numbers in every grid cell at the first time step
#' @param K_map [`SpatRaster`][terra::SpatRaster-class] object with one layer;
#' carrying capacity map (if K is constant across time) or maps (if K is
#' time-varying)
#' @param K_sd numeric vector of length 1 with value `>= 0` (default 0);
#' this parameter can be used if additional environmental stochasticity
#' is required; if `K_sd > 0`, a random numbers are generated from a log-normal
#' distribution with the mean `K_map` and standard deviation `K_sd`
#' @param r numeric vector of length 1; intrinsic population growth rate
#' @param r_sd numeric vector of length 1 with value `>= 0` (default `0`);
#' if additional demographic stochasticity is required, `r_sd > 0` is
#' the standard deviation for a normal distribution around `r`
#' (defined for each time step)
#' @param growth character vector of length 1; the name of a population growth
#' function, either defined in [`growth`] or provided by
#' the user (case-sensitive, default [`"gompertz"`][growth])
#' @param A numeric vector of length 1; strength of the Allee effect
#' (see the [`growth`] function)
#' @param dens_dep character vector of length 1 specifying if the probability
#' of settling in a target grid cell is (case-sensitive, default `"K2N"`):
#' \itemize{
#' \item "none" - fully random,
#' \item "K" - proportional to the carrying capacity of a target cell,
#' \item "K2N" - density-dependent, i.e. proportional to the ratio of
#' carrying capacity of a target cell to the number of individuals
#' already present in a target cell
#' }
#' @param border character vector of length 1 defining how to deal
#' with borders (case-sensitive, default `"absorbing"`):
#' \itemize{
#' \item "reprising" - cells outside the study area are not allowed
#' as targets for dispersal
#' \item "absorbing" - individuals that disperse outside the study area
#' are removed from the population
#' }
#' @param kernel_fun character vector of length 1; name of a random number
#' generation function defining a dispersal kernel (case-sensitive,
#' default `"rexp"`)
#' @param ... any parameters required by `kernel_fun`
#' @param max_dist numeric vector of length 1; maximum distance of dispersal
#' to pre-calculate target cells
#' @param calculate_dist logical vector of length 1; determines if target cells
#' will be precalculated
#' @param dlist list; target cells at a specified distance calculated
#' for every cell within the study area
#' @param progress_bar logical vector of length 1; determines if progress bar
#' for calculating distances should be displayed
#' @param quiet logical vector of length 1; determines if messages should be displayed
#'
#' @return Object of class `sim_data` which inherits from `list`. This object
#' contains all necessary information to perform a simulation using
#' [`sim`] function.
#'
#' @export
#'
#'
#' @examples
#' \dontrun{
#'
#' # input maps
#' library(terra)
#'
#' n1_small <- rast(system.file("input_maps/n1_small.tif", package = "rangr"))
#' K_small <- rast(system.file("input_maps/K_small.tif", package = "rangr"))
#' K_small_changing <- rast(system.file("input_maps/K_small_changing.tif",
#' package = "rangr"))
#' n1_small_lon_lat <- rast(system.file("input_maps/n1_small_lon_lat.tif", package = "rangr"))
#' K_small_lon_lat <- rast(system.file("input_maps/K_small_lon_lat.tif", package = "rangr"))
#'
#' # basic example
#' sim_data_1 <- initialise(
#' n1_map = n1_small,
#' K_map = K_small,
#' r = log(2),
#' rate = 1 / 1e3
#' )
#'
#' # example with changing environment
#' K_interpolated <- K_get_interpolation(
#' K_small_changing,
#' K_time_points = c(1, 25, 50)
#' )
#'
#' sim_data_2 <- initialise(
#' n1_map = n1_small,
#' K_map = K_interpolated,
#' r = log(2),
#' rate = 1 / 1e3
#' )
#'
#' # example with lon/lat rasters
#' sim_data_3 <- initialise(
#' n1_map = n1_small_lon_lat,
#' K_map = K_small_lon_lat,
#' r = log(2),
#' rate = 1 / 1e3
#' )
#'
#' # example without progress bar and messages
#' sim_data_4 <- initialise(
#' n1_map = n1_small, K_map = K_small, K_sd = 5, r = log(5),
#' r_sd = 4, growth = "ricker", rate = 1 / 200,
#' max_dist = 5000, dens_dep = "K2N", progress_bar = FALSE, quiet = TRUE
#' )
#' }
#'
#'
#' @srrstats {G1.4} uses roxygen documentation
#' @srrstats {G2.0a} documented lengths expectation
#' @srrstats {G2.1a, G2.3, G2.3b, SP2.6} documented types expectation
#' @srrstats {SP1.0} specified domain of applicability
#' @srrstats {SP1.1} documented dimensional domain of applicability
#' @srrstats {SP2.0, SP2.0b} check if K_map and n1_map is
#' [`SpatRaster`][terra::SpatRaster-class] and error otherwise
#' @srrstats {SP2.3} load data in spatial formats
#' @srrstats {SP2.8, SP2.9} simple pre-processing routine that validates and
#' transforms input data while maintaining necessary metadata
#' @srrstats {SP4.0, SP4.0b} returns sim_data object
#' @srrstats {SP4.1} returned object has the same unit as the input
#' @srrstats {SP4.2} returned values are documented
#'
initialise <- function(
n1_map, K_map, K_sd = 0, r, r_sd = 0, growth = "gompertz", A = NA,
dens_dep = c("K2N", "K", "none"), border = c("reprising", "absorbing"),
kernel_fun = "rexp", ..., max_dist = NA, calculate_dist = TRUE,
dlist = NULL, progress_bar = TRUE, quiet = FALSE) {
# check if session is interactive
if(!interactive()) {
call_names <- names(match.call())
if(!("progress_bar" %in% call_names)) progress_bar <- FALSE
if(!("quiet" %in% call_names)) quiet <- TRUE
}
#' @srrstats {G2.0, G2.2, G2.13} assert input length
#' @srrstats {G2.1, G2.3, G2.3a, G2.6, SP2.7} assert input type
# Validation of arguments
## input maps
assert_that(inherits(K_map, "SpatRaster"))
assert_that(inherits(n1_map, "SpatRaster"))
changing_env <- nlyr(K_map) != 1
K_n1_map_check(K_map, n1_map, changing_env)
## K_sd
assert_that(length(K_sd) == 1)
assert_that(is.numeric(K_sd))
assert_that(K_sd >= 0)
## r
assert_that(length(r) == 1)
assert_that(is.numeric(r))
## r_sd
assert_that(length(r_sd) == 1)
assert_that(is.numeric(r_sd))
assert_that(r_sd >= 0)
## growth
assert_that(length(growth) == 1)
assert_that(is.character(growth))
## A
assert_that(length(A) == 1)
assert_that(
is.numeric(A) | is.na(A),
msg = "parameter A can be set either as NA or as a single number")
## dens_dep
dens_dep <- match.arg(dens_dep)
## border
border <- match.arg(border)
## kernel_fun
assert_that(length(kernel_fun) == 1)
assert_that(is.character(kernel_fun))
## max_dist
assert_that(length(max_dist) == 1)
assert_that(
(is.numeric(max_dist) && !(max_dist < 0)) || is.na(max_dist),
msg = "parameter max_dist can be set either as NA or as a single positive number") #nolint
## calculate_dist
assert_that(length(calculate_dist) == 1)
assert_that(is.logical(calculate_dist))
## dlist
assert_that(
is.list(dlist) || is.null(dlist),
msg = "parameter dlist can be set either as NULL or as a list with integers") #nolint
## progress_bar
assert_that(length(progress_bar) == 1)
assert_that(is.logical(progress_bar))
## quiet
assert_that(length(quiet) == 1)
assert_that(is.logical(quiet))
#' @srrstats {G2.16} Check for NaNs and convert them to Nas
# classify NaN to NA for input maps
if ((any(is.nan(values(n1_map))) || any(is.nan(values(K_map)))) && !quiet) {
message("NaN values were found in input maps and replaced with NA (cells outside the study area)") #nolint
n1_map <- classify(n1_map, cbind(NaN, NA))
K_map <- classify(K_map, cbind(NaN, NA))
}
# define ncells and id raster
ncells <- ncell(n1_map)
id <- n1_map
values(id) <- matrix(1:ncells, nrow(n1_map), ncol(n1_map))
# define population dynamic function and dispersal kernel
dynamics <- function(x, r, K, A) match.fun(growth)(x, r, K, A)
kernel <- function(n) match.fun(kernel_fun)(n, ...)
# apply environmental stochasticity if specified (space specific)
if (K_sd > 0) {
K_map <- app(K_map, function(x) {
suppressWarnings(rlnorm(length(x), log(x), log(K_sd)))
})
}
# define matrix witch data about each cell grid
data_table <- as.matrix(data.frame(
values(id),
xyFromCell(id, 1:ncells),
K_get_init_values(K_map, changing_env),
values(n1_map)
))
colnames(data_table) <- c("id", "x", "y", "K", "N")
data_table <- data_table[order(data_table[, "id"]), ]
# ids of cells within the study are
id_within <- data_table[!is.na(data_table[, "K"]), "id"]
# bool matrix - the study area
within_mask <- as.matrix(!is.na(n1_map), wide = TRUE)
# check if raster is planar
planar <- !is.lonlat(id)
# define dist_params: dist_bin, dist_resolution, max_avl_dist
if(!planar) {
# lon/lat input maps
# calculate and extract dist_params
dist_params <- calculate_dist_params(id, id_within, data_table, progress_bar, quiet)
dist_bin <- dist_params["dist_bin"]
dist_resolution <- dist_params["dist_resolution"]
max_avl_dist <- dist_params["max_avl_dist"]
} else {
# planar input maps
# get input map resolution
dist_resolution <- res(n1_map)
if(dist_resolution[1] != dist_resolution[2]) {
# not square grid cells
# calculate and extract dist_params
dist_params <- calculate_dist_params(id, ncells, data_table, progress_bar, quiet)
dist_bin <- dist_params["dist_bin"]
dist_resolution <- dist_params["dist_resolution"]
max_avl_dist <- dist_params["max_avl_dist"]
} else {
# square grid cells
# set default dist_params
dist_resolution <- dist_resolution[1]
dist_bin <- 0
max_avl_dist <- NULL
}
}
# define max_dist
max_dist <- ifelse(
is.na(max_dist),
round(quantile(kernel(1e4), 0.99, names = FALSE) /
dist_resolution) * dist_resolution,
max_dist
)
# define dlist
if (is.null(dlist)) {
dlist <- calc_dist(
calculate_dist, id, data_table, id_within, max_dist, dist_resolution,
dist_bin, progress_bar, quiet
)
}
# define number of cells at each distance for "absorbing" border
if(border == "reprising") {
ncells_in_circle <- NULL
} else if (border == "absorbing") {
if(planar) {
# planar input maps
ncells_in_circle <- ncell_in_circle_planar(id, dist_resolution)
}
else {
# lon/lat input maps
ncells_in_circle <- ncell_in_circle_lonlat(
id, dist_resolution, dist_bin, id_within, max_avl_dist, progress_bar, quiet)
}
}
# output list ---------------------------------------------
out <- list(
n1_map = as.matrix(n1_map, wide = TRUE),
id = wrap(id),
dist_bin = dist_bin,
dist_resolution = dist_resolution,
r = r,
r_sd = r_sd,
K_map = wrap(K_map),
K_sd = K_sd,
growth = growth,
A = A,
dynamics = dynamics,
dens_dep = dens_dep,
border = border,
planar = planar,
max_dist = max_dist,
max_avl_dist = max_avl_dist,
kernel_fun = kernel_fun,
kernel = kernel,
dlist = dlist,
data_table = data_table,
id_within = id_within,
within_mask = within_mask,
ncells = ncells,
ncells_in_circle = ncells_in_circle,
changing_env = changing_env,
call = get_initialise_call(match.call())
)
# set class
class(out) <- c("sim_data", class(out))
return(out)
}
# internal functions -----------------------------------------------------------
#' Validating K_map And n1_map
#'
#' This internal function checks if `K_map` and `n1_map` are correct (contain
#' only non-negative values or NAs) and corresponds to each other. If `K_map`
#' has more than one layer, object comparison is performed based on the first
#' layer.
#' In case of any mistake in given data, suitable error message is printed.
#'
#' @inheritParams initialise
#' @param changing_env logical vector of length 1; determines if carrying
#' capacity map is changing during the [`sim`]ulation
#'
#'
#' @srrstats {G1.4a} uses roxygen documentation (internal function)
#' @srrstats {G2.0a} documented lengths expectation
#'
#' @noRd
#'
K_n1_map_check <- function(K_map, n1_map, changing_env) {
# compare n1_map and K_map
compareGeom(n1_map, K_map)
# check NAs placement
assert_that(
ifelse(!changing_env,
all(is.na(values(n1_map)) == is.na(values(K_map))),
all(is.na(values(n1_map)) == is.na(values(subset(K_map, 1))))),
msg = "n1_map and K_map have NA values in different grid cells")
# check if values are non-negative
if (!all(values(n1_map) >= 0, na.rm = TRUE)) {
stop("n1_map can contain only non-negative values or NAs (which will be automatically reclassified to NA)") #nolint
}
if (!all(values(K_map) >= 0, na.rm = TRUE)) {
stop("K_map can contain only non-negative values or NAs (which will be automatically reclassified to NA)") #nolint
}
}
#' Get Carrying Capacity For The First Time Step
#'
#' [K_get_init_values] returns all values from map of carrying capacity
#' in the first time step.
#'
#' @inheritParams K_n1_map_check
#'
#' @return Numeric vector with values of all cells carrying capacity
#' in the first time step.
#'
#'
#' @srrstats {G1.4a} uses roxygen documentation (internal function)
#'
#' @noRd
#'
K_get_init_values <- function(K_map, changing_env) {
if (!changing_env) {
# get values from K_map
K_values <- values(K_map)
} else {
# get values from first K_map if changing environment
K_values <- values(subset(K_map, 1))
}
return(K_values)
}
#' Check For Precalculating Target Cells
#'
#' [calc_dist] checks if target cells should be precalculated
#' and if so calls [dist_list].
#'
#' @param id [`SpatRaster`][terra::SpatRaster-class]; contains all cells ids
#' @param data_table matrix; contains information about all cells in current
#' time points
#' @param dist_resolution integer vector of length 1; dimension of one side of
#' one cell of `id`; in case of an irregular grid or lon/lat raster it is calculated during [`initialisation`][`initialise`]
#' calculated by [`calculate_dist_params`]
#' @param dist_bin numeric vector of length 1 with value `>= 0`; in case of an irregular grid or lon/lat raster it is
#' calculated by [`calculate_dist_params`]
#' is equal to 0 if input maps are planar raster; if input maps are lon/lat it is calculated by
#' rasters, `dist_bin` is calculated by [`calculate_dist_params`]
#' @param id_within numeric vector; ids of cells inside the study area
#' @inheritParams initialise
#'
#' @return List of target cells ids for each target cells in any distance
#' within `max_dist`.
#'
#'
#' @srrstats {G1.4a} uses roxygen documentation (internal function)
#' @srrstats {G2.0a} documented lengths expectation
#'
#' @name calc_dist
#'
#' @noRd
#'
calc_dist <- function(
calculate_dist, id, data_table, id_within, max_dist, dist_resolution,
dist_bin, progress_bar, quiet) {
if (calculate_dist) {
# message for user
if (!quiet) cat("Calculating distances...", "\n")
# calculate dlist
dlist <- dist_list(
id, data_table, id_within, max_dist, dist_resolution,
dist_bin, progress_bar
)
} else {
# set dlist to NULL
dlist <- NULL
}
return(dlist)
}
#' Precalculating Target Cells For Dispersal
#'
#' `dist_list` prepares data for precalculation of target cells ids
#' and then uses [`target_ids`] for calculation
#'
#' @inheritParams calc_dist
#'
#' @return List of available target cells from each cell within the study area,
#' divided by distance
#'
#'
#' @srrstats {G1.4a} uses roxygen documentation (internal function)
#'
#' @noRd
#'
dist_list <- function(
id, data_table, id_within, max_dist, dist_resolution,
dist_bin, progress_bar) {
# prepare data table - add dist column anf fill it with NAs
data <- cbind(data_table[, c("id", "x", "y")], dist = NA)
# specify function and arguments (for clarity)
tfun <- function(x)
target_ids(x, id, data, min_dist_scaled = 1,
max_dist_scaled = max_dist / dist_resolution,
dist_resolution, dist_bin, id_within)
# calculate targets id with or without progress bar
if (!progress_bar) {
pbo <- pboptions(type = "none")
on.exit(pboptions(pbo))
}
out <- pblapply(id_within, tfun)
return(out)
}
#' Calculate distance parameters
#'
#' Calculates `dist_bin`, `dist_resolution` and `max_avl_dist`. These parameters
#' are necessary to dispersal process if input maps have cells in differen shape
#' than squares.
#'
#' @inheritParams calc_dist
#'
#' @return List of available target cells from each cell within the study area,
#'
#' @noRd
#'
calculate_dist_params <- function(id, id_within, data_table, progress_bar, quiet) {
# calculate distance to the closest neighbour from each grid cell
params_fun <- function(x) {
# get neighbours ids
neighbours <- adjacent(id, cells = x, directions = "queen")
neighbours <- neighbours[!is.na(neighbours)]
# get raster with distances from current id
xy <- vect(xyFromCell(id, x))
crs(xy) <- crs(id)
d <- values(distance(id, xy, progress = 0))
# select distances to the neighbours
neighbour_d <- round(d[neighbours])
# return minimum and maximum distance to neighbour and maximum available distance
return(c(min_neighbour = min(neighbour_d), max_neighbour = max(neighbour_d), max_avl_dist = max(d)))
}
# message to the user
if (!quiet) cat("Calculating distance parameters...", "\n")
# calculate dist params with or without progress bar
if (!progress_bar) {
pbo <- pboptions(type = "none")
on.exit(pboptions(pbo))
}
dist_params <- pbvapply(id_within, params_fun, numeric(3))
# calculate dist_resolution - min distance between each neighbours
dist_resolution <- round(min(dist_params["min_neighbour",]))
# calculate dist_bin - half of the maximum distance between each neighbours
dist_bin <- round(max(dist_params["max_neighbour",] / dist_resolution / 2))
# check if dist_bin isn't too small
if (dist_bin < 1) {
stop("Your input maps have too high resolution. Consider using terra::aggregate() to change it.")
}
# calculate max_avl_dist - max distance between any grid cells divided by dist_resolution plus dist_bin
max_avl_dist <- round(max(dist_params["max_avl_dist",]) / dist_resolution) + dist_bin
return(c(dist_bin = dist_bin, dist_resolution = dist_resolution, max_avl_dist = max_avl_dist))
}
#' Count Cells On Every Distance - planar raster
#'
#' This internal function counts how many cells are reachable on each distance
#' from any cells of template `r`. It takes raster's dist_resolution into account.
#'
#' @param template template [`SpatRaster`][terra::SpatRaster-class] object
#' @param dist_resolution parameter calculated by [`calculate_dist_params`] function
#'
#' @return numeric vector; numbers of target cells on every possible distance
#' range
#'
#'
#' @srrstats {G1.4a} uses roxygen documentation (internal function)
#'
#' @noRd
#'
ncell_in_circle_planar <- function(template, dist_resolution) {
# get the resolution
res <- res(template)
# get extents in both dimensions
a <- xmax(template) - xmin(template)
b <- ymax(template) - ymin(template)
# calculate maximum possible distance
max_dist <- round(sqrt(a**2 + b**2))
# prep new extent
e <- ext(
xmin(template) - max_dist, xmax(template) + max_dist,
ymin(template) - max_dist, ymax(template) + max_dist
)
# prep raster to calculate distances
d <- extend(template, e)
# prep the center (as vector)
center_point <- cbind((e[2] - e[1]) / 2 + e[1], (e[4] - e[3]) / 2 + e[3])
center_vect <- vect(center_point)
crs(center_vect) <- crs(template)
# calculate all distances
d <- distance(d, center_vect, progress = 0)
# adjust to dist_resolution
d_cell <- as.matrix(round(d / dist_resolution))
# number of cells on each distance
ncells_in_circle <- tabulate(d_cell)[1:(max_dist / dist_resolution)]
return(ncells_in_circle)
}
#' Count Cells On Every Distance - lon/lat raster
#'
#' This internal function counts how many cells are reachable on each distance
#' from any cells of template `r`. It takes raster's dist_resolution adn dist_bin
#' into account.
#'
#' @param template template [`SpatRaster`][terra::SpatRaster-class] object
#' @param max_avl_dist numeric vector of length 1; max distance between any grid
#' cells divided by dist_resolution plus dist_bin; describes max available
#' distance achievable in given input maps
#' @inheritParams calc_dist
#' @inheritParams initialise
#'
#' @return numeric matrix with number of columns corresponding to id_within
#' and number of rows equal to max_avl_dist; numbers of target cells on
#' every possible distance from each cell;
#'
#'
#'
#' @srrstats {G1.4a} uses roxygen documentation (internal function)
#'
#' @noRd
#'
ncell_in_circle_lonlat <- function(template, dist_resolution, dist_bin, id_within, max_avl_dist, progress_bar, quiet) {
# get extents in both dimensions
xdiff <- xmax(template) - xmin(template)
ydiff <- ymax(template) - ymin(template)
# prep new extent
multiplier <- 2
e <- ext(
xmin(template) - xdiff * multiplier, xmax(template) + xdiff * multiplier,
ymin(template) - ydiff * multiplier, ymax(template) + ydiff * multiplier
)
# prep raster to calculate distances
extended <- extend(template, e)
# calculates how many cells are reachable on each distance from each grid cell
circles_fun <- function(x) {
# get current grid cell location
xy <- vect(xyFromCell(template, x))
crs(xy) <- crs(template)
# calculate distance to each cell in extended raster from current grid cell
d <- distance(extended, xy, progress = 0)
d <- round(values(d / dist_resolution))
# get only distances that are reachable in the input maps
d <- d[d <= max_avl_dist]
# number of cells on each distance
d_table <- tabulate(d)
# available distances
d_avl <- (1:length(d_table))[d_table != 0]
# expand distances by dist_bin
bin_start <- ifelse(d_avl - dist_bin + 1 < 0, 0, d_avl - dist_bin + 1)
bin_stop <- d_avl + dist_bin
ds <- unlist(lapply(seq_len(length(d_avl)), function(x) {
rep(seq(bin_start[x], bin_stop[x]), times = d_table[x])
}))
# number of cells on each expanded distance
circle <- tabulate(ds)[1:max_avl_dist]
return(circle)
}
# message for user
if (!quiet) cat("Calculating number of cells on each distance...\nThis step may take some time. Consider using border = \"reprising\" with lon/lat rasters if possible.", "\n")
# calculate "circles" with or without progress bar
if (!progress_bar) {
pbo <- pboptions(type = "none")
on.exit(pboptions(pbo))
}
circles <- pbvapply(id_within, circles_fun, numeric(max_avl_dist))
# error if NA
if (any(is.na(circles))) {
stop("Your study area is too big for the \"absorbing\" border. Change the border parameter to \"reprising\".")
}
return(circles)
}
# helper function to get initialise call with info about dlist
get_initialise_call <- function(call) {
if ("dlist" %in% names(call)) {
call$dlist <- TRUE
}
return(call)
}