rm(list=ls())
library(ineq)
setwd('D:\\areas_protegidas\\xls_dbf')
# setwd("/home/marcos/grive/baldi_ecoaustral")
dat<-read.table('todas_las_variables_sam_simp_fk_1000km_-1990.txt',sep='\t',header=T)

dat$AREA_AL <-  rowSums(dat[,21:33]) # creo el area para todo el continente
dat$AP_AL <- rowSums(dat[,35:47]) # creo el area protegido para todo el continente
dat<-data.frame(dat[,c(1:34,49,35:48,50)])  # muevo las columnas

## POR PAIS, 21 a 34 Argentina a las tres guayanas juntas ------------------------------------------------
# datos_gini_AP <- matrix(nrow = 10,ncol = 14)
# datos_gini_FP <- matrix(nrow = 10,ncol = 14)
datos_gini_AP <- NULL
datos_gini_FP <- NULL
gini_AP <- NULL
gini_FP <- NULL
nom_gini <- NULL

paises <- c("AR","BO","BR","CL","CO","EQ","GF","GU","PE","PY","SU","UY","VZ","GI")
variables<-c("PPT","TMP","ELEV","SLOPE","POP","DIST_ROAD","ACCESS","CULTIVOS",
             "DIST_FRON","DIV_AN","BIOMASS","TEB","FLICKR_POP")

  
for (j in 8:20){
  dat_var <- data.frame(dat[,j], dat[,21:50]);
  for (i in 2:16){
    dat_simp <- data.frame(dat_var[,1], as.numeric(dat_var[,i]), dat_var[,(i+15)], stringsAsFactors = FALSE);
    dat_pais <- subset(dat_simp, dat_simp[,2]>0)
    names(dat_pais) <- c(names(dat)[j], names(dat_var)[c(i,(i+15))])
    print(head(dat_pais))
  
    limINF <- quantile(dat_pais[,1], 0.025)  
    limSUP <- quantile(dat_pais[,1], 0.975)  
  
    dat_pais <- cbind(dat_pais, simp = ifelse(dat_pais[,1] >= limSUP,limSUP,
                  ifelse(dat_pais[,1]<=limINF, limINF, dat_pais[,1])))
 
    cla = seq(limINF,limSUP,length.out = 11) #los histogramas siempre tienen 10 clases
    # cla=seq(floor(limINF), ceiling(limSUP), by=1)
    print(cla)
    names(dat_var)[i]
  
    nom_pais<-substr(names(dat_var)[i],6,7) # aca del campo AREA_PAIS se extrae el nombre del pais
    print(nom_pais)

    nom_variable<-names(dat)[j]
    print(nom_variable)
  
    histograma = hist(dat_pais$simp, breaks=cla, plot=F, include.lowest=T)
    cat = cut(dat_pais$simp, breaks=histograma$breaks, include.lowest=TRUE,labels=F)
    dat_pais  = cbind(dat_pais, cat)
  # tabla dinamica con categorias de las variables independientes 
  # y resumiendo por ellas el area de cada pais y el areas protegida
  # de cada pais
    areasclases <- tapply(as.numeric(dat_pais[,2]), INDEX=dat_pais$cat, FUN=sum)
    areasPclases <- tapply(as.numeric(dat_pais[,3]), INDEX=dat_pais$cat, FUN=sum)
    fraccionPclases <- areasPclases*100/areasclases
    
  # calculo del gini
    gini_areasPclases<- ineq(areasPclases, type = "Gini", na.rm = TRUE)
    gini_fraccionPclases<- ineq(fraccionPclases, type = "Gini", na.rm = TRUE)
    
    tran_nom <-paste(nom_pais, nom_variable, sep="-") 
    nom_gini <- c(nom_gini, tran_nom)
    
  # se van uniendo los datos
    datos_gini_AP <- cbind(datos_gini_AP, areasPclases) 
    datos_gini_FP <- cbind(datos_gini_FP, fraccionPclases)
    gini_AP <- cbind(gini_AP, gini_areasPclases)
    gini_FP <- cbind(gini_FP, gini_fraccionPclases)
    
    # histograma de areas del pais
    postscript(file=paste("salidas graficas/graficos sam/", nom_variable, "_",nom_pais, "_-1990",".ps"), width=600, height=400)
    # pdf(file=paste("salidas graficas/graficos sam/", nom_variable, "_",nom_pais,".pdf"))
    par(mar=c(4,5,1,5), cex.axis=1.8)
    
    namebar <- round(histograma$mids, 1)
    names(namebar)<- 1:length(namebar)    
    areasclases <- areasclases[match(names(namebar), names(areasclases))]
    names(areasclases) <- names(namebar)
    
    barplot(areasclases, axes=F,border="gray70",col="gray80", space=0,
          ylim=c(0, ceiling(max(areasclases, na.rm=T))), 
          names.arg=namebar)
    mtext(nom_variable, side=1, line=3, cex=1.8)

    # histograma de areas protegidas del pais 
    par(new=T)
    areasPclases <- areasPclases[match(names(namebar), names(areasPclases))]
    names(areasPclases) <- names(namebar)
    barplot(areasPclases, axes=F,border="gray70",col="gray30", 
    ylim=c(0, ceiling(max(areasclases, na.rm=T))),space=0, names.arg=F)
    axis(2,col="gray30", col.axis="gray30"); 
    area <- expression(paste("Área (1000 km"^"2",")"))
    mtext(text=area,side=2, line=3, col="gray30", cex=1.8)
    # areap <- expression(paste("Area protegida(km"^"2",")"))
    # mtext(text=areap,side=2, line=2, col="gray30", cex=1.8)
  
    # linea de fraccion de areas protegidas del pais  
    par(new=T)
    fraccionPclases <- fraccionPclases[match(names(namebar), names(fraccionPclases))]
    names(fraccionPclases) <- names(namebar)
    ymax <- ceiling(max(fraccionPclases, na.rm=T))
    plot(fraccionPclases~histograma$mids, xlab="", ylab="", ylim=c(0,ymax), xlim=range(histograma$breaks),pch=16, cex=1.5, type="o", col="red3",  axes=F)  
    axis(4,col="red", col.axis="red3"); 
    mtext(text="Protección (%)",side=4, line=3.2, col="red3", cex=2)
    dev.off()
    }
}
colnames(gini_AP)<- colnames(gini_FP)<- colnames(datos_gini_AP) <- colnames(datos_gini_FP) <- nom_gini
write.table(datos_gini_AP, "salidas txt/datos_gini_AP_-1990.txt",col.names=T, row.names=F)
write.table(datos_gini_FP, "salidas txt/datos_gini_FP_-1990.txt",col.names=T, row.names=F)
write.table(gini_AP, "salidas txt/gini_AP_-1990.txt",col.names=T, row.names=F)
write.table(gini_FP, "salidas txt/gini_FP_-1990.txt",col.names=T, row.names=F)

