#****************************************
# Creado por Magdiel Ablan  13-08-03
# Modificado por Mary Pili Vargas 01-03-05 
#****************************************

# Este algoritmo lee los valores de precipitación desde un data frame, en formato
# de 4 columnas, Año, Mes, Precipitación, Código; esta última columna indica el estatus
# del valor, es decir si estos valores de precipitación son posibles valores extremos o
# si es un valor estimado utilizando cualquiera de los métodos de estimación.

# Para el proyecto fueron consideradas 68 estaciones, ubicadas en los estados llaneros, 
# cada una de estas estaciones tiene un archivo asociado por lo que su marco de datos contiene 
# la  información mencionada anteriormente

# También se hace uso de la función outlier, esta función permite detectar los posibles
# valores extremos de precipitación. Esta función indica si hay valores extremos, si la
# desviación estándar es mayor a 4 veces la mediana de la desviación absoluta o mad un
# estimador  robusto de dispersión, este valor umbral de cuatro puede ser discutido.

# Si se detectan valores extremos se procede a generar un gráfico que muestra el año donde
# se ubica este valor, la ubicación de la media y la ubicación de 3 desviaciones estándares.
# Se consideran como posibles valores extremos aquellos valores que exceden este umbral.
# Este umbral, puede ser discutido y puede cambiarse.

# Función outlier. Esto es necesario ejecutarlo sólo la primera vez, ya que esta función
# luego de ser ejecutada se convierte en una función de R

outlier<-function(x) {
        if (sqrt(var(x)) > 4*mad(x)) T
        else F
}

# Se establece el marco de trabajo, este nombre cambiará para cada una de las
# estaciones

marco.nom <- readline("Nombre de la estacion ")
marco <- get(marco.nom)

# Esta función permite habilitar el marco de trabajo
attach(marco)

# Este el tipo de archivo donde se guarda el gráfico asociado a cada estación
jpeg("Estacion.jpeg")

# también podrían ser archivos pdf, ps ó png, dependiendo del formato que desee el usuario
# pdf("Estacion.pdf")
# postscript("Estacion.ps")
# png("Estacion.png")

# Repita para todos los meses
for(i in 1:12) {

	# La función  na.omit() permite extraer los valores omitiendo las filas del marco de datos
	# que estén incompletas
	# Extracción del subconjunto de los valores del mes i:

  	mes <- na.omit(marco[marco$M==i,])

 	# Detección de valores extremos

  	if(outlier(mes$PRE)) {

	# Se generan los gráficos donde hay posibles valores extremos

		plot(mes$A,mes$PRE,xlab="Año",ylab="PRECIPITACION (mm.)",pch=22,col="yellow",bg="red",col.axis="black")
   	 	title(marco.nom,col.main="black")
   	 	subtit <- as.character(i)
    		mtext(subtit,3)
    		media <- as.character("Media")
    		text(1999,(mean(mes$PRE)+2),labels=media,adj=2,col="black",font=2)
    		abline(h=mean(mes$PRE),col="red")
    		desv <-as.character("3 Dev Estandares")
    		text(1999,(mean(mes$PRE)+3*sqrt(var(mes$PRE))+2),labels=desv,adj=1.3,col="black",font=2)
    		abline(h=mean(mes$PRE) + 3*sqrt(var(mes$PRE)),lty=2,col="blue")
    	
	   	ext <- mes[mes$PRE > 3*sd(mes$PRE),]

 		# PE se escribe en la cuarta columna del marco de datos para indicar que se encontró un  posible valor extremo
   
   	 	PE <- as.character("PE")
    		marco[row.names(ext),4] <- PE
   	 	text(ext$A,ext$PRE,labels=ext$A,adj=1.2)
    		cat("Valores extremos en el año ",ext$A," mes ",i,"\n")
    		
 	 } # Fin del si
    else
    cat("no hay valores extremos \n")
    
 }# Fin del for

dev.off() # cerramos el archivo gráfico
detach(marco)