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","GI","PY","PE","UY","VZ","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_humanas",".PNG", sep=""), width=1400, height=3000, units = "px")
par(mfrow=c(12,5),mar=c(2,3,2,3), cex.axis=2, xpd=T)
pais_sel <- c(21:26,34,30,29,32,33,35)
var_sel <- c(3:7)  ####  Las variables que me interesan por posicion

for (j in pais_sel){ 
  dat_pais <- data.frame(dat[,j],dat[,j+15],dat[,c(20,17,13,14,16)])
  names(dat_pais) <- c(names(dat)[j], names(dat)[j+15], names(dat)[c(20,17,13,14,16)])
  dat_pais <- subset(dat_pais, dat_pais[,1]>0)
  
  for (i in var_sel){
    dat_var <- data.frame(dat_pais[,1:2], dat_pais[,i], stringsAsFactors = FALSE);
    names(dat_var) <- c(names(dat_var)[1:2], names(dat_pais)[i]) 
    print(head(dat_pais))
  
    limINF <- quantile(dat_var[,3], 0.025)  
    limSUP <- quantile(dat_var[,3], 0.975)  
  
    dat_var <- cbind(dat_var, simp = ifelse(dat_var[,3] >= limSUP,limSUP,
                  ifelse(dat_var[,3]<=limINF, limINF, dat_var[,3])))
 
    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]

    # esto lo tengo que revisar si quisiera poner los nombres
    # nom_var<-substr(names(dat_var)[j],6,7) # aca del campo AREA_PAIS se extrae el nombre del pais
    # print(nom_var)
    # 
    # nom_pais<-names(dat)[i]
    # print(nom_pais)
  
    histograma = hist(dat_var$simp, breaks=cla, plot=F, include.lowest=T)
    cat = cut(dat_var$simp, breaks=histograma$breaks, include.lowest=TRUE,labels=F)
    dat_var  = cbind(dat_var, 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_var[,1]), INDEX=dat_var$cat, FUN=sum)
    areasPclases <- tapply(as.numeric(dat_var[,2]), INDEX=dat_var$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(1000,1000, "textoxxxxxxxxxxxxxxxxxxxxxxxxxxxx") NO FUNCIONA
    } 
}
dev.off()