---
title: "Dung Beetle Inventory"
output: 
  flexdashboard::flex_dashboard:
    storyboard: true
    social: menu
    source: embed
---

```{r setup, include=FALSE}
library(flexdashboard)
```


```{r loadata, include=FALSE, warning=FALSE, message=FALSE}
library(readxl) # lee excel
library(tidyverse)# manipula datos eficientemente
library(sp) # spatial polygon
library(magrittr)
library(raster)
library(rasterVis)
library(sf) # new spatial package
library(mapview) # plotting maps in html
library(knitr)
library(colorspace)
library(DT)
library(leaflet)
library(lubridate)
library(tmap)
library(tmaptools)
library(vegan)
library(mobsim)

###################################
###
### Read Beetle Data
###
###################################

beetle_full <- as.data.frame(read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Tablas_GCS/data/GAICA_inventarios_Escarabajos_Ene_2019.xlsx", sheet = "Full"))

########## Remove NA row
ind_na <- which(is.na(beetle_full$Lat))
# beetle_full <- beetle_full[-ind_na,]

# put year
beetle_full$year <- year(beetle_full$Fecha)

# select by uniques
site_dept_loc_yr <- beetle_full %>% dplyr::select(Lon, Lat, Localidad, Departamento, year) %>%  distinct() # igual a site_dept_loc<-
# beetle_full <- filter(beetle_full, ORDEN == ordenes[1])


# View(site_dept_loc_yr)

```


### Inventories were carried out in farms part of the sustainable cattle ranching project using standardized methods. In this map we depict the farms and dates of sampling.

```{r map}
# make sp spatial point
coordinates(site_dept_loc_yr) = ~Lon+Lat
geo<- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj4string(site_dept_loc_yr) <-  geo

# make sf
sampling.sf = st_as_sf(site_dept_loc_yr)


mapview(sampling.sf, 
        zcol = c("year", "Localidad" ),
         map.types = c("OpenStreetMap.DE", "Esri.WorldImagery" ),
        leafletHeight = 8, 
        burst = TRUE, hide = TRUE)
 


```

***

Five ecoregions were sampled looking for dung beetles using

- Transects across the landscape with pitfall traps 

- The sampling covered forest, pastures and silvopastoral systems

- The main indicators to compare were richness and abundance per region


### Zoom to Meta. Here the samplings were made in ten farms and in years 2013, 2016-2017.

```{r}

sf_por_depto_q<- sampling.sf %>% filter(Departamento == "Meta") 

mapview(sf_por_depto_q, 
        zcol = c("year", "Localidad" ),
         map.types = c("OpenStreetMap.DE", "Esri.WorldImagery" ),
        leafletHeight = 8, 
        burst = TRUE, hide = TRUE)

```

***


- You can change the display options to show yeas or sampled fams.


### See the Species Acumulation Curve for Meta by Year

