1 Análisis exploratorio de diversidad

En este apartado se harán los análisis exploratorio de biodiversidad del ictioplancton de los años 1993, 2004, 2005 y 2019.

1.1 Carga de librerías necesarias para los análisis

Se cargan las librerías necesarias para el análisis de EDA de biodiversidad

require(dplyr)
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
require(ggplot2)
## Loading required package: ggplot2
require(reactable)
## Loading required package: reactable
require(iNEXT)
## Loading required package: iNEXT
require(readxl)
## Loading required package: readxl

1.2 Carga de los datos

Se cargan los datos necesarios para los análisis exploratorios.

data_quantity<-readxl::read_excel("../data/raw/biodiversity/PELAGDEMER_Especies.xlsx", sheet="biodiversidad_SNINA")

reactable(data_quantity, defaultPageSize = 50)

Histograma de la ddensidad 

densidad_total<-as.data.frame(colSums(data_quantity[, 10:286]))

colnames(densidad_total)<-c("Total_Density")
# Ordenar por la columna A de manera descendente
densidad_ordenado <- densidad_total %>%
  arrange(desc(Total_Density))

reactable(densidad_ordenado, defaultPageSize = 50)

Histograma de la ddensidad 

#densidad_ordenado$rango<-1:277
summary(densidad_total$Total_Density)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1      22      94    1921     473  120447
boxplot(densidad_total$Total_Density)
Histograma de la ddensidad

Histograma de la ddensidad 

hist(densidad_total$Total_Density,
     freq = TRUE)
Histograma de la ddensidad

Histograma de la ddensidad 

options(scipen=9999)
ggplot()+
  geom_point(aes(y=densidad_ordenado[["Total_Density"]], x=seq(1:277)))+
  scale_y_continuous(trans='log10', limits = function(x){c(min(x), ceiling(100000))})+
  labs(colour = "", subtitle = "Rango - Densidad", title="Abundancia Total",  x ="Rango", y =expression(paste("log(abundancia) [ind 10 ", m^-2, "]")),tag="A.")+
  theme_bw()+
  theme(legend.position="bottom", 
        axis.text.y =  element_text(size = 8),
        axis.title =  element_text(size = 8),
        plot.title = element_text(size = 8),
        plot.subtitle = element_text(size = 8),
        plot.caption= element_text(size = 8),
        plot.tag = element_text(size=8))

Total_Rank_data<- vegan::rad.lognormal(densidad_ordenado$Total_Density)

Total_Rank_radfit<- vegan::radfit(densidad_ordenado$Total_Density)
## Error in glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  : 
##   NA/NaN/Inf in 'x'
ggplot()+
  geom_point(aes(y=Total_Rank_radfit[["y"]], x=seq(1:277)))+
  geom_line(aes(y= Total_Rank_radfit[["models"]][["Null"]][["fitted.values"]],x=seq(1:277),color="Null" ))+
  geom_line(aes(y= Total_Rank_radfit[["models"]][["Preemption"]][["fitted.values"]],x=seq(1:277),color="Preemption"))+
  geom_line(aes(y= Total_Rank_radfit[["models"]][["Lognormal"]][["fitted.values"]],x=seq(1:277),color="Lognormal" ))+
  geom_line(aes(y= Total_Rank_radfit[["models"]][["Zipf"]][["fitted.values"]],x=seq(1:277),color="Zipf"))+
  geom_line(aes(y= Total_Rank_radfit[["models"]][["Mandelbrot"]][["fitted.values"]],x=seq(1:277),color="Mandelbrot" ))+
  scale_y_continuous(trans='log10')+
  scale_color_manual(values = c( "#d7191c","#fdae61", "#a6dba0", "#008837", "#2b83ba"),
                     #breaks=c("Jet Chocó", "Null", "Preemption", "Lognormal", "Zipf", "Mandelbrot"),
                     name = " ",
                     labels=c(  "Mandelbrot", "Zipf","Lognormal","Preemption", "Null")
  )+
  labs(colour = "", subtitle = "Rango - Densidad", title="Densidad total",  x ="Rango", y =expression(paste("log(Densidad) [ind 10 ", m^-2, "]")),tag="B.")+
  theme_bw()+
  theme(legend.position="bottom", 
        axis.text = element_text(size = 8),
        axis.title = element_text(size = 8),
        plot.title = element_text(size = 8),
        plot.subtitle = element_text(size = 8),
        plot.caption= element_text(size = 8),
        plot.tag = element_text(size=8))

#Diversidad####


D0_total<-hilldiv::hill_div(densidad_ordenado, qvalue = 0)
## Registered S3 methods overwritten by 'FSA':
##   method       from
##   confint.boot car 
##   hist.boot    car
D1_total<-hilldiv::hill_div(densidad_ordenado, qvalue = 1)
D2_total<-hilldiv::hill_div(densidad_ordenado, qvalue = 2)



Div_hill_Total<-base::list(D0_total,D1_total,D2_total)
Div_hill_Total_df<-base::as.data.frame(Div_hill_Total)
base::colnames(Div_hill_Total_df)<-c("q0", "q1","q2")

Div_hill_Total_df_trans<-as.data.frame(t(Div_hill_Total_df))
Div_hill_Total_df_trans$Orden<-c(0,1,2)


Rango_DivPlot <- ggplot2::ggplot() +
  geom_line(data = Div_hill_Total_df_trans, aes(x = Orden, y = Total_Density ), color = "black", lwd = 0.5, linetype = 1) +
  geom_point(data = Div_hill_Total_df_trans, aes(x = Orden, y = Total_Density ), color = "black", size =3) +
   labs(
    title = "Perfiles de diversidad",
    x = "Orden",
    y = "Diversidad",
    color = "Variable") +
     theme_minimal()

