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

rm(list=ls())
library(dsims)
library(intercali)
library(ggplot2) 
library(cowplot)
library(Distance)
library(dsm)
library(units)
library(viridis)

Zone d'étude

Je vais réutiliser la carte déjà présente dans le package dssd (plus facile à gérer que celle que j'utilise d'habitude) :

# La carte de la baie de St Andrews
shapefile.name <- system.file("extdata", "StAndrew.shp", package = "dssd")

# Il faut créer l'objet région
region <- make.region(region.name = "St Andrews bay",
                      shape = shapefile.name,
                      units = "m")

plot(region)

Créer des distributions

Je peux générer différentes distributions :

1. Une distribution homogène

density1 <- make.density(region = region,
                         x.space = 500,
                         y.space = 500,
                         constant = 5) # number of animal per m²
plot(density1)

1. Une distribution avec un gradient de distribution :

density2 <- add.hotspot(object = density1,
                        centre = c(-163000, 6245000),
                        sigma = 10000,
                        amplitude = 6)

density2 <- add.hotspot(object = density2,
                        centre = c(-145000, 6275000),
                        sigma = 10000,
                        amplitude = -3)

plot(density2)

1. Une distribution inhomogène dans l'espace créée de manière aléatoire :

density3 <- random_density(region_obj = region,
                          grid_m = 500,
                          density_base = 10,
                          crs = 2154,
                          amplitude = c(-5, 5),
                          sigma = c(2000, 6000),
                          nb_simu = 15)

plot(density3)

On peut ensuite choisir l'abondance d'individus N dans la zone étudiée.

1. Une centaine d'individus

density <- density2

map1 <- extract_map(density_obj = density,
                    N = 100,
                    crs = 2154)

ind1 <- simulate_ind(map_obj = map1,
                     crs = 2154)

plot_obs(obs_obj = ind1,
         map_obj = map1)

1. Cinq fois plus : 500 individus

map2 <- extract_map(density_obj = density,
                    N = 500,
                    crs = 2154)

ind2 <- simulate_ind(map_obj = map2,
                     crs = 2154)

plot_obs(obs_obj = ind2,
         map_obj = map2)

1. Dix fois plus : 1000 individus

map3 <- extract_map(density_obj = density,
                    N = 1000,
                    crs = 2154)

ind3 <- simulate_ind(map_obj = map3,
                     crs = 2154)

plot_obs(obs_obj = ind3,
         map_obj = map3)

Echantillonage

On peut choisir un échantillonnage en zigzag :

transects1 <- create_transect(region_obj = region,
                             crs = 2154,
                             design = "eszigzag",
                             line.length = 400000,
                             design.angle = 30,
                             truncation = 400)

# Plot design de l'étude
plot_transects(transect_obj = transects1, 
               map_obj = region, 
               crs = 2154, 
               ifsegs = FALSE)

Ou un echantillonage parallele :

transects2 <- create_transect(region_obj = region,
                             crs = 2154,
                             design = "systematic",
                             line.length = 400000,
                             design.angle = 30,
                             truncation = 400)

# Plot design de l'étude
plot_transects(transect_obj = transects2, 
               map_obj = region, 
               crs = 2154, 
               ifsegs = FALSE)

On peut controler l'angle des transects et approximativement la longueur totale de transect souhaitée.

transects3 <- create_transect(region_obj = region,
                             crs = 2154,
                             design = "systematic",
                             line.length = 600000,
                             design.angle = 0,
                             truncation = 400)

# Plot design de l'étude
plot_transects(transect_obj = transects3, 
               map_obj = region, 
               crs = 2154, 
               ifsegs = FALSE)

Je continue avec le transect en zigzag

ind <- ind2
map <- map2
transects <- transects1

plot_transects(transect_obj = transects, 
               map_obj = map, 
               crs = 2154, 
               ifsegs = FALSE)

# Crop transects
transects <- crop_transect(transect_obj = transects,
                           map_obj = map)

plot_transects(transect_obj = transects, 
               map_obj = map, 
               crs = 2154, 
               ifsegs = FALSE)

# Segmentize transects
segs <- segmentize_transect(transect_obj = transects, 
                            length_m = 2000, 
                            to = "LINESTRING") 

# Crop segs
segs <- crop_transect(transect_obj = segs,
                      map_obj = map, 
                      ifsegs = TRUE)

# Plot segment
plot_transects(transect_obj = segs, 
               map_obj = map, 
               crs = 2154, 
               ifsegs = TRUE)

Simulation d'individus observés