```{r specacum, warning=FALSE, message=FALSE}
library(vegan)


# select by uniques
site_species <- beetle_full %>% dplyr::select(Lon, Lat, Localidad, Departamento, year,  Fecha, Nombre, Trampa, Conteo) %>%  distinct() # igual a site_dept_loc<-

site_species_q<- site_species %>% filter(Departamento == "Meta")

site_species_q_y1 <- site_species_q %>% filter(year == unique(site_species_q$year)[1]) # yr 2016

site_species_q_y2 <- site_species_q %>% filter(year == unique(site_species_q$year)[2]) # yr 2017

site_species_q_y3 <- site_species_q %>% filter(year == unique(site_species_q$year)[3]) # yr 2017

##### Localidad Trampa
site_species_q_y1$site <- paste(site_species_q_y1$Localidad, 
                             site_species_q_y1$Trampa)
##### Localidad Trampa
site_species_q_y2$site <- paste(site_species_q_y2$Localidad, 
                             site_species_q_y2$Trampa)
##### Localidad Trampa
site_species_q_y3$site <- paste(site_species_q_y3$Localidad, 
                             site_species_q_y3$Trampa)

########  acum by yr1
# get species
sp_qy1 <- unique (site_species_q_y1$Nombre)
# get sites
site_qy1 <- unique (site_species_q_y1$site)

# mat of site by sp

mat_qy1 <- site_species_q_y1 %>%
  group_by(site, Nombre) %>%
  summarize_at("Conteo", sum) %>% 
  spread(Nombre, Conteo, fill=0) %>% as.data.frame() 

########  acum by yr2
# get species
sp_qy2 <- unique (site_species_q_y2$Nombre)
# get sites
site_qy2 <- unique (site_species_q_y2$site)

# mat of site by sp

mat_qy2 <- site_species_q_y2 %>%
  group_by(site, Nombre) %>%
  summarize_at("Conteo", sum) %>% 
  spread(Nombre, Conteo, fill=0) %>% as.data.frame() 

########  acum by yr3
# get species
sp_qy3 <- unique (site_species_q_y3$Nombre)
# get sites
site_qy3 <- unique (site_species_q_y3$site)

# mat of site by sp

mat_qy3 <- site_species_q_y3 %>%
  group_by(site, Nombre) %>%
  summarize_at("Conteo", sum) %>% 
  spread(Nombre, Conteo, fill=0) %>% as.data.frame() 

# Plot by years
par(mfrow = c(1,2))

sp_y2 <- specaccum(mat_qy2[2:dim(mat_qy2)[2]])
plot(sp_y2, ci.type="poly", col="blue", lwd=2, ci.lty=0, ci.col="lightblue", 
     main = paste ("Accumulation curve", unique(site_species_q$year)[2],unique(site_species_q$year)[1]))

sp_y1 <- specaccum(mat_qy1[2:dim(mat_qy1)[2]])
plot(sp_y1, ci.type="poly", col="blue", lwd=2, ci.lty=0, ci.col="lightblue", 
     main = paste ("Accumulation curve", unique(site_species_q$year)[1]), add=T)



sp_y3 <- specaccum(mat_qy3[2:dim(mat_qy3)[2]])
plot(sp_y3, ci.type="poly", col="blue", lwd=2, ci.lty=0, ci.col="lightblue", 
     main = paste ("Accumulation curve", unique(site_species_q$year)[3]))#, add=T)

```

***

Species accumulation curve: species richness

- Total Species found for `r unique(site_species_q$year)[1]` was: `r length( names(mat_qy1))-1`. The species were: `r names(mat_qy1[2:dim(mat_qy1)[2]])`

- Total Species found for `r unique(site_species_q$year)[2]` was: `r length( names(mat_qy2))-1`. The species were: `r names(mat_qy2[2:dim(mat_qy2)[2]])`


### Rank bundance, accumulation and rarefaction in `r unique(site_species_q$year)[1]` for `r unique( site_species_q$Localidad)[1]` and `r unique( site_species_q$Localidad)[2]`

