---
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 Guajira Here the samplings were made in five farms and in years 2013 and 2017.

```{r}

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

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 Guajira 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 == "Guajira")

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

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

##### 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)

########  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() 

# Plot by years
par(mfrow = c(1,2))
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]))

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]))#, 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])


```

***

Some 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}

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}

## 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://plot.ly/ggplot2/

If you use ggplot2, `ggplotly()` converts your plots to an interactive, web-based version! It also provides sensible tooltips, which assists decoding of values encoded as visual properties in the plot.

plotly supports some chart types that ggplot2 doesn't (such as 3D surface, point, and line plots). You can create these (or any other plotly) charts using `plot_ly()`.