Rango_DivPlot

reactable(Div_hill_Total_df_trans)

2 Analisis relacionados con el ENOS

2.1 Gráficas Indice ONI

datos_indice_Climaticos <- readxl::read_excel("../data/processed/Indices_Total.xlsx", sheet = "indices")
#datos_indice_Climaticos <- readxl::read_excel("../data/processed/Indices_Total.xlsx", sheet = "crucero_analizados")
ONI <- datos_indice_Climaticos%>%
          subset(index_name == "ONI")



ONI$date <- as.Date(ONI$date)

cruceros<-as.data.frame(as.Date(c("1991-02-01", 
                                  "1991-09-01",
                                  "1991-12-01",
                                  "1993-01-01",
                                  "1993-11-01",
                                  "1994-04-01",
                                  "1994-07-01",
                                  "1994-12-01",
                                  "1995-06-01",
                                  "1996-05-01",
                                  "1996-11-01",
                                  "2008-12-01",
                                  "2009-12-01")))
ceros<-as.data.frame(rep(0, times=13))

cruceros<-cbind(ceros, cruceros)

colnames(cruceros) <-c("value.y", "date")


group.colors <- c(Neutro = "#808080", Niño = "#d7191c", Niña ="#2c7bb6")

ggplot() + 
  geom_col(data = ONI, aes(x = date, y = value, fill = event)) +
  
  scale_fill_manual(values = group.colors) +
  scale_x_date(limits = c(as.Date("1991-01-01"), as.Date("2019-12-31")), 
               date_labels = "%Y", date_breaks = "2 years") +
  geom_hline(yintercept = 0.5, color = "#d7191c", size = 1) +  # Línea roja en 0.5
  geom_hline(yintercept = -0.5, color = "#2c7bb6", size = 1) +  # Línea azul en -0.5
  labs(
    x = "Date", 
    y = "ONI",
    fill = "Evento", 
    title = "Índice Oceánico El Niño (ONI)", 
    subtitle = "3-month rolling mean of ERSST.v5 SST anomalies in the Niño 3.4 region"
  ) +
  theme_bw(base_size = 12) +
  geom_point(data = cruceros, aes(x=date, y =value.y ), size = 2)+
  theme(legend.position = "bottom")

datos_indice_Climaticos_cruceros <- readxl::read_excel("../data/processed/Indices_Total.xlsx", sheet = "pequenos_pelagicos")
datos_indice_Climaticos_cruceros<-datos_indice_Climaticos_cruceros%>%
  select(!index_description)
reactable(datos_indice_Climaticos_cruceros,  defaultPageSize = 20)
nino<-data_quantity%>%
        filter(ONI == "Niño")

neutro<-data_quantity%>%
  filter(ONI == "Neutro")

densidad_nino<-as.data.frame(colSums(nino[, 12:277]))

colnames(densidad_nino)<-c("Niño")

# Ordenar por la columna A de manera descendente
densidad_nino_ord <- densidad_nino %>%
  arrange(desc(Niño))%>%
  filter(Niño != 0)



densidad_neutro<-as.data.frame(colSums(neutro[, 12:277]))

colnames(densidad_neutro)<-c("Neutro")

# Ordenar por la columna A de manera descendente
densidad_neutro_ord <- densidad_neutro %>%
  arrange(desc(Neutro))%>%
  filter(Neutro != 0)

Diversidad_Total_Enos<-cbind(densidad_nino, densidad_neutro)
Diversidad_Total_Enos$Niño<-as.integer(Diversidad_Total_Enos$Niño)

Diversidad_Total_Enos$Neutro<-as.integer(Diversidad_Total_Enos$Neutro)


reactable(Diversidad_Total_Enos, defaultPageSize = 50)
#Diversidad ENOS####

#Niño

D0_nino<-hilldiv::hill_div(densidad_nino_ord, qvalue = 0)
D1_nino<-hilldiv::hill_div(densidad_nino_ord, qvalue = 1)
D2_nino<-hilldiv::hill_div(densidad_nino_ord, qvalue = 2)



Div_hill_nino<-base::list(D0_nino,D1_nino,D2_nino)
Div_hill_nino_df<-base::as.data.frame(Div_hill_nino)
base::colnames(Div_hill_nino_df)<-c("q0", "q1","q2")

Div_hill_nino_df_trans<-as.data.frame(t(Div_hill_nino_df))
Div_hill_nino_df_trans$Orden<-c(0,1,2)



#Neutro

D0_neutro<-hilldiv::hill_div(densidad_neutro_ord, qvalue = 0)
D1_neutro<-hilldiv::hill_div(densidad_neutro_ord, qvalue = 1)
D2_neutro<-hilldiv::hill_div(densidad_neutro_ord, qvalue = 2)



Div_hill_neutro<-base::list(D0_neutro,D1_neutro,D2_neutro)
Div_hill_neutro_df<-base::as.data.frame(Div_hill_neutro)
base::colnames(Div_hill_neutro_df)<-c("q0", "q1","q2")

