#****************************************
# Creado por Magdiel Ablan  24-09-04
# Modificado por Mary Pili Vargas 20-06-05 
#****************************************

# Este programa verifica la homogeneidad de las series formadas con los valores de 
# precipitación, utilizando el método de Thom (Thom H. 1966."Some Methods of Climatological
# Analysis". Secretaria of the World Meteorological Organization. Geneva. Switzerland)

# tabla_homogeneidad es la tabla estadística que determina si la serie es homogénea o no

tabla_homogeneidad <- read.table("./Archivos_txt/tabla.homogeneidad.txt",header = T, sep = ",", dec = ".", as.is = F, na.strings = "NA")

# El número de estaciones consideradas es nest
nest<-length(datos$Serial.C)

# Este es el encabezado del archivo texto que guarda los resultados.
	
cat("Nombre de la Estación", "Longitud de la serie","No. de series","Mediana","Homogenea?","\n",sep=",",file="./Resultados/homogeneidad")

# Habilitamos el archivo gráfico


pdf("./Graficos/ANOMALIAS.pdf")

# df.nom.cr contiene el nombre de las estaciones corregidas

# Repita para todas las estaciones:
for (j in 1:nest){


	#	Cargue los datos en est.clave
	marco<-get(df.nom.cr[j])

	cat("-----------------------------------------------\n")
	cat("Procesando la estacion ", df.nom.cr[j],"\n")
	cat("-----------------------------------------------\n")

	# Establezca el marco de datos de trabajo
	attach(marco)


	# Se supone que cualquier valor faltante ya ha sido reemplazado por su
	# estimado. Por ahora utilizar na.rm=T para garantizar que no hay
	# valores faltantes

	# homo es la variable que indica si la serie es o no homogénea
	homo<-""

    	# tome en consideración sólo los valores a partir de 1960
    	precipitacion<-marco$PRE[marco$A>1959]
    	anos<-marco$A[marco$A>1959]

	# Verifique homogeneidad a nivel anual
	tot.anual<-aggregate(precipitacion,list(anos),FUN=sum)
	PREanual<-tot.anual$x

	# Realice un gráfico temporal para ver si es posible detectar alguna tendencia
	# solo por inspección visual

	# x11()
	# Graficamos las anomalías
	mediana<-median(PREanual,na.rm=T)
	anoma<-(PREanual-mediana)/sd(PREanual)
	marco.ts <- ts(anoma,start=marco$A[1])

	ts.plot(marco.ts,col="blue")
	abline(mediana,0)
	title(df.nom.cr[j],col.main="black")

	# longitud total de la serie
	tpre<-length(PREanual)

	# Inicializa el vector
	seriepre<-c(1:tpre)

	# seriepre es 1 si el valor es menor a la mediana
	# 2 si el valor es igual y 3 si el valor es mayor
	# El siguiente lazo asigna los valores correspondientes
	# a seriepre
	for (i in 1:tpre){
		if (PREanual[i] < mediana)
			seriepre[i]<-1
		else
			if (PREanual[i]==mediana)
				seriepre[i]<-2
			else
				seriepre[i]<-3
	}

	# Ahora vamos a contar el número de series
	i<-1

	# Los Ni son el número de series con valor menor, igual o superior
	# a la mediana
	N1<-0
	N2<-0
	N3<-0

	# Se recorre ahora toda la serie
	while (i <= tpre){
		# Los Li representan la longitud de las series
		L1<-0
		L2<-0
		L3<-0

		# Si estamos en una serie de valores menores a la mediana
		if (seriepre[i]==1)
			while ((seriepre[i]==1) & (i<=tpre))
			{
				L1 <-L1+1
				i<-i+1
			}
		else
		# Si estamos en una serie de valores iguales a la mediana
			if (seriepre[i]==2)
				while((seriepre[i]==2)&(i<=tpre))
				{
					L2 <-L2+1
					i<-i+1
				}
			else
		# La serie es superior a la mediana
				while ((seriepre[i]==3)&(i<=tpre))
				{
					L3 <-L3+1
					i<-i+1
				}
	# Este conjunto de if identifica de que tipo de serie se trata
	if (max(L1,L2,L3)==L1) N1<-N1 + 1
	if (max(L1,L2,L3)==L2) N2<-N2 + 1
	if (max(L1,L2,L3)==L3) N3<-N3 + 1
}

# Si hay series iguales a la mediana su número se reparte entre los
# dos anteriores

if (N2!=0){
	N1 <-N1 + round(N2/2,0)
	N3<- N3 + round(N2/2,0)
}

#  u es el número total de series
	u <- N1 + N3


# Se reportan las series por debajo y por arriba de la mediana
cat("Longitud de la serie = ",tpre,"\n")
cat("La mediana es ",mediana,"\n")
cat("Menor a la mediana = ",N1,"\n")
cat("Mayor a la mediana = ",N3,"\n")
cat("No. total de series u = ",N1 + N3,"\n")

# Se buscan los resultados en la tabla estadística
	indice<-round(tpre/2,0)
	refer<-tabla_homogeneidad[tabla_homogeneidad$NAn==indice,]
	if ((u >= refer$P10) & (u <= refer$P90)) {
		cat("ES HOMOGENEA \n")
		homo<-"si"}
	else {
		cat("NO ES HOMOGENEA \n")
		homo<-"no" }

# Se anotan los resultados en un archivo de texto, homogeneidad
cat(df.nom.cr[j],round(tpre/2),u,mediana,homo,"\n",sep=",",file="./Resultados/homogeneidad",append=T)


detach("marco")
}
dev.off()