rm(list=ls())

# install.packages("ggrepel")
# install.packages("gridExtra")
library(ggplot2)
library(ggrepel)
library(gridExtra)

setwd('D:\\areas_protegidas\\xls_dbf\\salidas txt\\')
paises <- c("All countries","Argentina","Bolivia","Brazil","Chile","Colombia",
            "Ecuador","Guayana Francesa","Guayanas","Guyana","Peru",
            "Paraguay","Suriname","Uruguay","Venezuela")

gini_AP<- read.table("td_AP_pais.txt")
# gini_AP <- data.frame(pais=rownames(gini_AP), gini_AP)
colnames(gini_AP)<-"AP_biof"
gini_AP$pais<-paises
gini_AP <- gini_AP[order(gini_AP$pais),]
print(gini_AP)
class(gini_AP)

gini_FP<-read.table('td_FP_pais.txt')
# gini_FP <- data.frame(pais=rownames(gini_FP), gini_FP)
colnames(gini_FP)<-"FP_biof"
gini_FP$pais<-paises
gini_FP <- gini_FP[order(gini_FP$pais),]
print(gini_FP)
class(gini_FP)

gini_eco_pais<-read.table("gini_eco_pais.txt",sep=';',header=T)  
gini_eco_pais <- data.frame(pais=rownames(gini_eco_pais), gini_eco_pais)
gini_eco_pais$giniAP_1a6<-NULL
gini_eco_pais$giniFP_1a6<-NULL
gini_eco_pais <- gini_eco_pais[-(7),] # saco FK
gini_eco_pais <- gini_eco_pais[order(gini_eco_pais$pais),] # cambio de lugar GI
colnames(gini_eco_pais)<-c("pais","AP_eco","FP_eco")
print(gini_eco_pais)
class(gini_eco_pais)

gini_union<-cbind(gini_AP, gini_FP,gini_eco_pais)
gini_union<-gini_union[-c(1,8,10,13),] #-c(7,9,12) cuando no esta SAM
print(gini_union)
names(gini_union)
class(gini_union)

# Invierto el Gini para que sea una medidad de representatividad
gini_union$AP_biof_inv<- 1-gini_union$AP_biof
gini_union$FP_biof_inv<- 1-gini_union$FP_biof
gini_union$AP_eco_inv<- 1-gini_union$AP_eco
gini_union$FP_eco_inv<- 1-gini_union$FP_eco

# para generar un campo para la ecuacion sin UY
gini_union$FP_biof_inv_sUY<- 1-gini_union$FP_biof
gini_union$FP_biof_inv_sUY[gini_union$pais=="Uruguay"] <- NA

gini_union$AP_eco_inv[gini_union$pais=="Uruguay"] <- NA #porque tiene una sola ecoregion
gini_union$FP_eco_inv[gini_union$pais=="Uruguay"] <- NA

##############################################################
# GRAFICOS XY

# First plot
p1 <- ggplot(data=gini_union, aes(x=AP_biof_inv, y=FP_biof_inv)) +
      geom_point(size=1.5) +
      geom_text_repel(aes(label=pais), size = 2.5)+
      geom_point(aes(x=0.425074237,y=0.773586204), size=1.5, colour = "royalblue4") +
      annotate("text", x=0.465,y=0.773586204, label = "All countries", size = 2.5) +
      theme_grey(base_size = 8, base_family = "")+
      labs(x = expression(QR~ biophysical~(1-G)), y=expression(FR~biophysical~gradients~(1-G)), size = 2.5)+
      geom_abline(intercept = 0, slope = 1, colour = "red4", linetype = "dashed")
p1
ggsave("xy_representatividades_p1_1-G.png",plot = p1, scale = 1, width = 7, height = 7, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')

# Second plot
p2 <- ggplot(data=gini_union, aes(x=AP_eco_inv, y=FP_eco_inv)) +
  geom_point(size=1.5) +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = expression(QR~biogeographical~units~(1-G)), y= expression(FR~biogeographical~units~(1-G)), size = 2.5)+
  geom_abline(intercept = 0, slope = 1, colour = "red4", linetype = "dashed")
p2

# Third plot
p3 <- ggplot(data=gini_union, aes(x=FP_eco_inv, y=FP_biof_inv)) +
  geom_point(size=1.5) +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x =expression(FR~biogeographical~units~(1-G)), y= expression(FR~biophysical~gradients~(1-G)), size = 2.5)+
  geom_abline(intercept = 0, slope = 1, colour = "red4", linetype = "dashed")
