Édité le 2/11/2019 : rajout des paramètres audiomoth.

Édité le 8/11/2019 : reformulations et vérification.

Conditions

Enregistreurs

Disposés en croix à ~ 30 mètres les uns des autres. Coordonnées WGS84 :

library(OpenStreetMap)
library(maptools)
waypt <- tibble(device = 1:9,
                lat = c(47.21193, 47.21207, 47.21223, 47.21200, 47.21213,
                        47.21166, 47.21146, 47.21175, 47.21156),
                lon = c(5.66267, 5.66293, 5.66318, 5.66227, 5.66196, 5.66240,
                        5.66214, 5.66295, 5.66327))

latrange <- extendrange(waypt$lat, f = 3)
lonrange <- extendrange(waypt$lon, f = 3)
ul <- c(latrange[1], lonrange[2])
lr <- c(latrange[2], lonrange[1])
map <- openmap(ul, lr, type = "bing")
map.latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

start <- c(lon = extendrange(waypt$lon, f = 0.3)[1],
           lat = extendrange(waypt$lat, f = 0.3)[1])
end <- gcDestination(lon = start[1], lat = start[2], bearing = 90,
                     dist = 0.03, dist.units = "km", model = "WGS84")

autoplot(map.latlon) +
  geom_label(data = waypt, aes(x = lon, y = lat, label = device)) +
  geom_segment(x = start[1], y = start[2], xend = end[1], yend = end[2], size = 2) +
  geom_label(x = (start[1] + end[1]) / 2, y = (start[2] + end[2]) / 2, label = "30 m") +
  coord_sf(crs = 4326, xlim = extendrange(waypt$lon, f = 0.5),
             ylim = extendrange(waypt$lat, f = 0.5))
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Malgré une précision affichée au GPS (Garmin GPSMAP 66s) de 3 mètres, la précision pour certains points est moins bonnes (~ 5 mètres). En particulier, 1 devrait être plus proche de l’intersection, à 5 mètres au sud-ouest, juste au niveau de la pointe forestière).

library(sf)
## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
st_as_sf(waypt, coords = c("lon", "lat"), crs = 4326) %>%
  st_distance()
## Units: [m]
##           [,1]     [,2]      [,3]     [,4]      [,5]     [,6]      [,7]
##  [1,]  0.00000 25.10367  51.03956 31.28525  58.20054 36.32346  65.89639
##  [2,] 25.10367  0.00000  25.98237 50.60006  73.78400 60.74308  90.44724
##  [3,] 51.03956 25.98237   0.00000 73.52611  93.08662 86.64407 116.34111
##  [4,] 31.28525 50.60006  73.52611  0.00000  27.57491 39.06133  60.83688
##  [5,] 58.20054 73.78400  93.08662 27.57491   0.00000 61.97843  75.72506
##  [6,] 36.32346 60.74308  86.64407 39.06133  61.97843  0.00000  29.70424
##  [7,] 65.89639 90.44724 116.34111 60.83688  75.72506 29.70424   0.00000
##  [8,] 29.16131 35.60824  56.13641 58.53298  86.07754 42.84986  69.31612
##  [9,] 61.30287 62.27527  74.79861 90.17585 117.74581 66.83812  86.32244
##           [,8]      [,9]
##  [1,] 29.16131  61.30287
##  [2,] 35.60824  62.27527
##  [3,] 56.13641  74.79861
##  [4,] 58.53298  90.17585
##  [5,] 86.07754 117.74581
##  [6,] 42.84986  66.83812
##  [7,] 69.31612  86.32244
##  [8,]  0.00000  32.15352
##  [9,] 32.15352   0.00000

Les distances entre enregistreurs au GPS vont de 25 à 36 mètres, avec des distances totales (de 3 à 7 et de 5 à 9) de 116 et 118 mètres.

Mesures

Traitement automatique des enregistrements