```{r rank_abun_s1_2y1}

mat_q2y1 <- site_species_q_y1 %>%
  group_by(Localidad, Nombre) %>%
  summarize_at("Conteo", sum) %>% 
  spread(Nombre, Conteo, fill=0) %>% as.data.frame() 

############## Site 1
site_species_qs1 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[1])

census=(data.frame(x=site_species_qs1$Lon, y=site_species_qs1$Lat, species=site_species_qs1$Nombre))

sim_dat1 <-  list(census=census, x_min_max=c(min(site_species_qs1$Lon, max(site_species_qs1$Lon))), y_min_max=c(min(site_species_qs1$Lat, max(site_species_qs1$Lat))))

class(sim_dat1) = "community"
abund1 <- community_to_sad(sim_dat1)


############## Site 2
site_species_qs2 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[2])

census=(data.frame(x=site_species_qs2$Lon, y=site_species_qs2$Lat, species=site_species_qs2$Nombre))

sim_dat2 <-  list(census=census, x_min_max=c(min(site_species_qs2$Lon, max(site_species_qs2$Lon))), y_min_max=c(min(site_species_qs2$Lat, max(site_species_qs2$Lat))))

class(sim_dat2) = "community"
abund2 <- community_to_sad(sim_dat2)


############## Site 3
site_species_qs3 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[3])

census=(data.frame(x=site_species_qs3$Lon, y=site_species_qs3$Lat, species=site_species_qs3$Nombre))

sim_dat3 <-  list(census=census, x_min_max=c(min(site_species_qs3$Lon, max(site_species_qs3$Lon))), y_min_max=c(min(site_species_qs3$Lat, max(site_species_qs3$Lat))))

class(sim_dat3) = "community"
abund3 <- community_to_sad(sim_dat3)


############## Site 4
site_species_qs4 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[4])

census=(data.frame(x=site_species_qs4$Lon, y=site_species_qs4$Lat, species=site_species_qs4$Nombre))

sim_dat4 <-  list(census=census, x_min_max=c(min(site_species_qs4$Lon, max(site_species_qs4$Lon))), y_min_max=c(min(site_species_qs4$Lat, max(site_species_qs4$Lat))))

class(sim_dat4) = "community"
abund4 <- community_to_sad(sim_dat4)


############## Site 5
site_species_qs5 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[5])

census=(data.frame(x=site_species_qs5$Lon, y=site_species_qs5$Lat, species=site_species_qs5$Nombre))

sim_dat5 <-  list(census=census, x_min_max=c(min(site_species_qs5$Lon, max(site_species_qs5$Lon))), y_min_max=c(min(site_species_qs5$Lat, max(site_species_qs5$Lat))))

class(sim_dat5) = "community"
abund5 <- community_to_sad(sim_dat5)


#######################################
######### Plot ####################
#######################################3
par(mfrow = c(2,2))
# plot(abund1, method = "rank") 
plot(abund1, method = "rank", sub=unique(site_species_q$Localidad)[1])
spec_curves1 <- spec_sample_curve(sim_dat1, method = c("accumulation", "rarefaction"))
plot(spec_curves1, sub=unique(site_species_q$Localidad)[1])

plot(abund2, method = "rank", sub=unique(site_species_q$Localidad)[2])
spec_curves2 <- spec_sample_curve(sim_dat2, method = c("accumulation", "rarefaction"))
plot(spec_curves2, sub=unique(site_species_q$Localidad)[2])


```

***

Explanatory text here.


### Rank bundance, accumulation and rarefaction in `r unique(site_species_q$year)[1]` for `r unique( site_species_q$Localidad)[3]` and `r unique( site_species_q$Localidad)[4]`

```{r rank_abun_s3_4y1}

#######################################
######### Plot ####################
#######################################3
par(mfrow = c(2,2))
# plot(abund1, method = "rank") 
plot(abund3, method = "rank", sub=unique(site_species_q$Localidad)[3])
spec_curves3 <- spec_sample_curve(sim_dat3, method = c("accumulation", "rarefaction"))
plot(spec_curves3, sub=unique(site_species_q$Localidad)[3])

plot(abund4, method = "rank", sub=unique(site_species_q$Localidad)[4])
spec_curves4 <- spec_sample_curve(sim_dat4, method = c("accumulation", "rarefaction"))
plot(spec_curves4, sub=unique(site_species_q$Localidad)[4])




```

***

Some text here.



### Rank bundance, accumulation and rarefaction in `r unique(site_species_q$year)[1]` for `r unique( site_species_q$Localidad)[5]`

```{r rank_abun_s5y1}

#######################################
######### Plot ####################
#######################################3
par(mfrow = c(1,2))

plot(abund5, method = "rank", sub=unique(site_species_q$Localidad)[5])
spec_curves5 <- spec_sample_curve(sim_dat5, method = c("accumulation", "rarefaction"))
plot(spec_curves5, sub=unique(site_species_q$Localidad)[5])



```



***

Some text here.