p3

# Fourth plot
p4 <- ggplot(data=gini_union, aes(x=AP_eco_inv, y=AP_biof_inv)) +
  geom_point(size=1.5) +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x =expression(QR~biogeographical~units~(1-G)), y= expression(QR~ biophysical~(1-G)), size = 2.5)+
  geom_abline(intercept = 0, slope = 1, colour = "red4", linetype = "dashed")
p4

xy<-grid.arrange(p1, p2, p3, p4, ncol=2)
ggsave("xy_representatividades_1-G.png",plot = xy, scale = 1, width = 14, height = 14, units = "cm",
        dpi = 300, path = '..\\salidas graficas\\graficos sam')


#########################################################################################

datos_basicos<- read.table("datos_basicos.txt", sep="\t", header=T)
# datos_basicos<-datos_basicos[-c(7,9,12,15),] cuando no estaba SAM
datos_basicos<-datos_basicos[-c(1,8,10,13,16),]
print(datos_basicos)
class(datos_basicos)

union<-cbind(gini_union,datos_basicos)
names(union)

union$FR_FP_bio<-union$FP_biof_inv/union$FP
union$QR_FP_bio<-union$AP_biof_inv/union$FP
union$FR_FP_eco<-union$FP_eco_inv/union$FP
union$QR_FP_eco<-union$AP_eco_inv/union$FP

# esto no tendría que estar porque lo puedo limitar en el scale_size_area()
union$FR_FP_bio<-ifelse(union$FR_FP_bio<1.5,union$FR_FP_bio,1.5) 
union$QR_FP_bio<-ifelse(union$QR_FP_bio<1.5,union$QR_FP_bio,1.5)
union$FR_FP_eco<-ifelse(union$FR_FP_eco<1.5,union$FR_FP_eco,1.5)
union$QR_FP_eco<-ifelse(union$QR_FP_eco<1.5,union$FR_FP_bio,1.5)

# para comparar con las clases de Sturges
write.table(union$FP_biof_inv, "10C_FP_biof_inv.txt",col.names=T, row.names=F)

##############################################################
# GRAFICOS RELACIONES ENTRE REPRESENTATIVIDAD Y FRACCION PROTEGIDA