On calcule la puissance du signal pour chaque tranche de 10 kHz et pour chaque seconde. On regarde aussi s’il y a de la saturation sur chaque seconde.

library(tuneR)

FFT_wave <- function(sndObj) {
  # Inspired from http://samcarcagno.altervista.org/blog/basic-sound-processing-r/
  # select left audio track
  s1 <- sndObj@left
  # convert sound array to floating point values ranging from -1 to 1
  s1 <- s1 / 2 ^ (sndObj@bit - 1)
  n <- length(s1)
  # Fast Fourier Transform algorithm to get the frequencies
  p <- fft(s1)
  # select just the first half since the second half is a mirror image
  # of the first
  nUniquePts <- ceiling((n + 1) / 2)
  p <- p[1:nUniquePts]
  # take the absolute value, or the magnitude
  p <- abs(p)
  # scale by the number of points so that the magnitude does not
  # depend on the length of the signal or on its sampling frequency
  p <- p / n
  # square it to get the power
  p <- p^2
  # multiply by two (see technical document for details) odd nfft
  # excludes Nyquist point
  if (n %% 2 > 0) {
    p[2:length(p)] <- p[2:length(p)] * 2
  } else {
    p[2:(length(p) - 1)] <- p[2:(length(p) - 1)] * 2
  }
  # create the frequencies array
  freqArray <- (0:(nUniquePts - 1)) * (sndObj@samp.rate / n)
  return(list(freq = freqArray, power = p))
}

contact_length <- 1
dirname <- "recording"
nb_devices <- 9

freq <- list()
# for (device in seq_len(nb_devices)) {
for (device in NULL) {
  files <- list.files(path = paste(dirname, device, sep = "/"),
                      pattern = "*.wav",
                      ignore.case = TRUE)
  for (file in files) {
    filename <- paste(dirname, device, file, sep = "/")
    print(paste("Treating instance", filename))
    sndObj <- readWave(filename)
    # Length of each considered contact
    split_len <- sndObj@samp.rate * contact_length
    iteration <- 0
    while (length(sndObj@left) > 0) {
      cutObj <- sndObj
      cutObj@left <- head(sndObj@left, split_len)
      cutFFT <- FFT_wave(cutObj)
      pow01 <- cutFFT$power[0e4 <= cutFFT$freq & cutFFT$freq < 1e4]
      pow12 <- cutFFT$power[1e4 <= cutFFT$freq & cutFFT$freq < 2e4]
      pow23 <- cutFFT$power[2e4 <= cutFFT$freq & cutFFT$freq < 3e4]
      pow34 <- cutFFT$power[3e4 <= cutFFT$freq & cutFFT$freq < 4e4]
      pow45 <- cutFFT$power[4e4 <= cutFFT$freq & cutFFT$freq < 5e4]
      pow56 <- cutFFT$power[5e4 <= cutFFT$freq & cutFFT$freq < 6e4]
      pow67 <- cutFFT$power[6e4 <= cutFFT$freq & cutFFT$freq < 7e4]
      pow78 <- cutFFT$power[7e4 <= cutFFT$freq & cutFFT$freq < 8e4]
      pow89 <- cutFFT$power[8e4 <= cutFFT$freq & cutFFT$freq < 9e4]
      pow90 <- cutFFT$power[9e4 <= cutFFT$freq & cutFFT$freq < 10e4]
      powRest <- cutFFT$power[10e4 <= cutFFT$freq]
      freq[[length(freq) + 1]] <-
        list(device = device, file = file,
             second = iteration * contact_length,
             sumPow01 = sum(pow01), sumPow12 = sum(pow12),
             sumPow23 = sum(pow23), sumPow34 = sum(pow34),
             sumPow45 = sum(pow45), sumPow56 = sum(pow56),
             sumPow67 = sum(pow67), sumPow78 = sum(pow78),
             sumPow89 = sum(pow89), sumPow90 = sum(pow90),
             sumPowRest = sum(powRest),
             sat = sum(cutObj@left == -2 ^ (cutObj@bit - 1) |
                       cutObj@left == 2 ^ (cutObj@bit - 1) - 1))
      iteration <- iteration + 1
      sndObj@left <- tail(sndObj@left, -split_len)
    }
  }
}
freq <- bind_rows(freq)