### Rank bundance, accumulation and rarefaction in `r unique(site_species_q$year)[2]` for `r unique( site_species_q$Localidad)[1]` and `r unique( site_species_q$Localidad)[2]`

```{r rank_abun_s1_2y2, eval=FALSE}

mat_q2y2 <- site_species_q_y2 %>%
  group_by(Localidad, Nombre) %>%
  summarize_at("Conteo", sum) %>% 
  spread(Nombre, Conteo, fill=0) %>% as.data.frame() 

############## Site 1
site_species_y2s1 <- site_species_q_y2 %>%
  filter(Localidad==mat_q2y1$Localidad[1])

census=(data.frame(x=site_species_y2s1$Lon, y=site_species_y2s1$Lat, species=site_species_y2s1$Nombre))

sim_dat1 <-  list(census=census, x_min_max=c(min(site_species_y2s1$Lon, max(site_species_y2s1$Lon))), y_min_max=c(min(site_species_y2s1$Lat, max(site_species_y2s1$Lat))))

class(sim_dat1) = "community"
abund1 <- community_to_sad(sim_dat1)


############## Site 2
site_species_y2s2 <- site_species_q_y2 %>%
  filter(Localidad==mat_q2y1$Localidad[2])

census=(data.frame(x=site_species_y2s2$Lon, y=site_species_y2s2$Lat, species=site_species_y2s2$Nombre))

sim_dat2 <-  list(census=census, x_min_max=c(min(site_species_y2s2$Lon, max(site_species_y2s2$Lon))), y_min_max=c(min(site_species_y2s2$Lat, max(site_species_y2s2$Lat))))

class(sim_dat2) = "community"
abund2 <- community_to_sad(sim_dat2)


############## Site 3
site_species_y2s3 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[3])

census=(data.frame(x=site_species_y2s3$Lon, y=site_species_y2s3$Lat, species=site_species_y2s3$Nombre))

sim_dat3 <-  list(census=census, x_min_max=c(min(site_species_y2s3$Lon, max(site_species_y2s3$Lon))), y_min_max=c(min(site_species_y2s3$Lat, max(site_species_y2s3$Lat))))

class(sim_dat3) = "community"
abund3 <- community_to_sad(sim_dat3)


############## Site 4
site_species_y2s4 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[4])

census=(data.frame(x=site_species_y2s4$Lon, y=site_species_y2s4$Lat, species=site_species_y2s4$Nombre))

sim_dat4 <-  list(census=census, x_min_max=c(min(site_species_y2s4$Lon, max(site_species_y2s4$Lon))), y_min_max=c(min(site_species_y2s4$Lat, max(site_species_y2s4$Lat))))

class(sim_dat4) = "community"
abund4 <- community_to_sad(sim_dat4)


############## Site 5
site_species_y2s5 <- site_species_q_y1 %>%
  filter(Localidad==mat_q2y1$Localidad[5])

census=(data.frame(x=site_species_y2s5$Lon, y=site_species_y2s5$Lat, species=site_species_y2s5$Nombre))

sim_dat5 <-  list(census=census, x_min_max=c(min(site_species_y2s5$Lon, max(site_species_y2s5$Lon))), y_min_max=c(min(site_species_y2s5$Lat, max(site_species_y2s5$Lat))))

class(sim_dat5) = "community"
abund5 <- community_to_sad(sim_dat5)


#######################################
######### Plot ####################
#######################################3
par(mfrow = c(2,2))
# plot(abund1, method = "rank") 
plot(abund1, method = "rank", sub=unique(site_species_q$Localidad)[1])
spec_curves1 <- spec_sample_curve(sim_dat1, method = c("accumulation", "rarefaction"))
plot(spec_curves1, sub=unique(site_species_q$Localidad)[1])

plot(abund2, method = "rank", sub=unique(site_species_q$Localidad)[2])
spec_curves2 <- spec_sample_curve(sim_dat2, method = c("accumulation", "rarefaction"))
plot(spec_curves2, sub=unique(site_species_q$Localidad)[2])


```