# 1 version, FR y QR bio y eco lm
# First plot
l1 <- ggplot(data=union, aes(x=FP, y=FP_biof_inv)) +
  geom_smooth(linetype=1, size=0.5, colour = "gray50", method = "lm", se=F) +
  geom_point(size=1.5, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Protection extent (%)", y=expression(FR~biophysical~gradients~(1-G)), size = 2.5)
  # geom_abline(model="lm", colour = "dodgerblue4", linetype = "dashed")
l1

# Second plot
l2 <- ggplot(data=union, aes(x=FP, y=AP_biof_inv)) +
  geom_smooth(linetype=1, size=0.5, colour = "gray50", method = "lm", se=F) +
  geom_point(size=1.5, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Protection extent (%)", y=expression(QR~ biophysical~(1-G)), size = 2.5)
# geom_abline(model="lm", colour = "dodgerblue4", linetype = "dashed")
l2

# Third plot
l3 <- ggplot(data=union, aes(x=FP, y=FP_eco_inv)) +
  geom_smooth(linetype=1, size=0.5, colour = "gray50", method = "lm", se=F) +
  geom_point(size=1.5, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Protection extent (%)", y=expression(FR~biogeographical~units~(1-G)), size = 2.5)
# geom_abline(model="lm", colour = "dodgerblue4", linetype = "dashed")
l3

# Fourth plot
l4 <- ggplot(data=union, aes(x=FP, y=AP_eco_inv)) +
  geom_smooth(linetype=1, size=0.5, colour = "gray50", method = "lm", se=F) +
  geom_point(size=1.5, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Protection extent (%)", y=expression(QR~biogeographical~units~(1-G)), size = 2.5)
# geom_abline(model="lm", colour = "dodgerblue4", linetype = "dashed")
l4

dependencia_FP<-grid.arrange(l1, l2, l3, l4, ncol=2)
ggsave("dependencia_FP_1-G.png",plot = dependencia_FP, scale = 1, width = 14, height = 14, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')


# 2 version, FR y QR bio y eco lm, con la relación FP y representatividad
# First plot
l5 <- ggplot(data=union) +
  geom_smooth(aes(x=FP, y=FP_biof_inv), linetype=1, size=1, method = "lm", se=F) +
  geom_point(aes(x=FP, y=FP_biof_inv, size=union$FR_FP_bio), colour = "red4") +
  scale_size_area(breaks=c(0.2,0.4,0.6,0.8,1), max_size=6) +
  geom_text_repel(aes(x=FP, y=FP_biof_inv, label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  theme(legend.position=c(0.9, 0.18), legend.key.size = unit(.4, "cm"),
        legend.title = element_text(size=7)) +
  labs(x = "Protection extent (%)", y=expression(FR~biophysical~gradients~(1-G)), size = "FR:PF")
l5

ggsave("dependencia_FP_version_ln1_1-G.png", plot = l5, scale = 1, width = 7, height = 7, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')

# Second plot
l6 <- ggplot(data=union) +
  geom_smooth(aes(x=FP, y=AP_biof_inv), linetype=1, size=1, method = "lm", se=F) +
  geom_point(aes(x=FP, y=AP_biof_inv, size=union$QR_FP_bio), colour = "gray25") +
  scale_size_area(breaks=c(0.2,0.4,0.6,0.8,1), max_size=6) +
  geom_text_repel(aes(x=FP, y=AP_biof_inv, label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  theme(legend.position=c(0.9, 0.18), legend.key.size = unit(.4, "cm"),
        legend.title = element_text(size=7)) +
  labs(x = "Protection extent (%)", y=expression(QR~ biophysical~(1-G)), size = "QR:PF")
l6

ggsave("dependencia_FP_version_ln2_1-G.png", plot = l6, scale = 1, width = 7, height = 7, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')

# Third plot
l7 <- ggplot(data=union) +
  geom_smooth(aes(x=FP, y=FP_eco_inv), linetype=1, size=1, method = "lm", se=F) +
  geom_point(aes(x=FP, y=FP_eco_inv, size=union$FR_FP_eco), colour = "red4") +
  scale_size_area(breaks=c(0.2,0.4,0.6,0.8,1), max_size=6) +
  geom_text_repel(aes(x=FP, y=FP_eco_inv, label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  theme(legend.position=c(0.9, 0.18), legend.key.size = unit(.4, "cm"),
        legend.title = element_text(size=7)) +
  labs(x = "Protection extent (%)", y=expression(FR~biogeographical~units~(1-G)), size = "FR:PF")
l7

ggsave("dependencia_FP_version_ln3_1-G.png", plot = l7, scale = 1, width = 7, height = 7, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')

# Fourth plot
l8 <- ggplot(data=union) +
  geom_smooth(aes(x=FP, y=AP_eco_inv), linetype=1, size=1, method = "lm", se=F) +
  geom_point(aes(x=FP, y=AP_eco_inv, size=union$QR_FP_eco), colour = "gray25") +
  scale_size_area(breaks=c(0.2,0.4,0.6,0.8,1), max_size=6) +
  geom_text_repel(aes(x=FP, y=AP_eco_inv, label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  theme(legend.position=c(0.9, 0.18), legend.key.size = unit(.4, "cm"),
        legend.title = element_text(size=7)) +
  labs(x = "Protection extent (%)", y=expression(QR~biogeographical~units~(1-G)), size = "QR:PF")
l8

ggsave("dependencia_FP_version_ln4_1-G.png", plot = l8, scale = 1, width = 7, height = 7, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')

dependencia_FP_ln<-grid.arrange(l5, l6, l7, l8, ncol=2)
ggsave("dependencia_FP_version_ln_1-G.png", plot = dependencia_FP_ln, scale = 1, width = 14, height = 14, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')


# 3 version, FR bio y eco lm

# GET EQUATION AND R-SQUARED AS STRING
# SOURCE: http://goo.gl/K4yh

# First plot
lm_eqn_l9 <- function(union){
  m <- lm(FP_biof_inv_sUY ~ FP, union);
  eq <- substitute(italic(paste("G", "'", sep="")) == a + b %.% italic(PEx)*","~~italic(r)^2~"="~r2, 
                   list(a = format(coef(m)[1], digits = 2), 
                        b = format(coef(m)[2], digits = 1), 
                        r2 = format(summary(m)$r.squared, digits = 2)))
  as.character(as.expression(eq));                 
}
l9 <- ggplot(data=union, aes(x=FP, y=FP_biof_inv)) +
  geom_smooth(linetype=1, size=0.5, method = "lm", colour= "gray50", se=F) +
  geom_text(x = 12, y = 0.25, label = lm_eqn_l9(union), parse = TRUE, size = 2.5, colour= "gray50")+  
  geom_point(size=1.5, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5, force = 10)+
  geom_point(aes(x=6.7229968,y=0.773586204), size=1.5, colour = "royalblue4") +
  geom_point(aes(x=0.0879089,y=0.2099528), size=1.5, colour = "royalblue4") +
  annotate("text", x=9.8,y=0.79, label = "All countries", size = 2.5) +
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Protection extent  (%)", y="Protection physical equality  (G')", size = "FR:PF")+
  annotate("text", size=3, x=0.8, y=0.79, label= "(b)")
l9

# Second plot
lm_eqn_l10 <- function(union){
  m <- lm(FP_eco_inv ~ FP, union);
  eq <- substitute(italic(paste("G", "'", sep="")) == a + b %.% italic(PEx)*","~~italic(r)^2~"="~r2, 
                   list(a = format(coef(m)[1], digits = 2), 
                        b = format(coef(m)[2], digits = 1), 
                        r2 = format(summary(m)$r.squared, digits = 1)))
  as.character(as.expression(eq));                 
}
l10 <- ggplot(data=union, aes(x=FP, y=FP_eco_inv)) +
  geom_smooth(linetype=1, size=0.5, method = "lm", colour= "gray50", se=F) +
  geom_text(x = 12, y = 0.25, label = lm_eqn_l10(union), parse = TRUE, size = 2.5, colour= "gray50")+  
  geom_point(size=1.5, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5, force = 22)+
  geom_point(aes(x=6.7229968,y=0.338171301), size=1.5, colour = "royalblue4") +
  annotate("text", x=7.5, y=0.33, label = "All countries", size = 2.5) +
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Protection extent  (%)", y="Protection biogeographical equality (G')", size = "FR:PF")+
  annotate("text", size=3, x=0.8, y=0.47, label= "(a)")
l10

dependencia_FP_ln_FR<-grid.arrange(l10, l9, ncol=2)
ggsave("dependencia_FP_version_ln_FR_1-G.png", plot = dependencia_FP_ln_FR, scale = 1, width = 14, height = 7, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')


##############################################################
# GRAFICOS CONTEXTO ANOS

c1 <- ggplot(data=union, aes(x=anios, y=FP_biof_inv)) +
  scale_x_reverse()+
  scale_y_continuous(position = "right")+
  geom_smooth(linetype=1, size=1, colour = "red4") +
  geom_point(size=1, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Years since first PA", y=expression(FR~biophysical~gradients~(1-G)), size = 2.5)

c2 <- ggplot(data=union, aes(x=anios, y=AP_biof_inv)) +
  scale_x_reverse()+
  scale_y_continuous(position = "right")+
  geom_smooth(linetype=1, size=1, colour = "red4") +
  geom_point(size=1, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Years since first PA", y=expression(QR~biophysical~gradients~(1-G)), size = 2.5)

c3 <- ggplot(data=union, aes(x=anios, y=FP_eco_inv)) +
  scale_x_reverse()+
  scale_y_continuous(position = "right")+
  geom_smooth(linetype=1, size=1, colour = "red4") +
  geom_point(size=1, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Years since first PA", y=expression(FR~biogeographical~units~(1-G)), size = 2.5)

c4 <- ggplot(data=union, aes(x=anios, y=AP_eco_inv)) +
  scale_x_reverse()+
  scale_y_continuous(position = "right")+
  geom_smooth(linetype=1, size=1, colour = "red4") +
  geom_point(size=1, colour = "red4") +
  geom_text_repel(aes(label=pais), size = 2.5)+
  theme_grey(base_size = 8, base_family = "")+
  labs(x = "Years since first PA", y=expression(FR~biogeographical~units~(1-G)), size = 2.5)

dependencia_anios<-grid.arrange(c1, c2, c3, c4, ncol=2)
ggsave("dependencia_anios_1-G.png",plot = dependencia_anios, scale = 1, width = 14, height = 14, units = "cm",
       dpi = 300, path = '..\\salidas graficas\\graficos sam')