##############################################################################
# gini_AP

gini_AP_df<-as.data.frame(gini_AP)
gini_AP_df[76:195]<-NULL # dejar solo las variables biofisicas
gini_AP_df<-as.data.frame(t(gini_AP_df))
gini_AP_df$PA_VA<-row.names.data.frame(gini_AP_df)
gini_AP_df$PA<-substr(gini_AP_df$PA_VA,1,2) # extraer el pais
gini_AP_df$VA<-substr(gini_AP_df$PA_VA,4,6) # extraer la variable

# tabla dinamica, con el promedio para el pais&variable y el pais
AP_pais_variable <- tapply(as.numeric(gini_AP_df$V1), INDEX=gini_AP_df$PA_VA, FUN=mean)
AP_pais <- tapply(as.numeric(gini_AP_df$V1), INDEX=gini_AP_df$PA, FUN=mean)

write.table(AP_pais_variable, "salidas txt/td_AP_pais_variable_-1990.txt",col.names=T, row.names=T)
write.table(AP_pais, "salidas txt/td_AP_pais_-1990.txt",col.names=T, row.names=T)


# gini_FP
gini_FP_df<-as.data.frame(gini_FP)
gini_FP_df[76:195]<-NULL # dejar solo las variables biofisicas
gini_FP_df<-as.data.frame(t(gini_FP_df))
gini_FP_df$PA_VA<-row.names.data.frame(gini_FP_df)
gini_FP_df$PA<-substr(gini_FP_df$PA_VA,1,2) # extraer el pais
gini_FP_df$VA<-substr(gini_FP_df$PA_VA,4,6) # extraer la variable

# tabla dinamica, con el promedio para el pais&variable y el pais
FP_pais_variable <- tapply(as.numeric(gini_FP_df$V1), INDEX=gini_FP_df$PA_VA, FUN=mean)
FP_pais <- tapply(as.numeric(gini_FP_df$V1), INDEX=gini_FP_df$PA, FUN=mean)

write.table(FP_pais_variable, "salidas txt/td_FP_pais_variable_-1990.txt",col.names=T, row.names=T)
write.table(FP_pais, "salidas txt/td_FP_pais_-1990.txt",col.names=T, row.names=T)