Camera trap assessment of terrestrial mammals from Alto Fragua Indi-Wasi National Park, Colombia
Diego J. Lizcano, Angélica Díaz-Pulido, Javier García, Luisa F. Martinez, Gerlando Delgado-Bermeo, Olber Leonardo Llanos, Arleth Gonzales-Perez, Edwin García Jaramillo, María Fernanda Ordoñez
20/5/2020
yes
Load packages
library (rgdal) # basic maps and coords
library (maptools) # tools to convert
library (sf) # spatial maps
library (mapview) # nice html maps
library (spatstat) # distance maps
library (raster) # raster
library (readxl) # read excel data
library (tidyverse) # Data handling
library (dplyr) # Data frames handling
library (ggmap) # maps in ggplot
library (knitr) # make documents
library (xtable) # tables
library (kableExtra) # Table html
library (stargazer) # tables
library (tmap) # nice plain maps
library (tmaptools) # mas mapas
library (osmdata) # read osm
library (OpenStreetMap) # osm maps
library (grid) # mix maps
library (GADMTools) # subset GADM# Read data
datos.raw <- read.csv("D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/data/mammal_data.csv")PNN_AltoFragua <- st_read( "D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/shp/PNN_AltoFragua.shp")Reading layer PNN_AltoFragua' from data sourceD:Sync_Caqueta_AltoFragua.shp’ using driver `ESRI Shapefile’ Simple feature collection with 1 feature and 28 fields Geometry type: POLYGON Dimension: XY Bounding box: xmin: -76.29621 ymin: 1.18821 xmax: -75.88896 ymax: 1.590837 Geodetic CRS: WGS 84
Lugar del muestreo
datos.raw$Lon <- as.numeric(datos.raw$Longitude)
datos.raw$Lat <- as.numeric(datos.raw$Latitude)
# make a large sf object
datos.raw_sf = st_as_sf(datos.raw, coords = c("Longitude", "Latitude"),
crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
camaras <- datos.raw_sf
camaras_fragua <- camaras %>% filter (Sampling_Event =="PNN Alto Fragua Indi Wasi")
tmap_mode("plot") ## tmap mode set to plotting
# Make bounding boxes
bb <- c(-76.05, 1.45 , -75.97, 1.495) #Caqueta
bb2 <- c(left = -76.05, bottom = 1.445, right = -75.97, top = 1.495)#Fragua
map <- get_map(c(left = -76.05, bottom = 1.445, right = -75.97, top = 1.495))
# get backgroung from osm
fragua_area <- read_osm(bb2, type="stamen-terrain", mergeTiles = TRUE)
colombia <- gadm_sf_loadCountries("COL", level=1, basefile="./")
collimit <- gadm_sf_loadCountries("COL", level=0, basefile="./")
deptos <- gadm_subset(colombia, regions=c("Caquetá"))#, "Huila"))
caqueta_osm1 <- read_osm(bb, type="stamen-terrain", mergeTiles = TRUE) # type puede ser tambien bing, osm
depto_window <- qtm(caqueta_osm1) +
tm_shape(camaras) + # tm_symbols (col="red", size = 0.25) +
tm_dots(col = "Camera_Trap_Name", palette = "Set1", size = 0.5,
shape = 16, title = "Camara", legend.show = TRUE,
legend.is.portrait = TRUE,
legend.z = NA)
dep_map <- tm_shape(deptos$sf, "NAME_1") +
tm_polygons() + tm_text("NAME_1") +
tm_shape(PNN_AltoFragua) +
tm_polygons(col = "green", border.col = "green") +
tm_shape(camaras_fragua) +
tm_symbols(shape = 0, col = "red", size = 0.3)
col_map <- tm_shape(collimit$sf) + tm_polygons() + tm_shape(deptos$sf) + tm_polygons()
ventana <- qtm (fragua_area) +
tm_shape(camaras_fragua) +
tm_dots(col = "red", size = 0.5,
shape = 16, title = "Camara", legend.show = TRUE,
legend.is.portrait = TRUE) + tm_compass(position=c("RIGHT", "TOP")) +
tm_scale_bar() +
tm_graticules(col="blue", alpha=0.2)
ventana
print(dep_map, vp = viewport(0.32, 0.325, width = 0.5, height = 0.5))
print(col_map, vp = viewport(0.21, 0.755, width = 0.38, height = 0.38))
Read Data
| species | events | phothos | RAI | naiveoccu | |
|---|---|---|---|---|---|
| Dasypus novemcinctus | Dasypus novemcinctus | 19 | 19 | 3.51851851851852 | 0.666666666666667 |
| Puma yagouaroundi | Puma yagouaroundi | 10 | 10 | 1.85185185185185 | 0.833333333333333 |
| Pecari tajacu | Pecari tajacu | 10 | 10 | 1.85185185185185 | 0.5 |
| Dasyprocta fuliginosa | Dasyprocta fuliginosa | 10 | 10 | 1.85185185185185 | 0.666666666666667 |
| Dasypus pastasae | Dasypus pastasae | 10 | 10 | 1.85185185185185 | 0.666666666666667 |
| Cuniculus paca | Cuniculus paca | 9 | 9 | 1.66666666666667 | 0.5 |
| Leopardus pardalis | Leopardus pardalis | 6 | 6 | 1.11111111111111 | 0.5 |
| Puma concolor | Puma concolor | 4 | 4 | 0.740740740740741 | 0.5 |
| Cabassous unicinctus | Cabassous unicinctus | 4 | 4 | 0.740740740740741 | 0.333333333333333 |
| Nasua nasua | Nasua nasua | 4 | 4 | 0.740740740740741 | 0.666666666666667 |
| Panthera onca | Panthera onca | 4 | 4 | 0.740740740740741 | 0.666666666666667 |
| Canis lupus familiaris | Canis lupus familiaris | 4 | 4 | 0.740740740740741 | 0.666666666666667 |
| Eira barbara | Eira barbara | 4 | 4 | 0.740740740740741 | 0.333333333333333 |
| Tremarctos ornatus | Tremarctos ornatus | 3 | 3 | 0.555555555555556 | 0.5 |
| Procyon cancrivorus | Procyon cancrivorus | 2 | 2 | 0.37037037037037 | 0.166666666666667 |
RAI es: Relative Abundance Index
Distribución posterior de la riqueza de especies de mamíferos
Por camaras
Riqueza de especies y acumulación, modelando la ocurrencia y la detectabilidad. Este análisis sigue el método de Dorazio et al. (2006).
############################################################
## Distribucion posterior de la riqueza de especies
############################################################
# Species richness and acumulation by combining occurence and detectability
# Following Dorazio et al. (2006). method
# 2 minutes running in a decent laptop
# Get the code
##########
source("D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/R/MultiSpeciesSiteOcc_JAGS.R")
X1 = as.matrix(row.per.sp) # col.per.sp por dias y row.per.sp por sitios (camaras)
nrepls = 90 #dias
especies = MultiSpeciesSiteOcc_JAGS(nrepls=nrepls, X=X1) #### run JAGS##
## Processing function input.......
##
## Done.
##
## Beginning parallel processing using 3 cores. Console output will be suppressed.
##
## Parallel processing completed.
##
## Calculating statistics.......
##
## Done.
## Posterior computed in 1.82755810022354 minutesprint(especies) ## $fit
## JAGS output for model 'MultiSpeciesSiteOccModel', generated by jagsUI.
## Estimates based on 3 chains of 55000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 5000 iterations and thin rate = 30,
## yielding 4998 total samples from the joint posterior.
## MCMC ran in parallel for 1.826 minutes at time 2021-09-17 02:18:18.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## alpha 0.820 0.845 -1.254 0.859 2.346 TRUE 0.903 1.024 121
## beta -4.078 0.281 -4.652 -4.072 -3.507 FALSE 1.000 1.012 173
## rho -0.435 0.459 -0.976 -0.556 0.696 TRUE 0.820 1.008 273
## sigma.u 1.502 1.751 0.264 0.862 6.699 FALSE 1.000 1.068 111
## sigma.v 0.683 0.286 0.355 0.628 1.366 FALSE 1.000 1.026 292
## omega 0.157 0.049 0.087 0.149 0.279 FALSE 1.000 1.013 339
## N 17.356 4.056 15.000 16.000 30.000 FALSE 1.000 1.045 123
## deviance 207.781 15.718 179.705 206.888 240.751 FALSE 1.000 1.000 4998
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 123.5 and DIC = 331.315
## DIC is an estimate of expected predictive error (lower is better).
##
## $data
## $data$n
## [1] 15
##
## $data$nzeroes
## [1] 100
##
## $data$J
## [1] 6
##
## $data$K
## [1] 90
##
## $data$X
## V1 V2 V3 V4 V5 V6
## Puma yagouaroundi 3 2 0 2 2 1
## Puma concolor 1 2 0 1 0 0
## Pecari tajacu 7 0 1 0 0 2
## Cuniculus paca 6 1 2 0 0 0
## Dasypus novemcinctus 5 0 2 11 0 1
## Leopardus pardalis 3 1 2 0 0 0
## Cabassous unicinctus 2 0 0 0 0 2
## Dasyprocta fuliginosa 3 0 2 3 0 2
## Nasua nasua 1 1 1 0 0 1
## Panthera onca 1 1 0 1 0 1
## Dasypus pastasae 2 0 2 2 0 4
## Canis lupus familiaris 1 0 1 1 0 1
## Eira barbara 2 2 0 0 0 0
## Tremarctos ornatus 0 1 0 1 1 0
## Procyon cancrivorus 0 0 0 0 0 2
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
## 0 0 0 0 0 0
##
##
## $X
## V1 V2 V3 V4 V5 V6
## Puma yagouaroundi 3 2 0 2 2 1
## Puma concolor 1 2 0 1 0 0
## Pecari tajacu 7 0 1 0 0 2
## Cuniculus paca 6 1 2 0 0 0
## Dasypus novemcinctus 5 0 2 11 0 1
## Leopardus pardalis 3 1 2 0 0 0
## Cabassous unicinctus 2 0 0 0 0 2
## Dasyprocta fuliginosa 3 0 2 3 0 2
## Nasua nasua 1 1 1 0 0 1
## Panthera onca 1 1 0 1 0 1
## Dasypus pastasae 2 0 2 2 0 4
## Canis lupus familiaris 1 0 1 1 0 1
## Eira barbara 2 2 0 0 0 0
## Tremarctos ornatus 0 1 0 1 1 0
## Procyon cancrivorus 0 0 0 0 0 2sp.mean <- round(mean(especies$fit$sims.list$N), digits=2)
sp.median <- median(especies$fit$sims.list$N)
sp.mode <- as.numeric(names(sort(table(especies$fit$sims.list$N), decreasing=TRUE))[1])
#hist(especies$fit$sims.list$N, breaks=150, xlab="Species Richness", main=paste("main"))#,
alpha.post = especies$fit$sims.list$alpha
sigmaU.post = especies$fit$sims.list$sigma.u
N.post = especies$fit$sims.list$N
nsites = 90
cum_sp <- CumNumSpeciesPresent (nsites=nsites, alpha=alpha.post, sigmaU=sigmaU.post, N=N.post)
# #histogram of posteriors
hist(especies$fit$sims.list$N ,breaks = 85,
xlab="Number of mammal species",
ylab="Relative frecuency", main="",
xlim=c(14,45))
abline(v=length(row.per.sp[,1]),col="blue", lty = 1, lwd=3) # -> lines.histogram(*) observadas
abline(v=sp.median,col="red", lty = 2,lwd=3) # -> esperadas por detectabilidad
abline(v=sp.mode,col="orange", lty = 2,lwd=3) # -> esperadas por detectabilidad
sp.mean #### Media## [1] 17.36sp.median #### Mediana## [1] 16sp.mode #### Moda## [1] 15#######################################################
library(ggplot2);
library(reshape2);
############## ############## ##############
# S is # sp in Super Comunity
############## ############## ##############
S <- max(especies$fit$sims.list$N)
freq <- rep(0,S);
N <- especies$fit$sims.list$N
for (i in 1:length(N))
freq[N[i]] <- freq[N[i]] + 1;
freq <- freq / length(N);
dat <- data.frame(freq);
dat <- cbind(1:S, dat);
colnames(dat)[1] <- "N";
dat <- dat[15:50, ]; ###### adjust here OJO
N_bar_plot <- ggplot(data=dat, aes(x=N, y=freq)) +
geom_bar(stat="identity") +
geom_vline(xintercept = length(row.per.sp[,1]), colour = "blue", linetype = "dashed") + # observed
geom_vline(xintercept = sp.median, colour = "red", linetype = "dashed") + # expected
geom_vline(xintercept = sp.mode, colour = "purple", linetype = "dashed") + # moda
annotate("text", x = 19.5, y=0.08, label = "observed", colour = "blue", size=3 , angle=90) +
annotate("text", x = 25.5, y=0.08, label = "expected", colour = "red", size=3 , angle=90) +
annotate("text", x = 22.5, y=0.085, label = "mode", colour = "purple", size=3 , angle=90) +
scale_x_continuous(name="Mammal species (N)",
breaks=c(20,25,30,35,40,45,50)) +
scale_y_continuous(name="relative frequency",
breaks=(0.1 * (0:6))) +
ggtitle("Posterior: Number of Species (N)");
plot(N_bar_plot)
Especies observadas en azul = 15. Especies esperadas en rojo = 16-18 El número esperado esta corrregido por la detectabilidad.
Mammal per genus
## Select Mammals
##### remove some species.
datos.raw.mam.PNN <- datos.raw.mam %>% filter (Class=="MAMMALIA" & Sampling_Event=="PNN Alto Fragua Indi Wasi")
datos.raw.mam.PNN$sp <- paste(datos.raw.mam.PNN$Genus, datos.raw.mam.PNN$Species, sep=" ")
library("doBy")
count.mammal1<-summaryBy(datos.raw.mam.PNN~sp+Camera_Trap_Name+Lon+Lat+Photo_Date, data=datos.raw.mam.PNN, FUN=length)
count.mammal1$events <- 1
count.mammal<-summaryBy(count.mammal1~sp+Camera_Trap_Name+Lon+Lat+event, data=count.mammal1, FUN=length)
###
## Google key
register_google(key = "AIzaSyDJmcl2kNxd_-qbmajoHobxa6G-RFR_f8s", write = TRUE)
overlay <- stat_density2d(
aes(x = Lon, y = Lat, fill = ..level.., alpha = ..level..),
bins = 10, geom = "polygon",
data = datos.raw.mam.PNN)
### add camera traps
map <- get_map(c(-75.98785, 1.469082), zoom = 14,
source = 'google', maptype = "terrain", color ="bw")
# ggmap(map) #, fullpage = TRUE)
g3 <- ggmap(map, extent = 'panel') +
geom_point(data = count.mammal,
aes(x = Lon, y = Lat, size = events.length), colour = "red" ) +
facet_wrap (~ sp, ncol = 4, labeller = label_bquote(col = italic(.(sp)))) +
# guides(x = guide_axis(angle = 90)) +
scale_size(breaks = c(1,2,3,5,6,7,9,11)) +
# scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
theme_linedraw() + theme(axis.text.x = element_text(angle = 90, hjust = 1))
# theme(legend.position="bottom") +
# guides(x = guide_axis(angle = 90))
ggsave("sp_plot.jpg",plot=g3, width = 30, height = 30, units = "cm", dpi="print")
g3
Sesion Info.
print(sessionInfo(), locale = FALSE)## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] doBy_4.6.10 jagsUI_1.5.2 plyr_1.8.6
## [4] reshape2_1.4.4 GADMTools_3.8-2 classInt_0.4-3
## [7] OpenStreetMap_0.3.4 osmdata_0.1.5 tmaptools_3.1-1
## [10] tmap_3.3-2 stargazer_5.2.2 kableExtra_1.3.4
## [13] xtable_1.8-4 knitr_1.33 ggmap_3.0.0
## [16] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7
## [19] purrr_0.3.4 readr_1.4.0 tidyr_1.1.3
## [22] tibble_3.1.2 ggplot2_3.3.5 tidyverse_1.3.1
## [25] readxl_1.3.1 raster_3.4-13 spatstat_2.2-0
## [28] spatstat.linnet_2.2-1 spatstat.core_2.2-0 rpart_4.1-15
## [31] nlme_3.1-152 spatstat.geom_2.2-0 spatstat.data_2.1-0
## [34] mapview_2.10.0 sf_1.0-0 maptools_1.1-1
## [37] rgdal_1.5-23 sp_1.4-5
##
## loaded via a namespace (and not attached):
## [1] backports_1.2.1 systemfonts_1.0.2 lwgeom_0.2-6
## [4] splines_4.0.3 crosstalk_1.1.1 leaflet_2.0.4.1
## [7] digest_0.6.27 htmltools_0.5.1.1 fansi_0.5.0
## [10] magrittr_2.0.1 tensor_1.5 modelr_0.1.8
## [13] svglite_2.0.0 spatstat.sparse_2.0-0 jpeg_0.1-8.1
## [16] colorspace_2.0-2 rvest_1.0.0 haven_2.4.1
## [19] xfun_0.22 leafem_0.1.6 microbenchmark_1.4-7
## [22] crayon_1.4.1 jsonlite_1.7.2 glue_1.4.2
## [25] stars_0.5-3 polyclip_1.10-0 gtable_0.3.0
## [28] webshot_0.5.2 abind_1.4-5 scales_1.1.1
## [31] DBI_1.1.1 Rcpp_1.0.6 curry_0.1.1
## [34] viridisLite_0.4.0 units_0.7-2 foreign_0.8-81
## [37] proxy_0.4-26 stats4_4.0.3 htmlwidgets_1.5.3
## [40] httr_1.4.2 RColorBrewer_1.1-2 wk_0.4.1
## [43] ellipsis_0.3.2 farver_2.1.0 pkgconfig_2.0.3
## [46] XML_3.99-0.6 rJava_1.0-4 sass_0.4.0
## [49] dbplyr_2.1.1 deldir_0.2-10 utf8_1.2.1
## [52] labeling_0.4.2 tidyselect_1.1.1 rlang_0.4.10
## [55] ggspatial_1.1.5 munsell_0.5.0 cellranger_1.1.0
## [58] tools_4.0.3 cli_2.5.0 generics_0.1.0
## [61] broom_0.7.8 evaluate_0.14 yaml_2.2.1
## [64] goftest_1.2-2 leafsync_0.1.0 fs_1.5.0
## [67] s2_1.0.6 satellite_1.0.2 RgoogleMaps_1.4.5.3
## [70] xml2_1.3.2 compiler_4.0.3 rstudioapi_0.13
## [73] curl_4.3.2 png_0.1-7 e1071_1.7-7
## [76] spatstat.utils_2.2-0 reprex_2.0.0 bslib_0.2.5.1
## [79] stringi_1.6.2 highr_0.9 rgeos_0.5-5
## [82] lattice_0.20-44 Matrix_1.3-4 vctrs_0.3.8
## [85] rosm_0.2.5 pillar_1.6.1 lifecycle_1.0.0
## [88] jquerylib_0.1.4 bitops_1.0-7 R6_2.5.0
## [91] KernSmooth_2.23-20 gridExtra_2.3 rjags_4-10
## [94] codetools_0.2-18 dichromat_2.0-0 MASS_7.3-54
## [97] assertthat_0.2.1 rjson_0.2.20 withr_2.4.2
## [100] Deriv_4.1.3 mgcv_1.8-36 parallel_4.0.3
## [103] hms_1.1.0 coda_0.19-4 class_7.3-19
## [106] rmarkdown_2.9 lubridate_1.7.10 base64enc_0.1-3