# Walter & Lieth diagram plot for R
# From R package climatol <http://stat.ethz.ch/CRAN/web/packages/climatol/index.html>
# © 2009 José A. Guijarro <jaguijarro@inm.es>
# © 2010 Xavier Robin <http://xavier.robin.name>
#
# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with this program.  If not, see <http://www.gnu.org/licenses/>.
# Usage:
# diagwl(dat, ...)
# with dat a 2x12 matrix, first row the precipitation (in mm), second row the temperature (in °C) for the 12 months of the year.
# for more details, see diagwl function in package climatol <http://stat.ethz.ch/CRAN/web/packages/climatol/index.html>


#Diagrama de Walter y Lieth
diagwl <- function( dat, est="", alt=NA, per="", margen=c(2,5,5,5),
    mlab="", pcol="#005ac8", tcol="#e81800", pfcol="cyan",
    sfcol="#0eb6d7", drycol"yellow", hwcol="mediumspringgreen", ...) {
  plot.new()
  old.par <- par(no.readonly = TRUE)
  on.exit(par(old.par, new=FALSE))
  par(mar=margen, pty="s", las=1, new=TRUE)
  nr <- nrow(dat) #no. de filas de datos mensuales
  #etiquetas de los meses
  if(mlab=="es") mlab=c("E","F","M","A","M","J","J","A","S","O","N","D")
  else if(mlab=="en") mlab=c("J","F","M","A","M","J","J","A","S","O","N","D")
  else mlab=c(1:12) #etiquetas numéricas
  dat <- as.matrix(dat)
  p <- dat[1,] #precipitaciones medias mensuales
  tm <- dat[2,]
  pmax <- max(p) #precipitación máxima
  ymax <- 60  #máxima ordenada por defecto
  if(pmax > 300) ymax <- 50 + 10*floor((pmax+100)/200)
  ymin <- min(-1.5,min(tm)-1.5) #mínima ordenada sin redondear
  #ejes:
  if(ymin < -1.5) {
    ymin=floor(ymin/10)*10 #mínima ordenada redondeada
    labT <- paste(ymin)
    labP <- ""
    if(ymin < -10) {
      for(i in (ymin/10+1):-1) {
        labT <- c(labT,i*10)
        labP <- c(labP,"")
      }
    }
    labT <- c(labT,"0","10","20","30","40","50","")
    labP <- c(labP,"0","20","40","60","80","100","300")
  }
  else {
    labT <- c("0","10","20","30","40","50","")
    labP <- c("0","20","40","60","80","100","300")
  }
  if(ymax > 60) {
    for(i in 6:(ymax/10-1)) {
      labT <- c(labT,"")
      labP <- c(labP,100*(2*i-7))
    }
  }
  plot(0:13-0.5,c(tm[12],tm[1:12],tm[1]),xlim=c(0,12),ylim=c(ymin,ymax),type="l",col=tcol,xaxs="i",yaxs="i",xaxp=c(0,12,12),xlab="",ylab="",xaxt="n",yaxt="n",bty="n",...)
  lmin <- ymin #mínima ordenada a rotular
  if(lmin==-1.5) lmin=0
  axis(2,((lmin/10):(ymax/10))*10,labels=labT,col.axis=tcol)
  axis(4,((lmin/10):(ymax/10))*10,labels=labP,col.axis=pcol)
  mtext("°C",2,col=tcol,las=1,line=3,adj=0,at=55)
  mtext("mm",4,col=pcol,las=1,line=3,adj=1,at=55)
  abline(0,0)
  abline(50,0)
  #rótulos:
  if(is.na(alt)) mtext(est,line=2,adj=0)
  else mtext(paste(est," (",alt," m)",sep=""),line=2,adj=0)
  mtext(per,line=1,adj=0)
  mtext(paste(round(mean(tm*10))/10,"°C    ",round(sum(p))," mm",sep=""),line=1,adj=1)
  x <- 0:13-0.5
  p2 <- c(p[12],p[1:12],p[1])
  plotted.pmax <- FALSE
  plotted.pmin <- FALSE
  if(pmax<=100) {
    xl <- x
    yl <- c(p[12],p[1:12],p[1])/2
    n2 <- 14
  }
  else { #cambio de escala en prec. > 100
    xp <- numeric(30)
    yp <- numeric(30)
    xl <- numeric(25)
    yl <- numeric(25)
    n <- 0
    n2 <- 0
    gr <- FALSE
    if(p2[1]>100) { #primer punto
      n <- n+1
      xp[n] <- x[1]
      yp[n] <- 50
      n <- n+1
      xp[n] <- x[1]
      yp[n] <- 50 + (p2[1]-100)/20
      n2 <- n2+1
      xl[n2] <- x[1]
      yl[n2] <- 50
      gr <- TRUE
    }
    else {
      n2 <- n2+1
      xl[n2] <- x[1]
      yl[n2] <- p2[1]/2
    }
    for(i in 2:14) {  #demás puntos
      if(gr) {  #si p anterior > 100
        n <- n+1
        if(p2[i]>100) {
          xp[n] <- x[i]
          yp[n] <- 50 + (p2[i]-100)/20
        }
        else {
          xp[n] <- x[i-1] + (100-p2[i-1])/(p2[i]-p2[i-1])
          yp[n] <- 50
          n2 <- n2+1
          xl[n2] <- xp[n]
          yl[n2] <- 50
          n <- n+1
          xp[n] <- NA
          yp[n] <- NA
          n2 <- n2+1
          xl[n2] <- x[i]
          yl[n2] <- p2[i]/2
          gr <- FALSE
        }
      }
      else {  # p anterior <=100
        if(p2[i]>100) { #si p > 100
          n <- n+1
          xp[n] <- x[i-1] + (100-p2[i-1])/(p2[i]-p2[i-1])
          yp[n] <- 50
          if(xl[n2]!=x[i-1]) {  #no repetir puntos!
            n2 <- n2+1
            xl[n2] <- x[i-1]
            yl[n2] <- p2[i-1]/2
          }
          n2 <- n2+1
          xl[n2] <- xp[n]
          yl[n2] <- 50
          n <- n+1
          xp[n] <- x[i]
          yp[n] <- 50 + (p2[i]-100)/20
          gr <- TRUE
        }
        else { # p <=100
          n2 <- n2+1
          xl[n2] <- x[i]
          yl[n2] <- p2[i]/2
        }
      }
    }
    if(!is.na(yp[n])) {  #cerrar último polígono
      n <- n+1
      xp[n] <- xp[n-1]
      yp[n] <- 50
      n2 <- n2+1
      xl[n2] <- 12.5
      yl[n2] <- 50
    }
    polygon(xp[1:n],yp[1:n],col=pcol,border=pcol)
    mtext(formatC(max(as.matrix(dat[1,])),digits=1,format="f"),4,las=1, line=2.5,at=max(yp[1:n], na.rm=TRUE))
    plotted.pmax <- TRUE
    if (min(as.matrix(dat[2,])) > 100) {
      mtext(formatC(min(as.matrix(dat[1,])),digits=1,format="f"),4,las=1, line=2.5,at=min(yp[1:n], na.rm=TRUE))
      plotted.pmin <- TRUE
    }
  }
  #tramas:
  pi <- approx(xl[1:n2],yl[1:n2],n=66)$y
  ti <- approx(x,c(tm[12],tm[1:12],tm[1]),n=66)$y
  d <- pi - ti
  xi <- (1:66)/5-0.7
  xw <- subset(xi,d>0) #periodo húmedo
  y1 <- subset(pi,d>0)
  y2 <- subset(ti,d>0)
  if(length(xw)>0) segments(xw,y1,xw,y2,col=pcol,lty=1,lwd=1)
  xw <- subset(xi,d<0) #periodo seco
  y1 <- subset(pi,d<0)
  y2 <- subset(ti,d<0)
  if(length(xw)>0) segments(xw,y1,xw,y2,col=tcol,lty=3,lwd=2)
  #curvas de P y T:
  lines(xl[1:n2],yl[1:n2],col=pcol,lwd=2)
  lines(x,c(tm[12],tm[1:12],tm[1]),col=tcol,lwd=2)
  #yearly max temp
  mtext(formatC(max(as.matrix(dat[2,])),digits=1,format="f"),2,las=1, line=2.5,at=max(as.matrix(dat[2,])))
  #yearly min temp
  mtext(formatC(min(as.matrix(dat[2,])),digits=1,format="f"),2,las=1, line=2.5,at=min(as.matrix(dat[2,])))
  # yearly min rain
  if (!plotted.pmin)
    mtext(formatC(min(as.matrix(dat[1,])),digits=1,format="f"),4,las=1, line=2.5,at=min(dat[1,])/2)
  if (!plotted.pmax)
    mtext(formatC(max(as.matrix(dat[1,])),digits=1,format="f"),4,las=1, line=2.5,at=max(dat[1,])/2)
  #marcar límites de los meses:
  for(i in 0:13) segments(i,ymin+1.5,i,ymin)
  abline(h=ymin+1.5)
  #heladas seguras
  for(i in 1:12) if(dat[2,i]<=0) rect(i-1,ymin,i,ymin+1.5,col=sfcol)
  #heladas probables
  for(i in 1:12) if(dat[2,i]<=10 & dat[2,i] > 0) rect(i-1,ymin,i,ymin+1.5,col=pfcol)
  # Hot dry months
  for(i in 1:12) if(dat[2,i]>10 & dat[1,i] <= dat[2,i]*2) rect(i-1,ymin,i,ymin+1.5,col=drycol)
  # Hot humid months
  for(i in 1:12) if(dat[2,i]>10 & dat[1,i] > dat[2,i]*2) rect(i-1,ymin,i,ymin+1.5,col=hwcol)
  #rótulos meses:
  mtext(mlab,1,las=1,line=0.5,adj=0.5,at=x[2:13])
  #reset old.par (restablecemos parámetros gráficos anteriores):
  invisible()
}