Div_hill_neutro_df_trans<-as.data.frame(t(Div_hill_neutro_df))
Div_hill_neutro_df_trans$Orden<-c(0,1,2)
ggplot()+
  geom_point(aes(y=densidad_nino_ord[["Niño"]], x=seq(1:208)), color ="#d7191c")+
  
  geom_point(aes(y=densidad_neutro_ord[["Neutro"]], x=seq(1:215)), color ="#808080")+
  scale_y_continuous(trans='log10', limits = function(x){c(min(x), ceiling(100000))})+
  labs(colour = "", subtitle = "Rango - Abundancia", title="Abundancia Total",  x ="Rango", y =expression(paste("log(abundancia) [ind 10 ", m^-2, "]")),tag="A.")+
   geom_label(data = Div_hill_nino_df_trans[nrow(Div_hill_nino_df_trans), ],aes(x = 10, y = 700, label = "Niño"), color = "#d7191c", fill = "white")+
  geom_label(data = Div_hill_nino_df_trans[nrow(Div_hill_nino_df_trans), ], aes(x = 30, y = 5000, label = "Neutro"), color = "#808080", fill = "white")+
 
  theme_bw()+
  theme(legend.position="bottom", 
        axis.text.y =  element_text(size = 8),
        axis.title =  element_text(size = 8),
        plot.title = element_text(size = 8),
        plot.subtitle = element_text(size = 8),
        plot.caption= element_text(size = 8),
        plot.tag = element_text(size=8))

valor = 400000
 
Calculo_Inext_ENOS<-iNEXT(Diversidad_Total_Enos,q=c(0,1,2),datatype = "abundance",endpoint = valor)
data_quantity_size<-Calculo_Inext_ENOS[["iNextEst"]][["size_based"]]


subset_observed = data_quantity_size%>%
              subset(Method == "Observed")

subset_Extrapolation = data_quantity_size%>%
              subset(Method == "Extrapolation")

subset_Rarefaction = data_quantity_size%>%
              subset(Method == "Rarefaction")

subset_extrapolation_max = subset_Extrapolation%>%
    subset(m==valor)