# freq <- freq %>%
#   mutate(date = ymd_hms(str_extract(file, "20[0-9]{6}_[0-9]{6}")) +
#            seconds(second)) %>%
#   mutate(device = factor(device,
#                          levels = c(3, 2, 1, 6, 7, 4, 5, 8, 9)))
# 
# write_rds(freq, "intermediate_data.rds", "xz", compression = 9L)
freq <- read_rds("intermediate_data.rds")

Comptons le nombre de secondes saturées par enregistreurs :

freq %>%
  group_by(device) %>%
  summarise(sat = sum(sat != 0), total = n(), percent = sat / n() * 100)
## # A tibble: 9 x 4
##   device   sat total percent
##   <fct>  <int> <int>   <dbl>
## 1 3         18  5293   0.340
## 2 2         32  5349   0.598
## 3 1        113  5406   2.09 
## 4 6        265  5660   4.68 
## 5 7       1741  5695  30.6  
## 6 4        143  5567   2.57 
## 7 5        143  5602   2.55 
## 8 8        150  5471   2.74 
## 9 9        115  5500   2.09

L’enregistreur 3 était bien protégé de la pluie et du trafic de la départementale alors que le 7 était exposé à la pluie, d’où sa saturation plus importante (30 % des cas).

On filtrera par la suite ces saturations qui représentent au total environ 5 % des secondes considérées. Sachant que l’on perd aussi 2 secondes entre chaque pair d’enregistrements successifs, on traite en moyenne 55 secondes par minute (40 secondes seulement pour le 7).

Critère de classification

On va utiliser le critère comparant la puissance sur signal entre 30 et 40 kHz et celle du signal entre 50 et 60 kHz. L’objectif est de déterminer les secondes dans lesquelles des minioptères émettent (leur fréquence d’émission se situe entre 50 et 52 kHz pour le maximum d’intensité, et dans les fréquences supérieures pour le reste du signal).

freq %>%
  filter(date <= make_datetime(2019, 10, 27, 18)) %>%
  filter(sat == 0) %>%
  ggplot(aes(x = sumPow56 / sumPow34, color = factor(device))) +
  stat_ecdf() +
  xlim(0, 2)
## Warning: Removed 414 rows containing non-finite values (stat_ecdf).

On remarque que les caractéristiques des audiomoth varient (pour l’enregistreur 2 notamment).

On va garder seulement les contacts en fixant le seuil à 1.

threshold <- 1
contact <- freq %>%
  filter(sumPow56 / sumPow34 > threshold) %>%
  filter(sat == 0) %>%
  filter(date <= make_datetime(2019, 10, 27, 18))

Vérifions la pertinence du critère pour quelques cas extrêmes. Pour le contact ayant le plus fort signal dans la bande des 50 kHz (ce serait le cas le plus favorable pour obtenir un contact de minioptère) :

library(seewave)
plotSpectro <- function(sample) {
  require(seewave)
  require(tuneR)
  dirname <- "recording"
  contact_length <- 1
  filename <- paste(dirname, sample$device, sample$file, sep = "/")
  sndObj <- readWave(filename)
  split_len <- sndObj@samp.rate * contact_length
  sndObj@left <- head(tail(sndObj@left, -sample$second * split_len), split_len)
  spectro(sndObj, flim = c(0, 100), collevels = seq(-60, 0, 1), main = filename)
}
plotSpectro(contact %>% arrange(desc(sumPow56)) %>% slice(1))

Il s’agit d’un profil de type minioptère qui n’exclue pas les pipistrelles en première analyse (des analyses manuelles confirment la présence de minioptères en particulier pour la seconde 43 de 20191027_165142.WAV sur l’enregistreur 6).

