Inventario Mamíferos Alto Fragua-Indiwassi, Agroforestería para la Conservación
Diego J. Lizcano, Angelica Diaz-Pulido et al.
20/5/2020
library(rgdal) # basic maps and coords
library(maptools) # tools to convert
# library(unmarked)
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 GADMhome <- "D:/BoxFiles/"
tnc <- "C:/Users/diego.lizcano/"
##### extra functions
source(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/R/TEAM_code.R", sep=""))
source(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/R/calendar.R", sep=""))
source(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/R/MultiSpeciesSiteOcc_Stan.R", sep=""))
source(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/R/MultiSpeciesDayOcc_Stan.R", sep=""))
##
datos.raw <- read_excel(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/data/Datos_cameratrap_Caqueta_fixed.xlsx", sep=""))predios <- st_read(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/shp/prediosAgroforesteriaPiedemonte.shp", sep=""))Reading layer prediosAgroforesteriaPiedemonte' from data sourceD:Sync_Caqueta.shp’ using driver `ESRI Shapefile’ Simple feature collection with 1792 features and 26 fields geometry type: MULTIPOLYGON dimension: XY bbox: xmin: 1099776 ymin: 617840.4 xmax: 1141460 ymax: 670269.7 projected CRS: MAGNA_SIRGAS_Colombia_West_zone
Reading layer san_miguel' from data sourceD:Sync_Caqueta_miguel.shp’ using driver `ESRI Shapefile’ Simple feature collection with 292 features and 29 fields geometry type: MULTIPOLYGON dimension: XY bbox: xmin: 1089542 ymin: 614395 xmax: 1094099 ymax: 623402.5 projected CRS: MAGNA-SIRGAS / Colombia West zone
predios_sanMiguel <- st_read(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/shp/prediosAgroforesteriaPiedemonte_sanMiguel_2.shp", sep="")) Reading layer prediosAgroforesteriaPiedemonte_sanMiguel_2' from data sourceD:Sync_Caqueta_sanMiguel_2.shp’ using driver `ESRI Shapefile’ Simple feature collection with 2084 features and 26 fields geometry type: MULTIPOLYGON dimension: XY bbox: xmin: 1089542 ymin: 614395 xmax: 1141460 ymax: 670269.7 projected CRS: MAGNA_SIRGAS_Colombia_West_zone
def_2019 <- raster(paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/defore_raster/def_2019.tif", sep=""))
PNN_AltoFragua <- st_read( paste(home, "Box Sync/CodigoR/Biodiv_Caqueta/shp/PNN_AltoFragua.shp", sep=""))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 bbox: xmin: -76.29621 ymin: 1.18821 xmax: -75.88896 ymax: 1.590837 geographic CRS: WGS 84
# casas.sp1 <- as(casas, "Spatial") # Create Spatial* object
# casas.sp <- as(casas.sp1, "SpatialPoints")
# casas.ppp <- as(casas.sp, "ppp") # make ppp
# casas.dist <- distmap(casas.ppp) # make distance map
#
# Lugar del muestreo
# datos.raw
# get the centroid
centercoord<- c(mean(as.numeric( datos.raw$Longitude)),
mean(as.numeric( datos.raw$Latitude)))
elevation<-raster::getData("SRTM",lon=centercoord[1], lat=centercoord[2])
e<-extent (-76.4, -75.7, 1.05, 1.65) # make the extent
elevation.crop<-crop(elevation, e) # cut elevation to small window
## save shp
# p <- as(e, 'SpatialPolygons')
# crs(p) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
# shapefile(p, paste(tnc, "Box Sync/CodigoR/Biodiv_Caqueta/shp/extent.shp", sep=""))
# plot(elevation.crop)
# points(datos.raw$Longitude, datos.raw$Latitude, add=T)
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("Lon", "Lat"),
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")
# mapview (list(predios, camaras) ,
# zcol = list("COB_ACTUAL", "Sampling_Event"),
# cex= list(NA, 0.8),
# legend = list(TRUE, TRUE),
# col.regions = list(topo.colors, heat.colors), # terrain.colors
# alpha = list(0.5, 0.1),
# alpha.regions=list(0.5, 0.1)
# )
# tmap_mode("plot") ## tmap mode set to plotting
# fincas <- tm_shape(predios) + tm_polygons("COB_ACTUAL") +
# # tm_legend(show = FALSE) +
# tm_shape(camaras) +
# tm_symbols(col="red", size = 0.1)
#
# resgu <- tm_shape(san_miguel) + tm_polygons("COB_ACTUAL") +
# tm_shape(camaras) +
# tm_symbols(col="red", size = 0.1)
#
# tmap_arrange(resgu, fincas)
predios_sanMiguel.geo <- st_transform(predios_sanMiguel, "+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
tmap_mode("plot") ## tmap mode set to plotting
# # get fondo de osm
# bb <- c(-76.31, 1.05111 , -75.7, 1.65)
# caqueta_osm1 <- read_osm(bb, type="osm")
# qtm(CBS_osm1)
#
# tm_shape(predios_sanMiguel.geo) + tm_polygons("COB_ACTUAL", colorNA =NULL, border.col = "grey") +
# tm_shape(camaras) + # tm_symbols (col="red", size = 0.25) +
# tm_dots(col = "Sampling_Event", palette = "Dark2", size = 0.25,
# shape = 16, title = "Camara", legend.show = TRUE,
# legend.is.portrait = TRUE,
# legend.z = NA) +
# tm_layout(scale = .5,
# # legend.position = c(.78,.72),
# legend.outside.size = 0.1,
# legend.title.size = 1.6,
# legend.height = 0.9,
# legend.width = 1.5,
# legend.text.size = 1.2,
# legend.hist.size = 0.5) +
# tm_layout(frame=F) + tm_scale_bar()
# use osm data by bounding box
# bb <- c(-76.31, 1.05111 , -75.7, 1.65)
bb <- c(-76.05, 1.45 , -75.97, 1.495)
bb2 <- c(left = -76.05, bottom = 1.445, right = -75.97, top = 1.495)
map <- get_map(c(left = -76.05, bottom = 1.445, right = -75.97, top = 1.495))
# ggmap(map)
# get fondo de 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) # +
# tm_layout(scale = .5,
# # legend.position = c(.78,.72),
# legend.outside.size = 0.1,
# legend.title.size = 1.6,
# legend.height = 0.9,
# legend.width = 1.5,
# legend.text.size = 1.2,
# legend.hist.size = 0.5) +
# tm_layout(frame=F) +
# tm_scale_bar() + tm_style("classic")
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()
# ##### print ambos
# depto_window
# print(dep_map, vp = viewport(0.75, 0.20, width = 0.2, height = 0.2))
# print(col_map, vp = viewport(0.75, 0.418, width = 0.25, height = 0.25))
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))
##### con facets
# tm_shape(predios_sanMiguel.geo) +
# tm_polygons("COB_ACTUAL", palette = "Dark2", colorNA =NULL, border.col = "grey") +
# tm_shape(camaras) + # tm_symbols (col="red", size = 0.25) +
# tm_dots(col = "red", size = 0.25,
# shape = 16, title = "Camara", legend.show = TRUE,
# legend.is.portrait = TRUE,
# legend.z = NA) + tm_facets(by = "Sampling_Event", free.coords = TRUE)
# tm_layout(scale = .5,
# # legend.position = c(.78,.72),
# legend.outside.size = 0.1,
# legend.title.size = 1.6,
# legend.height = 0.9,
# legend.width = 1.5,
# legend.text.size = 1.2,
# legend.hist.size = 0.5) +
# tm_layout(frame=F) + tm_scale_bar()Duración del muestreo
Las trampas cámara permanecieron activas desde mediados de septiembre 2019 hasta comienzos de diciembre 2019.
Calendario fotografias, inicio y finalizacion de cada cámara
Especies registradas
Las especies registradas en todo em muestreo fueron 87 entre aves y mamíferos.
Mamíferos detectados en Parque
Se detectaron 34 especies de mamíferos
| species | events | phothos | RAI | naiveoccu | |
|---|---|---|---|---|---|
| Dasypus novemcinctus | Dasypus novemcinctus | 18 | 18 | 3.33333333333333 | 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 |
| Microsciurus unknown | Microsciurus unknown | 10 | 10 | 1.85185185185185 | 0.333333333333333 |
| Cuniculus paca | Cuniculus paca | 9 | 9 | 1.66666666666667 | 0.5 |
| Dasypus kappleri | Dasypus kappleri | 9 | 9 | 1.66666666666667 | 0.666666666666667 |
| Leopardus pardalis | Leopardus pardalis | 6 | 6 | 1.11111111111111 | 0.5 |
| Puma concolor | Puma concolor | 4 | 4 | 0.740740740740741 | 0.5 |
| Sciurus unknown | Sciurus unknown | 4 | 4 | 0.740740740740741 | 0.333333333333333 |
| Sciurus igniventris | Sciurus igniventris | 4 | 4 | 0.740740740740741 | 0.333333333333333 |
| 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 |
| Microsciurus flaviventer | Microsciurus flaviventer | 2 | 2 | 0.37037037037037 | 0.166666666666667 |
| Procyon cancrivorus | Procyon cancrivorus | 2 | 2 | 0.37037037037037 | 0.166666666666667 |
| Dasypus unknown | Dasypus unknown | 1 | 1 | 0.185185185185185 | 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
############################################################
# Riqueza de especies y acumulación, modelando la ocurrencia y la detectabilidad.
# Este análisis sigue el método de Dorazio et al. (2006).
# 40 minutos en PC de casa
x <- as.matrix(row.per.sp)
X1 = as.matrix(row.per.sp) # col.per.sp por dias y row.per.sp por sitios (camaras)
nrepls = 90 #dias
especies = MultiSpeciesSiteOcc(nrepls=nrepls, x=X1) #### run stan
## Posterior computed in 2.40911603371302 minutes## Inference for Stan model: 9ba1cd7743adcbf9a81658d6760a5218.
## 4 chains, each with iter=5000; warmup=2500; thin=1;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
## alpha 0.67 0.01 0.46 -0.02 0.37 0.61 0.90 1.76 1596
## beta -4.10 0.00 0.21 -4.53 -4.23 -4.09 -3.96 -3.72 2875
## Omega 0.92 0.00 0.06 0.78 0.89 0.93 0.96 0.99 15066
## sigma_uv[1] 0.69 0.02 0.49 0.07 0.33 0.60 0.94 1.89 445
## sigma_uv[2] 0.61 0.00 0.18 0.29 0.49 0.60 0.72 1.01 2148
## rho_uv -0.34 0.01 0.38 -0.90 -0.64 -0.41 -0.11 0.54 1704
## E_N 18.34 0.01 1.11 15.62 17.75 18.58 19.17 19.79 15066
## E_N_2 20.00 NaN 0.00 20.00 20.00 20.00 20.00 20.00 NaN
## lp__ -171.34 0.92 15.70 -196.96 -181.93 -173.57 -162.61 -134.30 289
## Rhat
## alpha 1.00
## beta 1.00
## Omega 1.00
## sigma_uv[1] 1.01
## sigma_uv[2] 1.00
## rho_uv 1.00
## E_N 1.00
## E_N_2 NaN
## lp__ 1.01
##
## Samples were drawn using NUTS(diag_e) at Mon Jun 01 22:49:58 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).# summary(especies$fit$sims.matrix)
# alpha.post = especies$fit$sims.matrix[,"alpha"]
# sigmaU.post = especies$fit$sims.matrix[,"sigma.u"]
# N.post = especies$fit$sims.matrix[,"N"]
# Stan Equivalent
alpha.post = especies$fit@sim$samples[[1]]$alpha
sigmaU.post = especies$fit@sim$samples[[1]]$`sigma_uv[1]`
N.post = especies$fit@sim$samples[[1]]$E_N
nsites = 6
cum_sp<-CumNumSpeciesPresent (nsites=nsites, alpha=alpha.post, sigmaU=sigmaU.post, N=N.post)
# #histogram of posteriors openwinbugs
# hist(especies$fit$sims.matrix[,"N"],breaks = c(16:30), xlab="Number of mammal species", ylab="Relative frecuency", main="")
# abline(v=length(row.per.sp[,1]),col="blue", lty = 2) # -> lines.histogram(*) observadas
# abline(v=median(especies$fit$sims.matrix[,"N"]),col="red", lty = 2) # -> esperadas por detectabilidad
#
# mean(especies$fit$sims.matrix[,"N"])
# median(especies$fit$sims.matrix[,"N"])
hist(especies$fit@sim$samples[[1]]$E_N , xlab="Number of mammal species", ylab="Relative frecuency", main="")
abline(v=length(row.per.sp[,1]),col="blue", lty = 2) # -> lines.histogram(*) observadas
abline(v=median(especies$fit@sim$samples[[1]]$E_N),col="red", lty = 2) # -> esperadas por detectabilidad
## [1] 18.33106## [1] 18.56285# ##Plot in ggplot
#
# sims <- extract(especies$fit);
#
# freq <- rep(0,S);
# N <- sims$E_N_2
# 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[80:90, ];
#
# N_bar_plot <-
# ggplot(data=dat, aes(x=N, y=freq)) +
# geom_bar(stat="identity") +
# scale_x_continuous(name="Mammal species (N)",
# breaks=c(0,34,35,36,37,38,39,40,41)) +
# scale_y_continuous(name="relative frequency",
# breaks=(0.1 * (0:6))) +
# ggtitle("Posterior: Number of Species (N)");
#
# plot(N_bar_plot);
#
# # Posterior Probabilities of Occurrence and Detection
#
# ilogit <- function(u) { return(1 / (1 + exp(-u))); }
#
# df_psi_theta <- data.frame(cbind(ilogit(sims$logit_psi_sim),
# ilogit(sims$logit_theta_sim)));
# colnames(df_psi_theta)[1] <- "psi";
# colnames(df_psi_theta)[2] <- "theta";
#
# psi_density_plot <-
# ggplot(df_psi_theta, aes(x=psi)) +
# geom_line(stat = "density", adjust=1) +
# scale_x_continuous(name="probability of occurrence",
# limits=c(0,1), breaks=(0.2 * (0:5))) +
# scale_y_continuous(name="probability density",
# limits=c(0,5), breaks=(0:4)) +
# ggtitle("Posterior: Occurrence (psi)");
#
# plot(psi_density_plot);
#
# #####################
# #### compare to vegan
# #####################
# library(vegan)
# sac <- specaccum(t(row.per.sp))
# plot(sac)
#
# sp1 <- specaccum(t(row.per.sp), method = "random")
# sp2 <- specaccum(t(row.per.sp), method = "rarefaction")
# sp3 <- specaccum(t(row.per.sp), method = "collector")
#
# # # summary(sp1)
# plot(sp1, ci.type="poly", col="blue", lwd=2, ci.lty=0, ci.col="lightblue")
# boxplot(sp1, col="yellow", add=TRUE, pch="+")
# #
# plot(sp2, ci.type="poly", col="blue", lwd=2, ci.lty=0, ci.col="lightblue") # rarefaction
# #plot(sp2)
# # boxplot(sp1, col="yellow", add=TRUE, pch="+")
#
# plot(sp3, ci.type="poly", col="blue", lwd=2, ci.lty=0, ci.col="lightblue")
# #plot(sp2)
#
# H <- diversity(t(row.per.sp))
# simp <- diversity(t(row.per.sp), "simpson")
# invsimp <- diversity(t(row.per.sp), "inv")
# r.2 <- rarefy(t(row.per.sp), 2)
# alpha <- fisher.alpha(t(row.per.sp))
# pairs(cbind(H, simp, invsimp, r.2, alpha), pch="+", col="blue")
#
# ## Species richness (S) and Pielou's evenness (J):
# S_1 <- specnumber(t(row.per.sp)) ## rowSums(BCI > 0) does the same...
# J <- H/log(S_1)
#
# rarecurve(row.per.sp)#,step = 20, sample = raremax)
# ## Rarefaction
# (raremax <- max(rowSums(row.per.sp)))
# Srare <- rarefy(t(row.per.sp), raremax)
# plot(S_1, Srare, xlab = "Observed No. of Species", ylab = "Rarefied No. of Species")
# abline(0, 1)
# rarecurve(row.per.sp)#, step = 2, sample = raremax, col = "blue", cex = 0.6)
# Especies observadas en azul = 34. Especies esperadas en rojo = 35.1-35.4. El número esperado esta corrregido por la detectabilidad.
############################################################
## Distribucion posterior de la riqueza de especies
############################################################
# Riqueza de especies y acumulación, modelando la ocurrencia y la detectabilidad.
# Este análisis sigue el método de Dorazio et al. (2006).
# 120 minutos en PC de casa
X2 <- as.matrix(col.per.sp)
# X1 = as.matrix(row.per.sp) # col.per.sp por dias y row.per.sp por sitios (camaras)
nrepls = 90 #dias
especies = MultiSpeciesDayOcc(nrepls=nrepls, x=X2) #### run stan
## Posterior computed in 103.743117018541 minutes## Inference for Stan model: 9ba1cd7743adcbf9a81658d6760a5218.
## 4 chains, each with iter=5000; warmup=2500; thin=1;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
## alpha 6.40 1.61 34.08 -0.68 0.36 1.33 3.92 35.71 446
## beta -6.81 0.04 0.37 -7.43 -7.11 -6.84 -6.52 -6.08 82
## Omega 0.92 0.00 0.05 0.78 0.89 0.93 0.96 0.99 1219
## sigma_uv[1] 1.82 0.22 2.63 0.08 0.78 1.29 2.00 7.50 141
## sigma_uv[2] 0.44 0.02 0.21 0.07 0.30 0.44 0.58 0.88 128
## rho_uv -0.06 0.02 0.44 -0.81 -0.41 -0.08 0.28 0.76 423
## E_N 18.34 0.03 1.09 15.61 17.77 18.57 19.15 19.78 1219
## E_N_2 20.00 NaN 0.00 20.00 20.00 20.00 20.00 20.00 NaN
## lp__ -443.04 3.08 23.31 -489.06 -456.40 -443.64 -430.09 -393.84 57
## Rhat
## alpha 1.01
## beta 1.05
## Omega 1.00
## sigma_uv[1] 1.03
## sigma_uv[2] 1.06
## rho_uv 1.01
## E_N 1.00
## E_N_2 NaN
## lp__ 1.08
##
## Samples were drawn using NUTS(diag_e) at Tue Jun 02 00:33:49 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).# summary(especies$fit$sims.matrix)
# alpha.post = especies$fit$sims.matrix[,"alpha"]
# sigmaU.post = especies$fit$sims.matrix[,"sigma.u"]
# N.post = especies$fit$sims.matrix[,"N"]
# Stan Equivalent
alpha.post = especies$fit@sim$samples[[1]]$alpha
sigmaU.post = especies$fit@sim$samples[[1]]$`sigma_uv[1]`
N.post = especies$fit@sim$samples[[1]]$E_N
nsites = 90 #days
cum_sp<-CumNumSpeciesPresent_day (nsites=nsites, alpha=alpha.post, sigmaU=sigmaU.post, N=N.post)
# #histogram of posteriors openwinbugs
# hist(especies$fit$sims.matrix[,"N"],breaks = c(16:30), xlab="Number of mammal species", ylab="Relative frecuency", main="")
# abline(v=length(row.per.sp[,1]),col="blue", lty = 2) # -> lines.histogram(*) observadas
# abline(v=median(especies$fit$sims.matrix[,"N"]),col="red", lty = 2) # -> esperadas por detectabilidad
#
# mean(especies$fit$sims.matrix[,"N"])
# median(especies$fit$sims.matrix[,"N"])
hist(especies$fit@sim$samples[[1]]$E_N , xlab="Number of mammal species", ylab="Relative frecuency", main="")
abline(v=length(row.per.sp[,1]),col="blue", lty = 2) # -> lines.histogram(*) observadas
abline(v=median(especies$fit@sim$samples[[1]]$E_N),col="red", lty = 2) # -> esperadas por detectabilidad
## [1] 18.33864## [1] 18.56444Mamíferos en las cámaras
Por géneros
## Select Mammals
datos.raw.mam.PNN <- datos.raw %>% 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)
g1<- ggmap(map, extent = 'device') + geom_point(data = datos.raw.mam.PNN,
aes(x=Lon, y=Lat),
colour = "red", size = I(2),
alpha = 1/10, na.rm = TRUE) + overlay +
scale_fill_gradient2("Records", low = "yellow", high = "red", midpoint = 30) + # overlay + scale_fill_gradient(low = "blue", high = "red") +
facet_wrap(~ Genus, ncol = 5)
g2<- ggmap(map, extent = 'device') +
stat_density2d(aes(x = Lon, y = Lat, fill = ..level.., alpha = ..level..),
size = 3, bins = 5, data = datos.raw.mam.PNN, geom = "polygon") +
scale_fill_gradient("Events") +
scale_alpha(range = c(.4, .75), guide = FALSE) +
guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10)) +
geom_point(data = datos.raw.mam.PNN,
aes(x=Lon, y=Lat),
colour = "blue", size = I(2),
alpha = 1/10, na.rm = TRUE) +
facet_wrap(~ Genus, ncol = 4)
g3 <- ggmap(map, extent = 'panel') +
geom_point(data = count.mammal,
aes(x = Lon, y = Lat, size = events.length), colour = "red" ) +
facet_wrap (~ sp, ncol = 5, 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
Aves en las cámaras
Se registraron 32 especies de aves.
Por géneros:
## Select Mammals
datos.raw.ave <- filter(datos.raw, Class=="AVES")
##
## 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.ave)
### add camera traps
map <- get_map(c(-76.058971, 1.335207), zoom = 10,
source = 'google', maptype = "terrain", color ="bw")
# ggmap(map) #, fullpage = TRUE)
g1<- ggmap(map, extent = 'device') + geom_point(data = datos.raw.ave,
aes(x=Lon, y=Lat),
colour = "red", size = I(2),
alpha = 1/10, na.rm = TRUE) +
overlay +
scale_fill_gradient(low = "blue", high = "red")
# con google
g3<- g1 + overlay + scale_fill_gradient(low = "blue", high = "red") +
facet_wrap(~ Genus, ncol = 8)
g3References
Información de la sesión en R.
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## attached base packages:
## [1] parallel grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] doBy_4.6.6 rstan_2.19.3 StanHeaders_2.21.0-3
## [4] plyr_1.8.6 reshape2_1.4.4 GADMTools_3.8-1
## [7] classInt_0.4-3 OpenStreetMap_0.3.4 osmdata_0.1.3
## [10] tmaptools_3.0 tmap_3.0 stargazer_5.2.2
## [13] kableExtra_1.1.0 xtable_1.8-4 knitr_1.28
## [16] ggmap_3.0.0 forcats_0.5.0 stringr_1.4.0
## [19] dplyr_1.0.0 purrr_0.3.4 readr_1.3.1
## [22] tidyr_1.1.0 tibble_3.0.1 ggplot2_3.3.1
## [25] tidyverse_1.3.0 readxl_1.3.1 raster_3.1-5
## [28] spatstat_1.64-1 rpart_4.1-15 nlme_3.1-148
## [31] spatstat.data_1.4-3 mapview_2.7.8 sf_0.9-3
## [34] maptools_1.0-1 rgdal_1.5-8 sp_1.4-2
##
## loaded via a namespace (and not attached):
## [1] leafem_0.1.1 colorspace_1.4-1 rjson_0.2.20
## [4] deldir_0.1-25 ellipsis_0.3.0 class_7.3-17
## [7] leaflet_2.0.3 satellite_1.0.2 base64enc_0.1-3
## [10] fs_1.4.1 dichromat_2.0-0 rstudioapi_0.11
## [13] Deriv_4.0 farver_2.0.3 fansi_0.4.1
## [16] lubridate_1.7.8 xml2_1.3.2 codetools_0.2-16
## [19] splines_3.6.1 rosm_0.2.5 polyclip_1.10-0
## [22] jsonlite_1.6.1 rJava_0.9-12 broom_0.5.6
## [25] dbplyr_1.4.3 rgeos_0.5-3 png_0.1-7
## [28] compiler_3.6.1 httr_1.4.1 backports_1.1.6
## [31] assertthat_0.2.1 Matrix_1.2-18 cli_2.0.2
## [34] prettyunits_1.1.1 htmltools_0.4.0 tools_3.6.1
## [37] gtable_0.3.0 glue_1.4.1 Rcpp_1.0.4.6
## [40] cellranger_1.1.0 vctrs_0.3.0 leafsync_0.1.0
## [43] crosstalk_1.1.0.1 lwgeom_0.2-3 xfun_0.13
## [46] ps_1.3.2 rvest_0.3.5 lifecycle_0.2.0
## [49] XML_3.99-0.3 goftest_1.2-2 MASS_7.3-51.6
## [52] scales_1.1.1 hms_0.5.3 spatstat.utils_1.17-0
## [55] inline_0.3.15 RColorBrewer_1.1-2 yaml_2.2.1
## [58] curl_4.3 gridExtra_2.3 loo_2.2.0
## [61] ggspatial_1.1.2 stringi_1.4.6 highr_0.8
## [64] e1071_1.7-3 pkgbuild_1.0.8 matrixStats_0.56.0
## [67] RgoogleMaps_1.4.5.3 rlang_0.4.6 pkgconfig_2.0.3
## [70] bitops_1.0-6 evaluate_0.14 lattice_0.20-41
## [73] tensor_1.5 labeling_0.3 htmlwidgets_1.5.1
## [76] processx_3.4.2 tidyselect_1.1.0 magrittr_1.5
## [79] R6_2.4.1 generics_0.0.2 DBI_1.1.0
## [82] pillar_1.4.4 haven_2.2.0 foreign_0.8-75
## [85] withr_2.2.0 mgcv_1.8-31 units_0.6-6
## [88] stars_0.4-1 abind_1.4-5 modelr_0.1.8
## [91] crayon_1.3.4 KernSmooth_2.23-17 rmarkdown_2.2
## [94] jpeg_0.1-8.1 callr_3.4.3 reprex_0.3.0
## [97] digest_0.6.25 webshot_0.5.2 stats4_3.6.1
## [100] munsell_0.5.0 viridisLite_0.3.0