In maudqueroue/intercali: Intercalibration Des Différents Types De Suivis Aeriens

library(testthat)

My function

#' Simulate data and calculate dsm
#'
#' @param map_obj dataframe. La carte utilisée.
#' @param transect_obj dataframe. Les transects utilisés.
#' @param N numeric. Le nombre d'individus dans la zone d'étude
#' @param crs numeric. Le systeme de projection.
#' @param key Character. Forme de la fonction de détection "hn" ou "unif"
#' @param esw_km numeric. Effective strip width (km). Utile que pour la demi normale, sinon NA. Par défaut NA
#' @param g_zero numeric. Le parametre de la loi uniforme. Par défaut NA.
#' @param truncation_m numeric. A partir de quelle distance aucun individu ne peut être détecté.
#' @param adjust character. Adjustment terms to use. "cos" gives cosine (default), "herm" gives Hermite polynomial and "poly" gives simple polynomial. A value of NULL indicates that no adjustments are to be fitted.
#' @param var boolean. TRUE : the variance associated to the abondance estimate is calculated. FALSE : the variance is not estimated. By default = TRUE.
#'
#' @importFrom dplyr filter
#' @importFrom Distance ds
#' @importFrom dsm dsm dummy_ddf
#' @importFrom stats predict
#' @importFrom purrr possibly
#'
#' @return list. Une liste avec toutes les informations utiles pour l'intercalibration des données.
#' @export


sim_and_calculate <- function(map_obj, transect_obj, N, crs, key, esw_km = NA, g_zero = NA, truncation_m, adjust, var = TRUE){


  # simuler les individus
  ind <- simulate_ind(map_obj = map_obj,
                      crs = crs)


  n_ind_simulated <- nrow(ind)


  # calcule les distances
  dist <- calculate_distance(obs_obj = ind, 
                             transect_obj = transect_obj, 
                             crs = crs)

  # fonction de dectection
  dist <- detection(dist_obj = dist,
                    key = key,
                    esw_km = esw_km,
                    g_zero = g_zero,
                    truncation_m = truncation_m)

  n_ind_detected <- dist %>%
    filter(detected == 1) %>%
    nrow()

  # Préparation des données
  list_dsm <- prepare_dsm(map_obj = map_obj,
                          dist_obj = dist, 
                          segs_obj = transect_obj)

  if(is.null(adjust)) {

    # Calcule processus de detection
    if(key == 'hn'){

      detect <- ds(data = list_dsm$dist_dsm,
                   truncation = max(list_dsm$dist_dsm$distance),
                   key = key,
                   adjustment = NULL)

      AIC_ds <- detect$ddf$criterion
    }

    if(key == 'unif'){

      detect <- dummy_ddf(object = list_dsm$obs_dsm$object,
                          size = list_dsm$obs_dsm$size,
                          width = truncation_m,
                          transect = "line")
      AIC_ds <- NA
    }

  }

  else{

    detect <- ds(data = list_dsm$dist_dsm,
                 truncation = max(list_dsm$dist_dsm$distance),
                 key = key,
                 adjustment = adjust)


    AIC_ds <- detect$ddf$criterion

  }


  # Density surface modelling
  # on a que X et Y pour ça
  dsm <- dsm(formula = count~s(X,Y),
             ddf.obj = detect,
             segment.data = list_dsm$segs_dsm,
             observation.data = list_dsm$obs_dsm,
             method="REML")

  # Prediction pour notre carte
  dsm_pred <- predict(object = dsm,
                      newdata = list_dsm$grid_dsm,
                      off.set = list_dsm$grid_dsm$area)

  if(isTRUE(var)) {

    # Récupération variance
    var_dsm <- get_var_dsm(grid_obj = list_dsm$grid_dsm,
                           dsm_obj = dsm)

    CI_2.5 <- var_dsm$CI[1]
    est_mean <- var_dsm$CI[2]
    CI_97.5 <- var_dsm$CI[3]

    in_95CI <- N > CI_2.5 & N < CI_97.5


    if(all(is.null(adjust), key == 'unif')){

      dsm_cv <- var_dsm$cv[1]
      dsm_se <- var_dsm$se[1]
    }

    else{

      dsm_cv <- var_dsm$cv[1, 1]
      dsm_se <- var_dsm$se[1, 1]
    }  

  }


  if(isFALSE(var)){

    CI_2.5 <- NA
    est_mean <- sum(dsm_pred)
    CI_97.5 <- NA

    in_95CI <- NA
    dsm_cv <- NA
    dsm_se <- NA
  }

  list_out <- list(est_mean = est_mean,
                   CI_2.5 = CI_2.5,
                   CI_97.5 = CI_97.5,
                   in_95CI = in_95CI,
                   dsm_cv = dsm_cv,
                   dsm_se = dsm_se,
                   n_ind_simulated = n_ind_simulated,
                   n_ind_detected = n_ind_detected,
                   AIC_ds = AIC_ds,
                   dsm_pred = dsm_pred)


  return(list_out)

}

# Run but keep eval=FALSE to avoid infinite loop
# Execute in the console directly
fusen::inflate(flat_file = "dev/flat_sim_and_calculate.Rmd", 
               vignette_name = "Simulate and Calculate",
               open_vignette = FALSE,
               check = FALSE,
               overwrite = TRUE)

maudqueroue/intercali documentation built on Oct. 8, 2022, 2:09 p.m.