On commence par calculer la distance entre les invidus et les transects.

# Distance pour chaque individu simulé
dist <- calculate_distance(obs_obj = ind, 
                           transect_obj = segs, 
                           crs = 2154)

Fonction de detection

On peut créer différentes fonctions de détection :

1. Une demi normale pour la détection visuelle.

dist1 <- detection(dist_obj = dist,
                   key = "hn",
                   esw_km = 0.16,
                   truncation_m = 400, 
                   g_zero = 1)


a <- ggplot(dist1, aes(x=distance_m, y=proba)) +
  geom_point(color = "indianred4") +
  xlim(0,500)

b <- ggplot(dist1, aes(x=distance_m, y=detected)) +
  geom_point(color = "indianred4") +
  xlim(0,500)

plot_grid(a,b)

1. On peut changer les paramètres : moins bonne détection, on diminue le esw :

dist2 <- detection(dist_obj = dist,
                   key = "hn",
                   esw_km = 0.05,
                   truncation = 250, 
                   g_zero = 1) 

a <- ggplot(dist2, aes(x=distance_m, y=proba)) +
  geom_point(color = "indianred4") +
  xlim(0,500)

b <- ggplot(dist2, aes(x=distance_m, y=detected)) +
  geom_point(color = "indianred4") +
  xlim(0,500)

plot_grid(a,b)

1. Une uniforme, de paramètre 1 pour la détection numérique.

dist3 <- detection(dist_obj = dist,
                   key = "unif",
                   g_zero = 1,
                   truncation = 250) 

a <- ggplot(dist3, aes(x=distance_m, y=proba)) +
  geom_point(color = "indianred4") +
  xlim(0,500)

b <- ggplot(dist3, aes(x=distance_m, y=detected)) +
  geom_point(color = "indianred4") +
  xlim(0,500)

plot_grid(a,b)

1. Une uniforme de paramètre 0.8 : tout n'est pas détecté :

dist4 <- detection(dist_obj = dist,
                   key = "unif",
                   g_zero = 0.8,
                   truncation = 250) 

a <- ggplot(dist4, aes(x=distance_m, y=proba)) +
  geom_point(color = "indianred4") +
  xlim(0,500) +
  ylim(0,1)

b <- ggplot(dist4, aes(x=distance_m, y=detected)) +
  geom_point(color = "indianred4") +
  xlim(0,500)

plot_grid(a,b)

On peut visualiser qui a été détecté

dist <- dist4

plot(plot_detect(dist_obj = dist, 
                 transect_obj = transects, 
                 map_obj = map, 
                 title = "Détection Visuelle"))

A partir des données simulées on peut calculer l'abondance et la distribution des individus dans la zone (avec Distance et dsm)

#test <- sim_and_calculate(map_obj = map, transect_obj= segs, N = 500, crs = 2154, key = 'unif', esw_km = NA, g_zero = 0.8, truncation_m = 250)

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


# Processus de détetcion
detect <- ds(data = list_dsm$dist_dsm, 
             truncation = max(list_dsm$dist_dsm$distance), 
             key ='unif', 
             adjustment = 'cos')

summary(detect)

detect

# Density surface modelling
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)
sum(dsm_pred)
# Plot
plot_dsm(dsm_pred_obj = dsm_pred,
         map_obj = map)

Estimation de l'abondance et de la 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]

cat(glue::glue("estimation de l'abondance : {round(est_mean)} CI95% {round(CI_2.5)}-{round(CI_97.5)}"))



#### BONUS 

# # Si l'on souhaite faire du soap film smoothing :
# list_dsm$segs_dsm <- list_dsm$segs_dsm %>%
#   rename(x=X) %>%
#   rename(y=Y)
# 
# list_dsm$grid_dsm <- list_dsm$grid_dsm %>%
#   rename(x=X) %>%
#   rename(y=Y)
# 
# # on veut les contours
# contour_coord <- region@region %>%
#   st_sf(crs = 2154) %>%
#   st_coordinates() %>%
#   as.data.frame %>%
#   select("X","Y")
# 
# test <- list(x = contour_coord$X,
#              y = contour_coord$Y)
# 
# soap.knots <- make.soapgrid(test,c(15,10))
# 
# 
# 
# dsm_soap < -dsm(formula = count~s(x, y, bs="so", k=15, xt=list(bnd=list(test))), 
#                 method="REML",
#                 ddf.object = detect, 
#                 segment.data = list_dsm$segs_dsm, 
#                 observation.data = list_dsm$obs_dsm, 
#                 knots = soap.knots)

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