Add soils (#178)

* add tests for multiple random seeds.

* add tests for random seeds

* fix error when creating random_seeds that somehow overwrote random_seed in config list.

* add tests to use quarantine_directions.

* fix error in create_random_seeds function where every anthropogenic dispersal was identical.

* remove unused packages from imports

* add tests for temp and precip sd.

* add combination of temp and precip sd correctly

* Add in boolean for using data for initial pests/pathogens in soil.

* add soil activation

* lint and add documentation for random_seeds to pops_model function

* fix notes

* return correct soil_reservoirs from cpp.

* add checks for weather size within r and remove check from cpp

.Rbuildignore

@@ -21,6 +21,7 @@
 ^CHANGELOG\.md$
 ^CODE_OF_CONDUCT\.md$
 ^CONTRIBUTING\.md$
+^.lintr
 
 ^inst/cpp/pops-core/contributing\.md$
 ^inst/cpp/pops-core/LICENSE$

DESCRIPTION

@@ -19,21 +19,20 @@ Authors@R:
              comment = c(ORCID = "0000-0002-5120-5538")))
 Depends: R (>= 4.0.0)
 Imports: 
-    raster (>= 3.5-15),
-    terra (>= 1.5-17),
-    Rcpp (>= 0.12.19),
+    terra,
+    Rcpp,
     stats,
     foreach,
     parallel,
     doParallel,
     landscapemetrics,
     lubridate,
-    sp,
     utils,
     MASS,
-    methods,
     aws.s3,
-    Metrics
+    Metrics,
+    raster,
+    methods
 License: GPL-3 | file LICENSE
 BugReports: https://github.com/ncsu-landscape-dynamics/rpops/issues
 URL: http://www.github.com/ncsu-landscape-dynamics/rpops
@@ -44,6 +43,7 @@ SystemRequirements:
 Suggests: 
     testthat,
     rgdal,
+    sp,
     rlist,
     knitr,
     rmarkdown,

NAMESPACE

@@ -62,5 +62,6 @@ importFrom(terra,vect)
 importFrom(terra,xres)
 importFrom(terra,yres)
 importFrom(utils,read.csv)
+importFrom(utils,read.table)
 importFrom(utils,write.csv)
 useDynLib(PoPS, .registration = TRUE)

R/RcppExports.R