Pour le contact avec le plus faible signal dans la bande des 50 kHz (ce serait le cas le plus défavorable pour obtenir un contact de minioptère) :

plotSpectro(contact %>% arrange(sumPow56) %>% slice(1))

On obtient bien un contact de type minioptère/pipistrelle, ce qui soutient le critère de classification utilisé. Mais cela signifie aussi que des signaux plus faibles (avec moins de différence entres les fréquences des gammes 30 et 50 kHz) ne seront pas conservé bien qu’ils contiennent aussi des contacts de chiroptères. Ce seuil est choisi ainsi pour ne conserver que des contacts clairs.

Pour le contact avec le plus fort signal dans la bande des 40 kHz :

plotSpectro(contact %>% arrange(desc(sumPow45)) %>% slice(1))

Il s’agit ici plutôt d’un profil pipistrelle commune car les fréquences de forte amplitude sont plus basses. On s’intéresse à la quantité de contacts concernés.

contact %>%
  ggplot(aes(x = sumPow56, y = sumPow45, color = device)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10() +
  scale_color_brewer(palette = "Set1")

contact %>% group_by(device) %>% arrange(desc(sumPow45)) %>% slice(1)
## # A tibble: 9 x 16
## # Groups:   device [9]
##   device file  second sumPow01 sumPow12 sumPow23 sumPow34 sumPow45 sumPow56
##   <fct>  <chr>  <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 3      2019…     35  1.98e-3  4.70e-6  1.33e-6  6.80e-7  1.95e-6  2.10e-6
## 2 2      2019…     46  1.63e-3  7.23e-5  3.97e-6  9.07e-7  1.79e-6  1.46e-6
## 3 1      2019…     21  6.71e-5  4.86e-6  1.88e-6  7.12e-7  2.43e-6  2.08e-6
## 4 6      2019…      0  1.84e-6  4.60e-6  3.32e-6  7.08e-7  1.10e-3  1.12e-4
## 5 7      2019…     11  4.41e-6  7.03e-6  1.43e-6  6.40e-7  1.95e-3  1.74e-6
## 6 4      2019…      0  2.73e-6  6.94e-6  4.82e-6  8.29e-7  8.18e-5  7.15e-6
## 7 5      2019…     57  1.85e-6  2.71e-6  1.04e-6  4.86e-7  7.46e-4  4.76e-5
## 8 8      2019…     42  2.44e-6  5.78e-6  4.03e-6  1.38e-6  3.30e-6  3.25e-6
## 9 9      2019…     59  1.90e-6  2.65e-6  1.25e-6  4.82e-7  1.13e-6  3.83e-6
## # … with 7 more variables: sumPow67 <dbl>, sumPow78 <dbl>, sumPow89 <dbl>,
## #   sumPow90 <dbl>, sumPowRest <dbl>, sat <int>, date <dttm>
c(contact %>% filter(sumPow45 > 1e-5) %>% nrow(), nrow(contact))
## [1]  23 734

Une vingtaine de contacts sur les ~700 semblent concerner plutôt des pipistrelles communes et ils sont tous localisés sur les enregistreurs 6, 7, 4 et 5.

Même démarche avec la bande des 30 kHz :

plotSpectro(contact %>% arrange(desc(sumPow34)) %>% slice(1))

On a plutôt un profil de fréquence modulée abrupte, il ne s’agit donc à priori pas d’un minioptère.

contact %>%
  ggplot(aes(x = sumPow56, y = sumPow34, color = device)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10() +
  scale_color_brewer(palette = "Set1")

contact %>% group_by(device) %>% arrange(desc(sumPow34)) %>% slice(1)
## # A tibble: 9 x 16
## # Groups:   device [9]
##   device file  second sumPow01 sumPow12 sumPow23 sumPow34 sumPow45 sumPow56
##   <fct>  <chr>  <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 3      2019…     35  1.98e-3  4.70e-6  1.33e-6  6.80e-7  1.95e-6  2.10e-6
## 2 2      2019…     41  3.72e-6  7.18e-6  2.29e-6  1.03e-6  1.41e-6  1.77e-6
## 3 1      2019…     38  5.79e-3  6.26e-4  1.07e-5  1.05e-6  1.19e-6  1.85e-6
## 4 6      2019…     43  6.78e-4  1.08e-4  3.26e-5  4.05e-6  2.19e-6  1.94e-5
## 5 7      2019…     31  4.02e-5  3.85e-6  1.21e-6  5.32e-6  4.02e-5  1.35e-5
## 6 4      2019…     19  3.18e-5  8.88e-6  6.36e-6  1.97e-6  1.19e-6  2.62e-5
## 7 5      2019…     59  8.52e-3  1.52e-3  1.02e-4  2.73e-6  2.15e-6  5.23e-6
## 8 8      2019…     41  2.46e-6  6.07e-6  3.90e-6  1.44e-6  2.23e-6  2.08e-6
## 9 9      2019…     45  1.61e-6  2.62e-6  1.24e-6  4.91e-7  8.73e-7  2.87e-6
## # … with 7 more variables: sumPow67 <dbl>, sumPow78 <dbl>, sumPow89 <dbl>,
## #   sumPow90 <dbl>, sumPowRest <dbl>, sat <int>, date <dttm>

L’analyse manuelle conduit à incriminer 4 contacts en FM abrupte, tous pour l’enregistreur 7.

Même démarche avec la bande des 20 kHz :

plotSpectro(contact %>% arrange(desc(sumPow23)) %>% slice(1))

Il s’agit ici d’un orthoptère.

contact %>%
  ggplot(aes(x = sumPow56, y = sumPow23, color = device)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10() +
  scale_color_brewer(palette = "Set1")

contact %>% group_by(device) %>% arrange(desc(sumPow23)) %>% slice(1)
## # A tibble: 9 x 16
## # Groups:   device [9]
##   device file  second sumPow01 sumPow12 sumPow23 sumPow34 sumPow45 sumPow56
##   <fct>  <chr>  <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 3      2019…     49  1.25e-5  5.63e-6  1.35e-6  4.74e-7  1.20e-6  6.71e-7
## 2 2      2019…     46  1.63e-3  7.23e-5  3.97e-6  9.07e-7  1.79e-6  1.46e-6
## 3 1      2019…     14  4.59e-3  3.98e-4  1.22e-5  9.44e-7  1.13e-6  1.08e-6
## 4 6      2019…     56  2.09e-6  1.13e-4  4.31e-5  2.86e-6  4.06e-6  3.32e-6
## 5 7      2019…     50  5.56e-6  5.56e-4  2.32e-4  7.98e-7  7.67e-7  9.27e-7
## 6 4      2019…     54  2.11e-6  1.86e-4  8.19e-5  7.29e-7  1.57e-6  7.92e-7
## 7 5      2019…     59  8.52e-3  1.52e-3  1.02e-4  2.73e-6  2.15e-6  5.23e-6
## 8 8      2019…     24  4.01e-6  8.15e-6  4.09e-6  1.32e-6  2.13e-6  2.60e-6
## 9 9      2019…     59  1.90e-6  2.65e-6  1.25e-6  4.82e-7  1.13e-6  3.83e-6
## # … with 7 more variables: sumPow67 <dbl>, sumPow78 <dbl>, sumPow89 <dbl>,
## #   sumPow90 <dbl>, sumPowRest <dbl>, sat <int>, date <dttm>
c(contact %>% filter(sumPow23 > 1e-5 & device == 7) %>% nrow(), nrow(contact))
## [1]  47 734

Les orthoptères apparaissent sur l’enregistreur 7 seulement sur une quarantaine de contacts. Celui-ci est en bordure entre une prairie et une zone arbustive, et le plus éloigné des axes routiers. On peut aussi remarquer qu’un contact survient très tôt par rapport aux autres premiers contacts et qui possède les même caractéristiques (mais ce serait plutôt une goutte d’eau).

contact %>% filter(date < make_datetime(2019, 10, 27, 16, 45))
## # A tibble: 1 x 16
##   device file  second sumPow01 sumPow12 sumPow23 sumPow34 sumPow45 sumPow56
##   <fct>  <chr>  <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 7      2019…     27  7.58e-5  3.34e-5  1.23e-5  9.54e-7  1.21e-6  1.49e-6
## # … with 7 more variables: sumPow67 <dbl>, sumPow78 <dbl>, sumPow89 <dbl>,
## #   sumPow90 <dbl>, sumPowRest <dbl>, sat <int>, date <dttm>

On enlève ces derniers contacts :

contact <- contact %>%
  filter(device != 7 | sumPow23 < 1e-5)

Conclusion sur le critère utilisé pour identifier les contacts de chiroptère:

  • il permet de détecter les contacts les moins bruités
  • il semble concerner essentiellement des minioptères/pipistrelles (> 95 %), quelques pipistrelles communes (~ 5 %) et des contacts marginaux d’autres espèces
  • il semble être précis (peu de faux positifs) : 10 à 20 vérifications manuelles sur des secondes aux caractéristiques extrêmes (celles pour lesquelles le critère est le plus à même d’être inadapté) montre la pertinence du critère

Ses limites :

  • la saturation (liée à la pluie) conduit à ignorer ~5% des secondes
  • les signaux les moins clairs sont ignorés (sans doute beaucoup de faux négatifs)

Analyse des contacts

Comptons le nombre de contacts par enregistreur :

contact %>%
  group_by(device) %>%
  summarise(n = n())
## # A tibble: 9 x 2
##   device     n
##   <fct>  <int>
## 1 3         78
## 2 2         87
## 3 1         76
## 4 6         84
## 5 7         95
## 6 4        122
## 7 5        139
## 8 8          4
## 9 9          2

8 et 9 ne sont pas dans une zone fréquentée. C’est peut-être lié à la proximité avec une route départementale et une zone forestière assez dense.

Analysons le nombre des contacts par minute pour chaque enregistreur :

contact %>%
  mutate(datemin = floor_date(date, "1 mins")) %>%
  group_by(device, datemin) %>%
  summarise(n = n()) %>%
  summarise(n_9 = quantile(n, 0.9), n_max = max(n))
## # A tibble: 9 x 3
##   device   n_9 n_max
##   <fct>  <dbl> <int>
## 1 3       5.20     9
## 2 2       4.9     10
## 3 1       5        7
## 4 6       6        7
## 5 7       5       10
## 6 4       6       11
## 7 5       7.80    10
## 8 8       1.8      2
## 9 9       1        1

On constate qu’il y a rarement plus de 6 secondes de contact par minute (dans 90 % des cas) et jamais plus de 11. Avec une moyenne de 55 secondes traitées par minute, on a rarement plus de 10 % de contact et toujours moins de 20 % (sauf peut-être pour le 7). La densité d’enregistrement est suffisamment faible pour émettre des hypothèses sur la localisation d’un individu en fonction de sa présence sur chaque enregistreur.

Heure de premier contact :

contact %>%
  group_by(device) %>%
  summarise(earliest = min(date))
## # A tibble: 9 x 2
##   device earliest           
##   <fct>  <dttm>             
## 1 3      2019-10-27 16:50:24
## 2 2      2019-10-27 16:50:26
## 3 1      2019-10-27 16:50:31
## 4 6      2019-10-27 16:50:36
## 5 7      2019-10-27 16:50:40
## 6 4      2019-10-27 16:51:26
## 7 5      2019-10-27 16:51:24
## 8 8      2019-10-27 16:59:20
## 9 9      2019-10-27 17:09:13

On constate une succession entre les enregistreurs 3, 2 (2 secondes plus tard), 1 (5 secondes plus tard), 6 (5 secondes plus tard) et 7 (4 secondes plus tard). En considérant une distance de 120 mètres entre l’enregistreur 3 et 7, une durée de transit de 16 secondes et en supposant qu’il s’agit du même individu, la vitesse estimée est de 27 km/h.

Vue globale des contacts :

contact %>%
  ggplot(aes(x = date, y = sumPow56 / sumPow34)) +
  geom_point() +
  facet_wrap(~ device, ncol = 1)

On constate un trafic important de 16h50 à 17h30 puis une période plus calme jusqu’à 18h (UTC).

Faisons un zoom sur les 10 premières minutes :

contact %>%
  ggplot(aes(x = date, y = sumPow56 / sumPow34)) +
  geom_point() +
  geom_vline(xintercept = make_datetime(2019, 10, 27, 16, 50) +
               minutes(0:10)) +
  facet_wrap(~ device, ncol = 1) +
  xlim(c(make_datetime(2019, 10, 27, 16, 50),
         make_datetime(2019, 10, 27, 17, 00)))
## Warning: Removed 516 rows containing missing values (geom_point).

Les successions entre 3, 2, 1, 6 et 7 suggèrent l’existence d’une route de vol. Calculons dans combien de cas des contacts s’enchaînent sur ces enregistreurs sur une durée inférieure à 40 secondes (~10 km/h):

contact_start <- contact %>% filter(device == 3)
flight <- list()
for (d_init in contact_start$date) {
  d <- d_init
  for (dev in c(2, 1, 6, 7)) {
    next_con <- contact %>%
      filter(device == dev & date > d) %>%
      arrange(date) %>%
      select(date) %>%
      slice(1)
    if (nrow(next_con) == 0)
      break
    d <- next_con$date
  }
  if (d == d_init)
    break
  flight[[length(flight) + 1]] <- list(start = as_datetime(d_init),
                                       end = d)
}
flight <- bind_rows(flight)
# Remove each flight with the same ending date
flight <- flight %>%
  group_by(end) %>%
  arrange(start) %>%
  summarise(start = last(start))
flight %>%
  mutate(duration = as.numeric(end - start)) %>%
  filter(duration <= 40) %>%
  mutate(speed = 120 / duration * 3.6) %>%
  mutate(n = n()) %>%
  select(start, duration, speed, n) %>%
  summary()
##      start                        duration        speed             n     
##  Min.   :2019-10-27 16:50:24   Min.   : 9.0   Min.   :13.94   Min.   :16  
##  1st Qu.:2019-10-27 16:55:10   1st Qu.:10.0   1st Qu.:27.00   1st Qu.:16  
##  Median :2019-10-27 16:58:21   Median :12.5   Median :34.62   Median :16  
##  Mean   :2019-10-27 16:59:23   Mean   :14.5   Mean   :33.79   Mean   :16  
##  3rd Qu.:2019-10-27 17:04:06   3rd Qu.:16.0   3rd Qu.:43.20   3rd Qu.:16  
##  Max.   :2019-10-27 17:11:09   Max.   :31.0   Max.   :48.00   Max.   :16

On obtient 16 trajets dans les 21 premières minutes à partir du premier contact et pour des durées de 9 à 31 secondes (vitesse de 14 à 48 km/h).

Il s’agit d’une approche conservative car on observe parfois une succession de contacts sauf pour un des enregistreurs (cas qui sont ici non-comptabilisés).

Une dernière visualisation pour mettre en évidence ces successions :

flight_select <- flight %>%
  mutate(duration = as.numeric(end - start)) %>%
  filter(duration <= 40)
contact %>%
  filter(device %in% c(3, 2, 1, 6, 7)) %>%
  ggplot(aes(x = date, y = factor(device, levels = rev(levels(device))))) +
  geom_segment(data = flight_select, aes(x = start, xend = end), y = 5,
               yend = 1, size = 2, colour = "green") +
  geom_jitter(width = 0, height = 0.1) +
  geom_vline(xintercept = make_datetime(2019, 10, 27, 16, 50) +
               minutes(0:10)) +
  scale_y_discrete("enregistreur") +
  xlim(c(make_datetime(2019, 10, 27, 16, 50),
         make_datetime(2019, 10, 27, 17, 00)))
## Warning: Removed 7 rows containing missing values (geom_segment).
## Warning: Removed 289 rows containing missing values (geom_point).

On suppose qu’au moins 9 individus sont en transit pendant les 10 premières minutes d’activité (il pourrait y en avoir identifier 6 autres).

Prenons l’exemple de la succession de contact entre 16:55:50 et 16:56:10. Il manque un contact en 2 pour avoir une succession parfaite. On peut visualiser l’enregistrement pour tester s’il ne s’agit pas d’un cri non-identifié par le traitement automatique (faux négatif).

(missing <- contact %>%
  select(device, date, file, second) %>%
  filter(device == 3 | device == 1) %>%
  filter(date > make_datetime(2019, 10, 27, 16, 55, 50) &
           date < make_datetime(2019, 10, 27, 16, 56, 10)))
## # A tibble: 5 x 4
##   device date                file                second
##   <fct>  <dttm>              <chr>                <dbl>
## 1 1      2019-10-27 16:55:56 20191027_165550.WAV      6
## 2 1      2019-10-27 16:55:57 20191027_165550.WAV      7
## 3 3      2019-10-27 16:55:51 20191027_165550.WAV      1
## 4 3      2019-10-27 16:55:55 20191027_165550.WAV      5
## 5 3      2019-10-27 16:56:06 20191027_165550.WAV     16
plotSpectro(data.frame(device = 2, file = missing %>% select(file) %>% slice(1),
                       second = 3.5))

Il y a effectivement quelques contacts mais faibles, ce qui explique pourquoi ils n’ont pas été identifiés automatiquement.

Conclusion

Le caractère exploratoire et semi-automatique de cette analyse laisse une légère incertitude sur les espèces contactées et leurs quantités, mais les résultats suggèrent un axe de vol pour au moins une dizaine d’individus (probablement des minioptères) du nord-est vers le sud-ouest au niveau de la route forestière et se poursuivant à la lisière de la prairie.

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] sf_0.7-6        stringr_1.4.0   readr_1.3.1     lubridate_1.7.4
## [5] ggplot2_3.2.0   dplyr_0.8.3    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.1         RColorBrewer_1.1-2 pillar_1.4.2      
##  [4] compiler_3.6.1     class_7.3-15       tools_3.6.1       
##  [7] zeallot_0.1.0      digest_0.6.20      evaluate_0.14     
## [10] tibble_2.1.3       gtable_0.3.0       pkgconfig_2.0.2   
## [13] rlang_0.4.0        cli_1.1.0          DBI_1.0.0         
## [16] yaml_2.2.0         xfun_0.8           e1071_1.7-2       
## [19] withr_2.1.2        knitr_1.23         vctrs_0.2.0       
## [22] hms_0.5.0          classInt_0.3-3     grid_3.6.1        
## [25] tidyselect_0.2.5   glue_1.3.1         R6_2.4.0          
## [28] fansi_0.4.0        rmarkdown_1.14     purrr_0.3.2       
## [31] magrittr_1.5       units_0.6-3        scales_1.0.0      
## [34] backports_1.1.4    htmltools_0.3.6    assertthat_0.2.1  
## [37] colorspace_1.4-1   labeling_0.3       utf8_1.1.4        
## [40] KernSmooth_2.23-16 stringi_1.4.3      lazyeval_0.2.2    
## [43] munsell_0.5.0      lwgeom_0.1-7       crayon_1.3.4