rm(list=ls())
library(ineq); library(lattice); library(png); 
library(gridExtra); library(ggplot2)
setwd('D:\\areas_protegidas\\xls_dbf')
# setwd("/home/marcos/Escritorio/Dropbox/gb")
dat<-read.table('todas_las_variables_sam_simp_fk_1000km.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

dat$DIST_ROAD_km <- dat$DIST_ROAD/1000
dat$DIST_FRON_km <- dat$DIST_FRON/1000

dat$DIST_ROAD <- dat$DIST_ROAD_km
dat$DIST_FRON <- dat$DIST_FRON_km

## POR PAIS, 21 a 34 Argentina a las tres guayanas juntas ------------------------------------------------

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

## ACA ESTA LO QUE MODIFIQUE....
png(file=paste("D:\\areas_protegidas\\figuras\\Sudamerica\\histogramas",".PNG", sep=""), width=3000, height=1250, units = "px")
par(mfrow=c(5,12),mar=c(2,3,2,3), cex.axis=2, xpd=NA)
var_sel <- c(20,17,13,14,16)  ####  Las variables que me interesan por posicion
pais_sel <- c(2:7,15,11,10,13,14,16)
# pais_sel <- c(1,2,3)
for (j in var_sel){ 
  dat_var <- data.frame(dat[,j], dat[,21:50]);
  for (i in pais_sel){
    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

    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("Area (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="Proteccion (%)",side=4, line=3.2, col="red3", cex=2)
    
    } 
}
# text(10,11, "textoxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
dev.off()