@@ -1,8 +1,8 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-pops_model_cpp <- function(random_seed, multiple_random_seeds, random_seeds, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, quarantine_directions, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, survival_rates, weather_coefficient, weather_coefficient_sd, bbox, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate = 0.0, mortality_time_lag = 2L, start_date = "2018-01-01", end_date = "2018-12-31", treatment_method = "ratio", natural_kernel_type = "cauchy", anthropogenic_kernel_type = "cauchy", percent_natural_dispersal = 0.0, natural_distance_scale = 21, anthropogenic_distance_scale = 0.0, natural_dir = "NONE", natural_kappa = 0, anthropogenic_dir = "NONE", anthropogenic_kappa = 0, frequencies_n_config = NULL, model_type_ = "SI", latency_period = 0L, establishment_probability = 0, dispersal_percentage = 0.99, survival_rate_month = 0L, survival_rate_day = 0L, overpopulation_config = NULL, network_config = NULL, network_data_config = NULL, weather_size = 0L, weather_type = "deterministic", dispersers_to_soils_percentage = 0) {
-    .Call(`_PoPS_pops_model_cpp`, random_seed, multiple_random_seeds, random_seeds, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, quarantine_directions, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, survival_rates, weather_coefficient, weather_coefficient_sd, bbox, res, rows_cols, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate, mortality_time_lag, start_date, end_date, treatment_method, natural_kernel_type, anthropogenic_kernel_type, percent_natural_dispersal, natural_distance_scale, anthropogenic_distance_scale, natural_dir, natural_kappa, anthropogenic_dir, anthropogenic_kappa, frequencies_n_config, model_type_, latency_period, establishment_probability, dispersal_percentage, survival_rate_month, survival_rate_day, overpopulation_config, network_config, network_data_config, weather_size, weather_type, dispersers_to_soils_percentage)
+pops_model_cpp <- function(random_seed, multiple_random_seeds, random_seeds, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, quarantine_directions, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, survival_rates, weather_coefficient, weather_coefficient_sd, bbox, res, rows_cols, soil_reservoirs, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate = 0.0, mortality_time_lag = 2L, start_date = "2018-01-01", end_date = "2018-12-31", treatment_method = "ratio", natural_kernel_type = "cauchy", anthropogenic_kernel_type = "cauchy", percent_natural_dispersal = 0.0, natural_distance_scale = 21, anthropogenic_distance_scale = 0.0, natural_dir = "NONE", natural_kappa = 0, anthropogenic_dir = "NONE", anthropogenic_kappa = 0, frequencies_n_config = NULL, model_type_ = "SI", latency_period = 0L, establishment_probability = 0, dispersal_percentage = 0.99, survival_rate_month = 0L, survival_rate_day = 0L, overpopulation_config = NULL, network_config = NULL, network_data_config = NULL, weather_size = 0L, weather_type = "deterministic", dispersers_to_soils_percentage = 0) {
+    .Call(`_PoPS_pops_model_cpp`, random_seed, multiple_random_seeds, random_seeds, lethal_temperature, lethal_temperature_month, infected, total_exposed, exposed, susceptible, total_populations, total_hosts, mortality_tracker, mortality, quarantine_areas, quarantine_directions, treatment_maps, treatment_dates, pesticide_duration, resistant, movements, movements_dates, temperature, survival_rates, weather_coefficient, weather_coefficient_sd, bbox, res, rows_cols, soil_reservoirs, reproductive_rate, spatial_indices, season_month_start_end, frequency_config, bool_config, mortality_rate, mortality_time_lag, start_date, end_date, treatment_method, natural_kernel_type, anthropogenic_kernel_type, percent_natural_dispersal, natural_distance_scale, anthropogenic_distance_scale, natural_dir, natural_kappa, anthropogenic_dir, anthropogenic_kappa, frequencies_n_config, model_type_, latency_period, establishment_probability, dispersal_percentage, survival_rate_month, survival_rate_day, overpopulation_config, network_config, network_data_config, weather_size, weather_type, dispersers_to_soils_percentage)
 }
 
 # Register entry points for exported C++ functions

R/calibrate.R

@@ -107,7 +107,8 @@
 #' @importFrom lubridate interval time_length mdy %within%
 #' @importFrom MASS mvrnorm
 #' @importFrom Metrics rmse
-#' @importFrom utils write.csv
+#' @importFrom utils write.csv read.table read.csv
+#' @importFrom methods is
 #'
 #' @return a dataframe of the variables saved and their success metrics for
 #' each run
@@ -197,9 +198,10 @@ calibrate <- function(infected_years_file,
                       dispersers_to_soils_percentage = 0,
                       quarantine_directions = "",
                       multiple_random_seeds = FALSE,
-                      random_seeds = NULL,
+                      file_random_seeds = NULL,
                       use_soils = FALSE,
-                      soil_starting_pest_file = "") {
+                      soil_starting_pest_file = "",
+                      start_with_soil_populations = FALSE) {
 
   # add all data to config list
   config <- c()
@@ -290,9 +292,10 @@ calibrate <- function(infected_years_file,
   config$precipitation_coefficient_sd_file <- precipitation_coefficient_sd_file
   config$dispersers_to_soils_percentage <- dispersers_to_soils_percentage
   config$multiple_random_seeds <- multiple_random_seeds
-  config$random_seeds <- random_seeds
+  config$file_random_seeds <- file_random_seeds
   config$use_soils <- use_soils
   config$soil_starting_pest_file <- soil_starting_pest_file
+  config$start_with_soil_populations <- start_with_soil_populations
 
   # call configuration function to perform data checks and transform data into
   # format used in pops c++
@@ -394,6 +397,7 @@ calibrate <- function(infected_years_file,
         reproductive_rate = reproductive_rate,
         spatial_indices = config$spatial_indices,
         season_month_start_end = config$season_month_start_end,
+        soil_reservoirs = config$soil_reservoirs,
         mortality_rate = config$mortality_rate,
         mortality_time_lag = config$mortality_time_lag,
         start_date = config$start_date,
@@ -438,8 +442,8 @@ calibrate <- function(infected_years_file,
         network_movement = config$network_movement,
         weather_size = config$weather_size,
         weather_type = config$weather_type,
-        dispersers_to_soils_percentage = config$dispersers_to_soils_percentage
-      )
+        dispersers_to_soils_percentage = config$dispersers_to_soils_percentage,
+        use_soils = config$use_soils)
       return(data)
     }
 

R/checks.R

@@ -599,8 +599,8 @@ random_seeds_file_checks <- function(x, number_of_iterations = 1) {
   }
 
   if (checks_passed) {
-    random_seeds <- read.csv(x)
-    if (NCOL(random_seeds) != 9 && NROW(random_seeds != number_of_iterations)) {
+    random_seeds <- read.table(x, sep = ",", header = TRUE)
+    if (base::ncol(random_seeds) != 9 || base::nrow(random_seeds) <= number_of_iterations) {
       checks_passed <- FALSE
       failed_check <- random_seeds_dimensions_error
     }

R/configuration.R

@@ -31,12 +31,14 @@ configuration <- function(config) {
   }
 
   if (config$multiple_random_seeds) {
-    if (!is.null(config$random_seeds)) {
+    if (!is.null(config$file_random_seeds)) {
       ## check random seed file
-      random_seeds_file_check <- random_seeds_file_checks(config$random_seeds)
+      random_seeds_file_check <- random_seeds_file_checks(config$file_random_seeds)
       if (!random_seeds_file_check$checks_passed) {
         config$failure <- random_seeds_file_check$failed_check
         return(config)
+      } else {
+        config$random_seeds <- random_seeds_file_check$random_seeds
       }
     } else {
       config$random_seeds <- create_random_seeds(config$number_of_iterations)
@@ -148,9 +150,9 @@ configuration <- function(config) {
     return(config)
   }
 
-  zero_matrix <- infected[[1]]
-  terra::values(zero_matrix) <- 0
-  zero_matrix <- terra::as.matrix(zero_matrix, wide = TRUE)
+  zero_rast <- infected[[1]]
+  terra::values(zero_rast) <- 0
+  zero_matrix <- terra::as.matrix(zero_rast, wide = TRUE)
 
   one_matrix <- infected[[1]]
   terra::values(one_matrix) <- 0
@@ -192,25 +194,33 @@ configuration <- function(config) {
 
   # check that soils raster has the same crs, resolutin, and extent.
   if (config$use_soils) {
-    if (config$function_name %in% aws_bucket_list) {
-      soils_check <-
-        secondary_raster_checks(
-          config$soil_starting_pest_file, infected, config$use_s3, config$bucket)
-    } else {
-      soils_check <- secondary_raster_checks(config$soil_starting_pest_file, infected)
+    config$soil_survival_steps <- ceiling(1 / config$dispersers_to_soils_percentage)
+    soil_reservoirs <- list(zero_matrix)
+    for (sr in 2:(config$soil_survival_steps)) {
+      soil_reservoirs[[sr]] <- zero_matrix
     }
-    if (soils_check$checks_passed) {
-      soil_pests <- soils_check$raster
-      config$soil_pests <- terra::as.matrix(soil_pests, wide = TRUE)
-    } else {
-      config$failure <- soils_check$failed_check
-      if (config$failure == file_exists_error) {
-        config$failure <- detailed_file_exists_error(config$soil_starting_pest_file)
+    if (config$start_with_soil_populations) {
+      if (config$function_name %in% aws_bucket_list) {
+        soils_check <-
+          secondary_raster_checks(
+            config$soil_starting_pest_file, infected, config$use_s3, config$bucket)
+      } else {
+        soils_check <- secondary_raster_checks(config$soil_starting_pest_file, infected)
+      }
+      if (soils_check$checks_passed) {
+        soil_pests <- soils_check$raster
+        soil_reservoirs[[config$soil_survival_steps]] <- terra::as.matrix(soil_pests, wide = TRUE)
+      } else {
+        config$failure <- soils_check$failed_check
+        if (config$failure == file_exists_error) {
+          config$failure <- detailed_file_exists_error(config$soil_starting_pest_file)
+        }
+        return(config)
       }
-      return(config)
     }
+    config$soil_reservoirs <- soil_reservoirs
   } else {
-    config$soil_pests <- zero_matrix
+    config$soil_reservoirs <- list(zero_matrix)
   }
 
   # check that survival_rates raster has the same crs, resolution, and extent
@@ -365,7 +375,10 @@ configuration <- function(config) {
 
       weather_coefficient_stack <- weather_coefficient_stack * precipitation_coefficient
       if (config$weather_type == "probabilistic") {
-        weather_coefficient_sd_stack <- weather_coefficient_sd_stack * precipitation_coefficient_sd
+        # compute sd from combined sd of the two rasters hard coded 10 years as our current
+        weather_coefficient_sd_stack <-
+          combined_sd(temperature_coefficient_sd, precipitation_coefficient_sd,
+                      temperature_coefficient, precipitation_coefficient, 10, 10)
       }
     }
   } else if (config$precip == TRUE) {
@@ -419,11 +432,29 @@ configuration <- function(config) {
 
   if (config$weather == TRUE) {
     config$weather_size <- terra::nlyr(weather_coefficient_stack)
+    if (config$weather_type == "deterministic") {
+      if (config$number_of_time_steps > config$weather_size) {
+        config$failure <- weather_size_deterministic_error
+        return(config)
+      }
+    }
+
     weather_coefficient <- list(terra::as.matrix(weather_coefficient_stack[[1]], wide = TRUE))
     for (i in 2:terra::nlyr(weather_coefficient_stack)) {
       weather_coefficient[[i]] <- terra::as.matrix(weather_coefficient_stack[[i]], wide = TRUE)
     }
+
     if (config$weather_type == "probabilistic") {
+      if (config$number_of_time_steps > config$weather_size) {
+        config$failure <- weather_size_probabilitic_error
+        return(config)
+      }
+
+      if (config$weather_size != terra::nlyr(weather_coefficient_sd_stack)) {
+        config$failure <- weather_sd_layer_error
+        return(config)
+      }
+
       weather_coefficient_sd <-
         list(terra::as.matrix(weather_coefficient_sd_stack[[1]], wide = TRUE))
       for (i in 2:terra::nlyr(weather_coefficient_sd_stack)) {

R/error_messages.R

@@ -146,6 +146,8 @@ infection_years_length_error <- function(num_layers_infected_years, number_of_ti
 
 success_metric_error <- "success_metric is not one of the listed options."
 
+
+
 initial_cond_uncert_error <-
   "use_initial_condition_uncertainty is TRUE but the number of layers in the infected file is not 2.
   This should be a raster file with 2 layers the first being the mean value and the second the
@@ -160,3 +162,15 @@ random_seeds_dimensions_error <-
   set in the model or the number of columns does not equal the number of unique random seeds"
 
 weather_type_error <- "Weather type is not one of 'probabilistic', 'deterministic', or 'none'"
+
+weather_size_deterministic_error <-
+  "Weather coeeficient number of layers with deterministic is not equal to the total number of time
+  steps."
+
+weather_size_probabilitic_error <-
+  "Weather coefficient number of layers with probablisitc is not equal to the total number of time
+  steps annual."
+
+weather_sd_layer_error <-
+  "weather coefficient sd file number of layers not equal to number of layers in weather coefficient
+  file"

R/helpers.R

@@ -231,7 +231,7 @@ create_random_seeds <- function(n) {
   random_seeds <-
     data.frame(disperser_generation = sample(1:999999999, n, replace = FALSE),
                natural_dispersal = sample(1:999999999, n, replace = FALSE),
-               anthropogenic_dispersal = sample(1:999999999999, 1, replace = FALSE),
+               anthropogenic_dispersal = sample(1:999999999999, n, replace = FALSE),
                establishment = sample(1:999999999, n, replace = FALSE),
                weather = sample(1:999999999, n, replace = FALSE),
                movement = sample(1:999999999, n, replace = FALSE),
@@ -264,3 +264,9 @@ output_from_raster_mean_and_sd <- function(x) {
   x2 <- suppressWarnings(terra::app(x, fun))
   return(x2)
 }
+
+# Combine two standard deviation spatRasters
+combined_sd <- function(v1, v2, m1, m2, n1, n2) {
+  (((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 1)) +
+    (((n1 * n2) * (m1 - m2)^2) / ((n1 + n2) * (n1 + n2 - 1)))
+}