***

Some text here.


### Rank bundance, accumulation and rarefaction in `r unique(site_species_q$year)[2]` for `r unique( site_species_q$Localidad)[3]` and `r unique( site_species_q$Localidad)[4]`

```{r rank_abun_s3_4y2}

#######################################
######### Plot ####################
#######################################3
par(mfrow = c(2,2))
# plot(abund1, method = "rank") 
plot(abund3, method = "rank", sub=unique(site_species_q$Localidad)[3])
spec_curves3 <- spec_sample_curve(sim_dat3, method = c("accumulation", "rarefaction"))
plot(spec_curves3, sub=unique(site_species_q$Localidad)[3])

plot(abund4, method = "rank", sub=unique(site_species_q$Localidad)[4])
spec_curves4 <- spec_sample_curve(sim_dat4, method = c("accumulation", "rarefaction"))
plot(spec_curves4, sub=unique(site_species_q$Localidad)[4])



```

***

Some text here.


### Rank bundance, accumulation and rarefaction in `r unique(site_species_q$year)[2]` for `r unique( site_species_q$Localidad)[5]`

```{r rank_abun_s5y2}

#######################################
######### Plot ####################
#######################################3
par(mfrow = c(1,2))

plot(abund5, method = "rank", sub=unique(site_species_q$Localidad)[5])
spec_curves5 <- spec_sample_curve(sim_dat5, method = c("accumulation", "rarefaction"))
plot(spec_curves5, sub=unique(site_species_q$Localidad)[5])


```

***

aditional text here

### iNEXT plot by Chao et al

```{r inext, eval=FALSE}

## Interpolation and extrapolation of Hill number with order q
library(iNEXT)

mat_q2 <- site_species_q %>%
  group_by(Localidad, Nombre) %>%
  summarize_at("Conteo", sum) %>% 
  spread(Nombre, Conteo, fill=0) %>% as.data.frame() 

BCI.test.no.zero1 <- as.numeric(mat_q2[2:27][1,])
BCI.test.no.zero2 <- as.numeric(mat_q2[2:27][2,])# as.vector(unlist(mat_q[2:27][11:21,]))
BCI.test.no.zero3 <- as.numeric(mat_q2[2:27][3,])
BCI.test.no.zero4 <- as.numeric(mat_q2[2:27][4,])
BCI.test.no.zero5 <- as.numeric(mat_q2[2:27][5,])

i.zero1 <- which(BCI.test.no.zero1 == 0)
BCI.test.no.zero1 <- BCI.test.no.zero1[-i.zero1]

i.zero2 <- which(BCI.test.no.zero2 == 0)
BCI.test.no.zero2 <- BCI.test.no.zero2[-i.zero2]

i.zero3 <- which(BCI.test.no.zero3 == 0)
BCI.test.no.zero3 <- BCI.test.no.zero3[-i.zero3]

i.zero4 <- which(BCI.test.no.zero4 == 0)
BCI.test.no.zero4 <- BCI.test.no.zero4[-i.zero4]

i.zero5 <- which(BCI.test.no.zero5 == 0)
BCI.test.no.zero5 <- BCI.test.no.zero4[-i.zero5]

out_list <- list(# sort(BCI.test.no.zero1, decreasing = TRUE),
                 sort(BCI.test.no.zero2, decreasing = TRUE),
                 sort(BCI.test.no.zero3, decreasing = TRUE),
                 sort(BCI.test.no.zero4, decreasing = TRUE),
                 sort(BCI.test.no.zero5, decreasing = TRUE))

names(out_list) <- c(# "Britania", 
                     "El_Cairo", "El_Cortijo", 
                     "La_Carelia", "Portugal")

out <- iNEXT(out_list, q=c(0), datatype="abundance")

ggiNEXT(out)#, type=1)#, facet.var="site")



```

***

https://github.com/AnneChao/iNEXT

De acuerdo a las curvas de acumulación, los muestreos fueron suficientes y consistentes.