grafica_extrapolation <- ggplot() +
  # Línea sólida para Rarefacción
  geom_line(data = subset_Rarefaction, aes(x = m, y = qD, color = Assemblage, linetype = "Rarefaction")) +  
  # Puntos para valores observados
  geom_point(data = subset_observed, aes(x = m, y = qD, color = Assemblage), size = 2) +  
  # Línea discontinua para Extrapolación
  geom_line(data = subset_Extrapolation, aes(x = m, y = qD, color = Assemblage, linetype = "Extrapolation")) +  
  # Bandas de confianza
  geom_ribbon(data = subset_Rarefaction, aes(x = m, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  
  geom_ribbon(data = subset_Extrapolation, aes(x = m, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  
  # Definir colores de línea
  scale_color_manual(values = c("Neutro" = "#808080", "Niño" = "#d7191c")) +  
  # Definir colores de bandas de confianza
  scale_fill_manual(values = c("Neutro" = "#808080", "Niño" = "#d7191c")) +  
  # Definir tipos de línea según método
  scale_linetype_manual(values = c("Rarefaction" = "solid", "Extrapolation" = "dashed")) +  
  # Separar por orden de diversidad
  facet_wrap(~ Order.q) +  
  # Estilo general
  theme_bw(base_size = 6) +  
  labs(
    x = "Número de Individuos",
    y = "Diversidad de especies"
  ) +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 9),
        axis.text = element_text(size = 9),
        axis.title = element_text(size = 9),
        strip.text = element_text(size = 9),
        plot.title = element_text(size = 9),
        plot.subtitle = element_text(size = 9),
        plot.caption = element_text(size = 9))

# Mostrar la gráfica
print(grafica_extrapolation)

data_coverage<-Calculo_Inext_ENOS[["iNextEst"]][["coverage_based"]]


subset_observed_coverage = data_coverage%>%
              subset(Method == "Observed")

subset_Extrapolation_coverage = data_coverage%>%
              subset(Method == "Extrapolation")

subset_Rarefaction_coverage = data_coverage%>%
              subset(Method == "Rarefaction")

subset_extrapolation_max_coverage = subset_Extrapolation_coverage%>%
    subset(m==valor)
grafica_extrapolation <- ggplot() +
  # Línea sólida para Rarefacción
  geom_line(data = subset_Rarefaction_coverage, aes(x = SC, y = qD, color = Assemblage, linetype = "Rarefaction")) +  
  # Puntos para valores observados
  geom_point(data = subset_observed_coverage, aes(x = SC, y = qD, color = Assemblage), size = 2) +  
  # Línea discontinua para Extrapolación
  geom_line(data = subset_Extrapolation_coverage, aes(x = SC, y = qD, color = Assemblage, linetype = "Extrapolation")) +  
  # Bandas de confianza
  geom_ribbon(data = subset_Rarefaction_coverage, aes(x = SC, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  
  geom_ribbon(data = subset_Extrapolation_coverage, aes(x = SC, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  
  # Definir colores de línea
  scale_color_manual(values = c("Neutro" = "#808080", "Niño" = "#d7191c")) +  
  # Definir colores de bandas de confianza
  scale_fill_manual(values = c("Neutro" = "#808080", "Niño" = "#d7191c")) +  
  # Definir tipos de línea según método
  scale_linetype_manual(values = c("Rarefaction" = "solid", "Extrapolation" = "dashed")) +  
  # Separar por orden de diversidad
  facet_wrap(~ Order.q) +  
  # Estilo general
  theme_bw(base_size = 6) +  
  labs(
    x = "Cobertura del muestreo",
    y = "Diversidad de especies"
  ) +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 9),
        axis.text = element_text(size = 9),
        axis.title = element_text(size = 9),
        strip.text = element_text(size = 9),
        plot.title = element_text(size = 9),
        plot.subtitle = element_text(size = 9),
        plot.caption = element_text(size = 9))

# Mostrar la gráfica
print(grafica_extrapolation)

ggplot2::ggplot() +
  geom_line(data = subset_observed, aes(x = Order.q, y = qD , color = Assemblage ),  lwd = 0.5, linetype = 1) +
  geom_point(data = subset_observed, aes(x = Order.q, y = qD , color = Assemblage ), size =3) +
    
  
  scale_color_manual(values = c("Neutro" = "#808080", "Niño" = "#d7191c")) +  
  geom_label(data = Div_hill_nino_df_trans[nrow(Div_hill_nino_df_trans), ], 
           aes(x = 0.1, y = 180, label = "Niño"), color = "#d7191c", fill = "white")+
  geom_label(data = Div_hill_nino_df_trans[nrow(Div_hill_nino_df_trans), ], 
           aes(x = 0.15, y = 217, label = "Neutro"), color = "#808080", fill = "white")+
  

  labs(
    title = "Perfiles de diversidad observada",
    x = "Orden",
    y = "Diversidad",
    color = "Variable") +
  theme_minimal()+
   theme(legend.position = "none")  # Elimina la leyenda

3 Análisis por Cruceros

datos_anos<-data_quantity[,c(8,10:286)]


datos_anos_trans<-as.data.frame(t(datos_anos %>%
  group_by(CODE) %>%
  summarise(across(everything(), sum, na.rm = TRUE))))

colnames(datos_anos_trans) <- datos_anos_trans[1, ]  # Usa la primera fila como nombres de columna
datos_anos_trans <- datos_anos_trans[-1, ]  # Elimina la primera fila repetida
datos_anos_trans <- datos_anos_trans %>%
  mutate(across(everything(), as.integer))

datos_anos_trans<-datos_anos_trans%>%select(  "1991_Febrero",   
                                              "1991_Septiembre",
                                              "1991_Diciembre", 
                                              "1993_Enero" ,     
                                              "1993_Noviembre",  
                                              "1994_Abril",  
                                              "1994_Julio",      
                                              "1994_Diciembre" ,     
                                              "1995_Junio", 
                                              "1996_Mayo",       
                                              "1996_Noviembre" , 
                                              "2009_Diciembre" )
reactable(datos_anos_trans)
feb_1991<-datos_anos%>%
        filter(CODE == "1991_Febrero")

sep_1991<-datos_anos%>%
        filter(CODE == "1991_Septiembre")
dic_1991<-datos_anos%>%
        filter(CODE == "1991_Diciembre")
ene_1993<-datos_anos%>%
        filter(CODE == "1993_Enero")
nov_1993<-datos_anos%>%
        filter(CODE == "1993_Noviembre")
abr_1994<-datos_anos%>%
        filter(CODE == "1994_Abril")
jul_1994<-datos_anos%>%
        filter(CODE == "1994_Julio")
dic_1994<-datos_anos%>%
        filter(CODE == "1994_Diciembre")
jun_1995<-datos_anos%>%
        filter(CODE == "1995_Junio")

may_1996<-datos_anos%>%
        filter(CODE == "1996_Mayo")
nov_1996<-datos_anos%>%
        filter(CODE == "1996_Noviembre")

dic_2009<-datos_anos%>%
        filter(CODE == "2009_Diciembre")



feb_1991_total<-as.data.frame(colSums(feb_1991[,2:278]))
sep_1991_total<-as.data.frame(colSums(sep_1991[,2:278]))
dic_1991_total<-as.data.frame(colSums(dic_1991[,2:278]))
ene_1993_total<-as.data.frame(colSums(ene_1993[,2:278]))
nov_1993_total<-as.data.frame(colSums(nov_1993[,2:278]))
abr_1994_total<-as.data.frame(colSums(abr_1994[,2:278]))
jul_1994_total<-as.data.frame(colSums(jul_1994[,2:278]))
dic_1994_total<-as.data.frame(colSums(dic_1994[,2:278]))
jun_1995_total<-as.data.frame(colSums(jun_1995[,2:278]))
may_1996_total<-as.data.frame(colSums(may_1996[,2:278]))
nov_1996_total<-as.data.frame(colSums(nov_1996[,2:278]))
dic_2009_total<-as.data.frame(colSums(dic_2009[,2:278]))



colnames(feb_1991_total)<-c("feb_1991")

colnames(sep_1991_total)<-c("sep_1991")
colnames(dic_1991_total)<-c("dic_1991")
colnames(ene_1993_total)<-c("ene_1993")
colnames(nov_1993_total)<-c("nov_1993")
colnames(abr_1994_total)<-c("abr_1994")
colnames(jul_1994_total)<-c("jul_1994")
colnames(dic_1994_total)<-c("dic_1994")
colnames(jun_1995_total)<-c("jun_1995")

colnames(may_1996_total)<-c("may_1996")
colnames(nov_1996_total)<-c("nov_1996")

colnames(dic_2009_total)<-c("dic_2009")


# Ordenar por la columna A de manera descendente
feb_1991_ord <- feb_1991_total %>%  arrange(desc(feb_1991))%>%  filter(feb_1991 != 0)

sep_1991_ord <- sep_1991_total %>%  arrange(desc(sep_1991))%>%  filter(sep_1991 != 0)
dic_1991_ord <- dic_1991_total %>%  arrange(desc(dic_1991))%>%  filter(dic_1991 != 0)
ene_1993_ord <- ene_1993_total %>%  arrange(desc(ene_1993))%>%  filter(ene_1993 != 0)
nov_1993_ord <- nov_1993_total %>%  arrange(desc(nov_1993))%>%  filter(nov_1993 != 0)
abr_1994_ord <- abr_1994_total %>%  arrange(desc(abr_1994))%>%  filter(abr_1994 != 0)
jul_1994_ord <- jul_1994_total %>%  arrange(desc(jul_1994))%>%  filter(jul_1994 != 0)
dic_1994_ord <- dic_1994_total %>%  arrange(desc(dic_1994))%>%  filter(dic_1994 != 0)
jun_1995_ord <- jun_1995_total %>%  arrange(desc(jun_1995))%>%  filter(jun_1995 != 0)

may_1996_ord <- may_1996_total %>%  arrange(desc(may_1996))%>%  filter(may_1996 != 0)
nov_1996_ord <- nov_1996_total %>%  arrange(desc(nov_1996))%>%  filter(nov_1996 != 0)

dic_2009_ord <- dic_2009_total %>%  arrange(desc(dic_2009))%>%  filter(dic_2009 != 0)



colnames(feb_1991_ord)<-c("feb_1991")

colnames(sep_1991_ord)<-c("sep_1991")
colnames(dic_1991_ord)<-c("dic_1991")
colnames(ene_1993_ord)<-c("ene_1993")
colnames(nov_1993_ord)<-c("nov_1993")
colnames(abr_1994_ord)<-c("abr_1994")
colnames(jul_1994_ord)<-c("jul_1994")
colnames(dic_1994_ord)<-c("dic_1994")
colnames(jun_1995_ord)<-c("jun_1995")

colnames(may_1996_ord)<-c("may_1996")
colnames(nov_1996_ord)<-c("nov_1996")

colnames(dic_2009_ord)<-c("dic_2009")
ggplot() +
  geom_point(aes(y = feb_1991_ord[["feb_1991"]], x = seq_len(152), color = "feb_1991")) +
 
  geom_point(aes(y = sep_1991_ord[["sep_1991"]], x = seq_len(122), color = "sep_1991")) +
  geom_point(aes(y = dic_1991_ord[["dic_1991"]], x = seq_len(115), color = "dic_1991")) +
  geom_point(aes(y = ene_1993_ord[["ene_1993"]], x = seq_len(62), color = "ene_1993")) +
  geom_point(aes(y = nov_1993_ord[["nov_1993"]], x = seq_len(46), color = "nov_1993")) +
  geom_point(aes(y = abr_1994_ord[["abr_1994"]], x = seq_len(13), color = "abr_1994")) +
  geom_point(aes(y = jul_1994_ord[["jul_1994"]], x = seq_len(50), color = "jul_1994")) +
  geom_point(aes(y = dic_1994_ord[["dic_1994"]], x = seq_len(29), color = "dic_1994")) +
  geom_point(aes(y = jun_1995_ord[["jun_1995"]], x = seq_len(120), color = "jun_1995")) +
 
  geom_point(aes(y = may_1996_ord[["may_1996"]], x = seq_len(100), color = "may_1996")) +
  geom_point(aes(y = nov_1996_ord[["nov_1996"]], x = seq_len(104), color = "nov_1996")) +
  
  geom_point(aes(y = dic_2009_ord[["dic_2009"]], x = seq_len(136), color = "dic_2009")) +

  scale_y_continuous(trans = 'log10', limits = c(NA, 100000)) +  # <- Asegúrate de cerrar correctamente

  labs(
    colour = "", 
    subtitle = "Rango - Abundancia", 
    title = "Abundancia Total",  
    x = "Rango", 
    y = expression(paste("log(abundancia) [ind 10 ", m^-2, "]")),
    tag = "A."
  ) +

  scale_color_manual(values = c(
    "feb_1991"  = "#e41a1c",
    
    "sep_1991"  = "#4daf4a",
    "dic_1991"  = "#984ea3",
    "ene_1993"  = "#ff7f00",
    "nov_1993"  = "#ffff33",
    "abr_1994"  = "#a65628",
    "jul_1994"  = "#f781bf",
    "dic_1994"  = "#999999",
    "jun_1995"  = "#66c2a5",
   
    "may_1996"  = "#8da0cb",
    "nov_1996"  = "#e78ac3",
    
    "dic_2009"  = "#ffd92f"
  )) +

  theme_bw() +
  theme(
    legend.position = "bottom", 
    axis.text.y = element_text(size = 8),
    axis.title = element_text(size = 8),
    plot.title = element_text(size = 8),
    plot.subtitle = element_text(size = 8),
    plot.caption = element_text(size = 8),
    plot.tag = element_text(size = 8)
  )

valor = 500000
 
Calculo_Inext_cruceros<-iNEXT(datos_anos_trans,q=c(0,1,2),datatype = "abundance",endpoint = valor)
rarefaccion_size_cruceros<-Calculo_Inext_cruceros[["iNextEst"]][["size_based"]]


cruceros_subset_observed = rarefaccion_size_cruceros%>%
              subset(Method == "Observed")

cruceros_subset_Extrapolation = rarefaccion_size_cruceros%>%
              subset(Method == "Extrapolation")

cruceros_subset_Rarefaction = rarefaccion_size_cruceros%>%
              subset(Method == "Rarefaction")

cruceros_subset_extrapolation_max = cruceros_subset_Extrapolation%>%
    subset(m==valor)
ggplot() +
  # Línea sólida para Rarefacción
  geom_line(data = cruceros_subset_Rarefaction, aes(x = m, y = qD, color = Assemblage, linetype = "Rarefaction")) +  
  # Puntos para valores observados
  geom_point(data = cruceros_subset_observed, aes(x = m, y = qD, color = Assemblage), size = 2) +  
  # Línea discontinua para Extrapolación
  geom_line(data = cruceros_subset_Extrapolation, aes(x = m, y = qD, color = Assemblage, linetype = "Extrapolation")) +  
  # Bandas de confianza
  geom_ribbon(data = cruceros_subset_Rarefaction, aes(x = m, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  
  geom_ribbon(data = cruceros_subset_Extrapolation, aes(x = m, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  

  # Escalas de color y relleno corregidas
  scale_color_manual(values = c(
    
    "1991_Febrero"  = "#1f78b4",
    "1991_Septiembre"  = "#b2df8a",
    "1991_Diciembre"  = "#a6cee3",
    "1993_Enero"  = "#33a02c",
    "1993_Noviembre"  = "#fb9a99",
    "1994_Abril"  = "#e31a1c",
    "1994_Julio"  = "#ff7f00",
    "1994_Diciembre"  = "#fdbf6f",
    "1995_Junio"  = "#cab2d6",
    "1996_Mayo"  = "#6a3d9a",
    "1996_Noviembre"  = "#ffff99",
    "2009_Diciembre"  = "#ffd92f"
  )) + 
  scale_fill_manual(values = c(
     "1991_Febrero"  = "#1f78b4",
    "1991_Septiembre"  = "#b2df8a",
    "1991_Diciembre"  = "#a6cee3",
    "1993_Enero"  = "#33a02c",
    "1993_Noviembre"  = "#fb9a99",
    "1994_Abril"  = "#e31a1c",
    "1994_Julio"  = "#ff7f00",
    "1994_Diciembre"  = "#fdbf6f",
    "1995_Junio"  = "#cab2d6",
    "1996_Mayo"  = "#6a3d9a",
    "1996_Noviembre"  = "#ffff99",
    "2009_Diciembre"  = "#ffd92f"
  )) + 

  # Tipos de línea
  scale_linetype_manual(values = c("Rarefaction" = "solid", "Extrapolation" = "dashed")) +  

  # Facetas por orden de diversidad
  facet_wrap(~ Order.q) +  

  # Estilo general
  theme_bw(base_size = 6) +  
  labs(
    x = "Número de Individuos",
    y = "Diversidad de especies"
  ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.text = element_text(size = 9),
    axis.text = element_text(size = 9),
    axis.title = element_text(size = 9),
    strip.text = element_text(size = 9),
    plot.title = element_text(size = 9),
    plot.subtitle = element_text(size = 9),
    plot.caption = element_text(size = 9)
  )

cruceros_data_coverage<-Calculo_Inext_cruceros[["iNextEst"]][["coverage_based"]]


cruceros_subset_observed_coverage = cruceros_data_coverage%>%
              subset(Method == "Observed")

cruceros_subset_Extrapolation_coverage = cruceros_data_coverage%>%
              subset(Method == "Extrapolation")

cruceros_subset_Rarefaction_coverage = cruceros_data_coverage%>%
              subset(Method == "Rarefaction")

cruceros_subset_extrapolation_max_coverage = subset_Extrapolation_coverage%>%
    subset(m==valor)
ggplot() +
  # Línea sólida para Rarefacción
  geom_line(data = cruceros_subset_Rarefaction_coverage, aes(x = SC, y = qD, color = Assemblage, linetype = "Rarefaction")) +  
  # Puntos para valores observados
  geom_point(data = cruceros_subset_observed_coverage, aes(x = SC, y = qD, color = Assemblage), size = 2) +  
  # Línea discontinua para Extrapolación
  geom_line(data = cruceros_subset_Extrapolation_coverage, aes(x = SC, y = qD, color = Assemblage, linetype = "Extrapolation")) +  
  # Bandas de confianza
  geom_ribbon(data = cruceros_subset_Rarefaction_coverage, aes(x = SC, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  
  geom_ribbon(data = cruceros_subset_Extrapolation_coverage, aes(x = SC, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  

   # Escalas de color y relleno corregidas
  scale_color_manual(values = c(
     "1991_Febrero"  = "#1f78b4",
    "1991_Septiembre"  = "#b2df8a",
    "1991_Diciembre"  = "#a6cee3",
    "1993_Enero"  = "#33a02c",
    "1993_Noviembre"  = "#fb9a99",
    "1994_Abril"  = "#e31a1c",
    "1994_Julio"  = "#ff7f00",
    "1994_Diciembre"  = "#fdbf6f",
    "1995_Junio"  = "#cab2d6",
    "1996_Mayo"  = "#6a3d9a",
    "1996_Noviembre"  = "#ffff99",
    "2009_Diciembre"  = "#ffd92f"
  )) + 
  scale_fill_manual(values = c(
      "1991_Febrero"  = "#1f78b4",
    "1991_Septiembre"  = "#b2df8a",
    "1991_Diciembre"  = "#a6cee3",
    "1993_Enero"  = "#33a02c",
    "1993_Noviembre"  = "#fb9a99",
    "1994_Abril"  = "#e31a1c",
    "1994_Julio"  = "#ff7f00",
    "1994_Diciembre"  = "#fdbf6f",
    "1995_Junio"  = "#cab2d6",
    "1996_Mayo"  = "#6a3d9a",
    "1996_Noviembre"  = "#ffff99",
    "2009_Diciembre"  = "#ffd92f"
  )) + 

  # Tipos de línea
  scale_linetype_manual(values = c("Rarefaction" = "solid", "Extrapolation" = "dashed")) +  

  # Facetas por orden de diversidad
  facet_wrap(~ Order.q) +  

  # Estilo general
  theme_bw(base_size = 6) +  
  labs(
    x = "Cobertura del muestreo",
    y = "Diversidad de especies"
  ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.text = element_text(size = 9),
    axis.text = element_text(size = 9),
    axis.title = element_text(size = 9),
    strip.text = element_text(size = 9),
    plot.title = element_text(size = 9),
    plot.subtitle = element_text(size = 9),
    plot.caption = element_text(size = 9)
  )

ggplot2::ggplot() +
  geom_line(data = cruceros_subset_observed, aes(x = Order.q, y = qD , color = Assemblage ),  lwd = 0.5, linetype = 1) +
  geom_point(data = cruceros_subset_observed, aes(x = Order.q, y = qD , color = Assemblage ), size =3) +
  
   scale_color_manual(values = c(
      "1991_Febrero"  = "#1f78b4",
    "1991_Septiembre"  = "#b2df8a",
    "1991_Diciembre"  = "#a6cee3",
    "1993_Enero"  = "#33a02c",
    "1993_Noviembre"  = "#fb9a99",
    "1994_Abril"  = "#e31a1c",
    "1994_Julio"  = "#ff7f00",
    "1994_Diciembre"  = "#fdbf6f",
    "1995_Junio"  = "#cab2d6",
    "1996_Mayo"  = "#6a3d9a",
    "1996_Noviembre"  = "#ffff99",
    "2009_Diciembre"  = "#ffd92f")) + 


  labs(
    title = "Perfiles de diversidad observada",
    x = "Orden",
    y = "Diversidad",
    color = "Crucero") +
  theme_minimal()+
   theme(legend.position = "bottom")  # Elimina la leyenda

ggplot2::ggplot() +
  geom_line(data = cruceros_subset_extrapolation_max, aes(x = Order.q, y = qD , color = Assemblage ),  lwd = 0.5, linetype = 1) +
  geom_point(data = cruceros_subset_extrapolation_max, aes(x = Order.q, y = qD , color = Assemblage ), size =3) +
  geom_ribbon(data = cruceros_subset_extrapolation_max, aes(x = Order.q, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +  
  geom_ribbon(data = cruceros_subset_extrapolation_max, aes(x = Order.q, ymin = qD.LCL, ymax = qD.UCL, fill = Assemblage), alpha = 0.2) +
  
   scale_color_manual(values = c(  "1991_Febrero"  = "#1f78b4",
    "1991_Septiembre"  = "#b2df8a",
    "1991_Diciembre"  = "#a6cee3",
    "1993_Enero"  = "#33a02c",
    "1993_Noviembre"  = "#fb9a99",
    "1994_Abril"  = "#e31a1c",
    "1994_Julio"  = "#ff7f00",
    "1994_Diciembre"  = "#fdbf6f",
    "1995_Junio"  = "#cab2d6",
    "1996_Mayo"  = "#6a3d9a",
    "1996_Noviembre"  = "#ffff99",
    "2009_Diciembre"  = "#ffd92f")) + 
   scale_fill_manual(values = c(  "1991_Febrero"  = "#1f78b4",
    "1991_Septiembre"  = "#b2df8a",
    "1991_Diciembre"  = "#a6cee3",
    "1993_Enero"  = "#33a02c",
    "1993_Noviembre"  = "#fb9a99",
    "1994_Abril"  = "#e31a1c",
    "1994_Julio"  = "#ff7f00",
    "1994_Diciembre"  = "#fdbf6f",
    "1995_Junio"  = "#cab2d6",
    "1996_Mayo"  = "#6a3d9a",
    "1996_Noviembre"  = "#ffff99",
    "2009_Diciembre"  = "#ffd92f")) + 


  labs(
    title = "Perfiles de diversidad estimada",
    x = "Orden",
    y = "Diversidad",
    color = "Crucero",
    fill = "Intervalo de Confianza") +
  theme_minimal()+
   theme(legend.position = "bottom")  # Elimina la leyenda

library(ggplot2)
library(dplyr)
library(tidyr)




  
# Crear dataframe de fechas y etiquetas personalizadas
cruceros <- data.frame(
  date = as.Date(c(
    
"1991-02-01",
"1991-09-01",
"1991-12-01",
"1993-01-01",
"1993-11-01",
"1994-04-01",
"1994-07-01",
"1994-12-01",
"1995-05-01",
"1996-05-01",
"1996-11-01",
"2009-12-01")))

# Filtrar datos por orden de diversidad
q0 <- cruceros_subset_observed %>% subset(Order.q == 0)
q1 <- cruceros_subset_observed %>% subset(Order.q == 1)
q2 <- cruceros_subset_observed %>% subset(Order.q == 2)

# Agregar datos de diversidad a 'cruceros'
cruceros$q0 <- q0$qD
cruceros$q1 <- q1$qD
cruceros$q2 <- q2$qD

# Convertir a formato largo para ggplot2
cruceros_long <- pivot_longer(cruceros, cols = c(q0, q1, q2), 
                               names_to = "Diversidad", values_to = "Valor")

# Asignar colores a cada índice de diversidad
colores_diversidad <- c("q0" = "#1b9e77", "q1" = "#d95f02", "q2" = "#7570b3")

# Gráfica con etiquetas personalizadas en el eje x

# Gráfica con leyenda para las líneas
ggplot(cruceros_long, aes(x = date, y = Valor, color = Diversidad)) + 
  geom_line(size = 1) +  
  geom_point(size = 2) +

  scale_x_date(limits = c(as.Date("1991-01-01"), as.Date("2010-12-31")), 
               date_labels = "%Y", date_breaks = "1 years") +
  
  scale_color_manual(values = colores_diversidad, 
                     labels = c("q0" = "Diversidad q0", 
                                "q1" = "Diversidad q1", 
                                "q2" = "Diversidad q2")) +
  
  labs(
    x = "Fecha", 
    y = "Índice de Diversidad",
    color = "Número de especies",  # Título de la leyenda
    title = "Comportamiento de la diversidad en cada crucero"
  ) +
  
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

3.0.1 Mapas de diversidad

data_larvas<-readxl::read_excel("../data/raw/biodiversity/PELAGDEMER_Especies.xlsx", sheet="Larvas")


breaks_LARVAS <- c(1, 10, 100, 1000, 10000)
data_larvas$LARVAS_log <- cut(data_larvas$Larvas,
                       breaks = c(-Inf, breaks_LARVAS, Inf), # Incluye -Inf y Inf para abarcar todos los valores posibles
                       labels = c("1", "10", "100", "1000", "10000", "100000"),
                       right = FALSE)  # right = FALSE significa que los intervalos son cerrados por la izquierda, abiertos por la derecha

data_larvas$LARVAS_log<-as.character(data_larvas$LARVAS_log )
data_larvas$LARVAS_log<-as.integer(data_larvas$LARVAS_log )
data_larvas$LARVAS_log<-as.factor(data_larvas$LARVAS_log )
class(data_larvas$LARVAS_log)
## [1] "factor"
data_larvas<-readxl::read_excel("../data/raw/biodiversity/PELAGDEMER_Especies.xlsx", sheet="Larvas")


breaks_LARVAS <- c(1, 10, 100, 1000, 10000)
data_larvas$LARVAS_log <- cut(data_larvas$Larvas,
                       breaks = c(-Inf, breaks_LARVAS, Inf), # Incluye -Inf y Inf para abarcar todos los valores posibles
                       labels = c("1", "10", "100", "1000", "10000", "100000"),
                       right = FALSE)  # right = FALSE significa que los intervalos son cerrados por la izquierda, abiertos por la derecha

data_larvas$LARVAS_log<-as.character(data_larvas$LARVAS_log )
data_larvas$LARVAS_log<-as.integer(data_larvas$LARVAS_log )
data_larvas$LARVAS_log<-as.factor(data_larvas$LARVAS_log )
class(data_larvas$LARVAS_log)
## [1] "factor"
#Carga de las capas para el mapa

CPC<-sf::st_read("../data/gis/GeoLayers.gpkg", layer="cuenca_pacifica") # read shapefile CPC
## Reading layer `cuenca_pacifica' from data source 
##   `/home/chrisbermudezr/IchthyoplanktonERFEN/data/gis/GeoLayers.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1 feature and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -86.00392 ymin: 1.439583 xmax: -77.01431 ymax: 7.218436
## Geodetic CRS:  WGS 84
Paises<-sf::st_read("../data/gis/GeoLayers.gpkg", layer="Continente") # read shapefile Countries
## Reading layer `Continente' from data source 
##   `/home/chrisbermudezr/IchthyoplanktonERFEN/data/gis/GeoLayers.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 23 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -94.33045 ymin: -6.996088 xmax: -64.48952 ymax: 18.93407
## Geodetic CRS:  WGS 84
MPA<-sf::st_read("../data/gis/GeoLayers.gpkg", layer="mpa2023_cpc") # read shapefile Countries
## Reading layer `mpa2023_cpc' from data source 
##   `/home/chrisbermudezr/IchthyoplanktonERFEN/data/gis/GeoLayers.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 13 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.76667 ymin: 1.43053 xmax: -77.09094 ymax: 6.152708
## Geodetic CRS:  WGS 84
####Larvas####
if(!require(readxl)) install.packages("readxl")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(dplyr)) install.packages("dplyr")

Mapas_graficas_LARVAS<-function(data, variable, titulo, subtitulo, leyenda) {
  ggplot() +
    geom_sf(data = CPC, color = "blue", linetype = 2, linewidth = 0.5, fill = "lightblue") +
    geom_sf(data = Paises, colour = "black", fill = "lightgrey") +
    geom_sf(data = MPA, color = "darkgreen", linetype = 1, linewidth = 0.5, , fill = "lightblue") +
    geom_point(data = data, aes(x = LONGITUD, y = LATITUD,  size = variable), colour="red") +
    coord_sf(xlim = c(-80, -77), ylim = c(1, 8), expand = FALSE) +
    
  scale_x_continuous(breaks = seq(-80, -77, by = 1)) +
    scale_y_continuous(breaks = seq(1, 8, by = 1)) +

    labs(
      title = titulo,
      subtitle = subtitulo,
      x = "Longitude",
      y = "Latitude",
      size = leyenda
    ) +
    scale_size_manual(values = c( `10` = 1, `100` = 2, `1000` = 3, `10000` = 4, `100000` = 6)) +
    # Ajusta los tamaños según las categorías
    #scale_size_continuous(range = c(1, 5),breaks = c(10, 100, 1000, 10000), labels = scales::label_number(accuracy = 1)) +
    theme_bw() +
    theme(
      plot.title = element_text(size = 12, face = "italic", color = "black"),
      axis.title = element_text(face = "bold", color = "black")
    )
}



options(scipen=9999)
Pelagicos_1991_02<-subset(data_larvas,data_larvas$CODE == "1991_Febrero")

Pelagicos_1991_02<-Pelagicos_1991_02%>%dplyr::filter( Larvas > 0)
Pelagicos_1991_02_LARVAS<- Mapas_graficas_LARVAS(Pelagicos_1991_02, 
                                             Pelagicos_1991_02$LARVAS_log, 
                                             expression(paste(bold("1991-02"))),
                                             expression(paste(bold(" "))), 
                                              expression(paste("[larvas.10m"^{-2}, "]")))

png(filename="Pelagicos_1991_02_LARVAS_log_10.png", height = 15, width =  10, units = "cm", res = 300, pointsize = 12)
Pelagicos_1991_02_LARVAS
dev.off()
## png 
##   2
Pelagicos_1991_02_LARVAS