Guia Para Usar Programas r y Clima

8/17/2019 Guia Para Usar Programas r y Clima

http://slidepdf.com/reader/full/guia-para-usar-programas-r-y-clima 1/44



II

Preámbulo

El paquete de R “Climatol” está mayoritariamente destinado al problema de la homogeneiza-ción de series climatológicas, es decir, a eliminar las perturbaciones producidas por cambiosen las condiciones de observación o en el entorno de la estación, para que las series reflejensólamente (hasta donde sea posible) las variaciones climáticas.

La documentación estándar del paquete se ciñe a las normas de R, y provee descripciones delas funciones y de sus parámetros, de modo que los usuarios pueden acudir a ella cuando lonecesiten. Esta guía, por otra parte, se ha escrito como un complemento, y está más enfocada a

explicar la metodología subyacente en los algoritmos del paquete, cómo llamar a sus funciones,y cómo interpretar y usar sus resultados.

La guía está estructurada en dos partes: una introducción rápida (en unas pocas páginas siguien-tes) para aquellos usuarios que deseen empezar a homogeneizar sus datos cuanto antes, y unaguía ampliada, en la que se tratan con más detalle los diferentes aspectos del paquete.

La mayor parte de los ejemplos de esta guía se pueden reproducir con los ficheros del archivoclimatol-dat.zip, que se puede descargar de http://webs.ono.com/climatol/climatol-dat.zip,y que contiene series reales de un área mediterránea, si bien los nombres y coordenadas de lasestaciones son ficticios.

Agradecimientos

Este paquete se ha beneficiado enormemente de las fructiferas discusiones mantenidas en elmarco de la Acción COST ES0601 (2006-2011), titulada Avances en los métodos de homoge-

neización de las series climáticas: una aproximación integrada (HOME). Mi agradecimiento a

todos los participantes, así como a la Fundación Europea de la Ciencia, por promover y financiarestos enriquecedores encuentros. También debo agradecer a la Agencia Estatal de Meteorologíade España (AEMET) por su continuado apoyo a mi participación en esta Acción.



II I

Introducción rápida

Lo primero que hemos de hacer es preparar los datos de entrada en dos ficheros de texto conlos formatos adecuados. En uno de ellos hay que relacionar las coordenadas y nombres de las

estaciones, incluyendo una línea de la formaX Y Z CÓDIGO NOMBRE

para cada estación, donde las coordenadas X e Y pueden estar en km (procedentes, por ejemplo,de una proyección UTM) o en coordenadas geográficas (longitud y latitud, en este orden), perono en forma de grados, minutos y segundos, sino en grados con decimales. Los otros paráme-tros son la altitud Z en m, un CÓDIGO identificativo de la estación, y su NOMBRE completo, quedebe estar encerrado entre comillas si está formado por más de una palabra. (Es aconsejableponer todos los nombres entre comillas para evitar errores). El nombre de este fichero debeser VAR_AINI-AFIN.est, donde VAR será una abreviatura de la variable climática que estemosanalizando, y AINI y AFIN los años inicial y final del periodo estudiado.

Los datos climáticos de esta variable irán en otro fichero, organizados por bloques, estación porestación, en el mismo orden en que aparecen en el fichero de estaciones. El nombre de ambosficheros será el mismo, distinguiéndose únicamente por su extensión, que en el caso del ficherode datos será dat.

Ejemplo: Supongamos que se quieren homogeneizar los datos mensuales medios de las tempe-raturas mínimas diarias de 1956 a 2005, y que se escoge Tmin como la abreviatura para estavariable. El fichero de estaciones sería Tmin_1956-2005.est, y podría comenzar, como en losdatos de ejemplo, por:

27.0 53.9 456 S03 "La Perla"31.8 26.5 123 S08 "El Palmeral"49.2 30.0 154 S11 "Miraflores"43.4 29.6 156 S13 "Torremar"... (etc)

Y el fichero de datos debería llamarse Tmin_1956-2005.dat, y sus primeras líneas podrían ser:

NA NA NA NA NA NA NA NA NA NA NA NA

-0.4 1.8 5.5 6.5 15.1 17.4 16.7 16.4 12.2 6.0 2.6 2.31.5 4.0 6.5 8.7 12.4 12.1 20.3 NA 14.7 11.0 3.2 0.5... (etc)

Estos serían los datos de la primera estación de nuestra red de observación1, en orden cronológi-co: enero a diciembre de 1956, lo mismo para 1957 en la segunda línea, 1958 en la tercera, etc.En este ejemplo faltan los datos de todo 1956 y de agosto de 1958, que se han substituido porNA ( Not Available), que es la representación estándar en R de los datos ausentes (aunque puedenusarse otras). Después de relacionar todos los datos de la primera estación, se continúa con los

1En realidad, este no es el comienzo de nuestro fichero de ejemplo, cuyas tres primeras líneas tienen todos losdatos completos. Aquí hemos introducido estas otras para ilustrar cómo proceder cuando nos falten datos.



IV

de la segunda, y así sucesivamente hasta completar los datos de la última estación. Es impor-tante tener en cuenta que todas las estaciones deben proveer datos para todos y cada uno de lostérminos (meses, estaciones del año, o la unidad temporal que estemos tratando) del periodo deestudio (1956-2005 en nuestro ejemplo), y de ahí la necesidad de incluir códigos para rellenarcualquier dato ausente. Por comodidad, hemos puesto 12 valores (un año completo) en cada

línea del fichero, pero los datos se pueden disponer de cualquier otro modo, en un formato libre,separados por espacios en blanco, puesto que se van a leer secuencialmente. (Nota importante:ningún término temporal (mes, etc) debe faltar simultáneamente en todas las estaciones, puestoque el proceso de relleno de datos ausentes no podría realizarse completamente y el programaacabaría dando un error).

Una vez preparados los ficheros de datos en nuestro directorio de trabajo, todo lo que tenemosque hacer para proceder a su homogeneización es arrancar R desde ese mismo directorio, cargarlas funciones de homogeneización con la orden

library(climatol)

Si se instaló este paquete desde R, o con

source("depurdat.R")

si se dispone de este fichero2 en el directorio de trabajo, y lanzar la orden de homogeneizaciónautomática, que para nuestro ejemplo sería:

homogen("Tmin", 1956, 2005)

Esta orden acepta otros parámetros opcionales, de los cuales cabe destacar los siguientes:

nm Número de datos por año en cada estación (12 por defecto: datos mensuales. Poner nm=1 sianalizamos datos anuales, nm=1 para datos estacionales, etc).

deg Ponerlo igual a TRUE (verdadero) si las coordenadas geográficas están en grados, o dejarloen su valor por defecto FALSE si están en km (la unidad de distancia usada internamentepor el paquete).

std Tipo de normalización. Por defecto, los datos se estandarizarán restándoles su media ydividiendo el resultado por su desviación típica, pero si la variable tiene un cero natural(como la precipitación), puede ser preferible usar std=2 (los datos sólo se dividirán porsu media). Otra opción es std=1, para que a los datos sólamente se les reste su media).

rtrans Transformación raíz a aplicar a los datos: 2 para raíz cuadrada, 3 para cúbica, etc(pueden usarse números no enteros). Útil si la distribución de la variable se aleja de lanormal, como sucede con la velocidad del viento, o con la precipitación de regiones ári-das).

na.strings Cadena de caracteres usada para los datos ausentes. Por defecto R usa ’NA’,pero se puede especificar cualquier otra, como por ejemplo: na.strings=’-999.0’.

Otro ejemplo para homogeneizar precipitaciones estacionales (cuatro datos por año) para elperiodo 1961-2005, con las coordenadas de las estaciones expresadas en grados geográficos, y

2El fichero depurdat.R contiene las funciones de homogeneización del paquete climatol.



V

aplicando una transformación raíz cúbica a los datos, será (no se proveen ficheros de datos paraeste ejemplo):

homogen("SsPrp", 1961, 2005, nm=4, deg=TRUE, rtrans=3)

La orden del primer ejemplo generaría los siguientes ficheros (en el directorio de trabajo):

Tmin_1956-2005.esh Fichero de estaciones después de la homogeneización. Tiene lamisma estructura que el fichero de entrada Tmin_1956-2005.est, pero con columnasadicionales (ver la guía ampliada) y, probablemente, nuevas líneas (cuando el procesodetecta un salto brusco en la media, la serie se corta, creando una nueva con las mis-mas coordenadas y añadiendo un sufijo numérico incremental al nombre y código de laestación).

Tmin_1956-2005.dah Datos homogeneizados, con todos los datos ausentes rellenados,análogo al fichero de entrada Tmin_1956-2005.dat.

Tmin_1956-2005.txt Fichero de bitácora del proceso, con todos los mensajes que hanido saliendo por pantalla (incluyendo los resúmenes finales).

Tmin_1956-2005.pdf Fichero con una (potencialmente larga) colección de gráficos dediagnóstico generados durante el proceso.

Los archivos de gráficos y de bitácora pueden sugerir repetir el proceso con diferentes pará-metros (ver la guía ampliada para más información), mientras que los ficheros con los datoshomogeneizados se pueden tratar con la función dahstat. Por ejemplo, si queremos una re-lación de los valores medios para el periodo 1971-2000 de las temperaturas que acabamos dehomogeneizar, podemos obtenerla en un archivo llamado Tmin_1971-2000.med con la orden:

dahstat("Tmin", 1956, 2005, 1971, 2000)

Como puede observarse, los parámetros son el nombre de la variable, el primer y el último añodel periodo de estudio, y el primer y último año del periodo para el que queremos calcular lasmedias (estos pueden omitirse si queremos las medias de todo el periodo de estudio). Otrosparámetros de esta función son:

out Tipo de salida (el fichero tendrá la extensión correspondiente):"med" para medias de los datos (la salida por defecto)."mdn" para medianas.

"max" para valores máximos."min" para valores mínimos."std" para desviaciones típicas."q" para cuantiles (ver el parámetro prob)."tnd" para tendencias.Con cualquier otra opción no reconocida la función se limitará a leer los datos homoge-neizados, facilitando así al usuario su posterior análisis.

vala Valor anual calculado en la tabla de salida. Puede ponerse 0 (para no calcular ningúnvalor anual), 1 (para la suma de los valores mensuales of de la periodicidad subanual queestemos manejando), 2 (para la media, que es la opción por defecto), 3 (para el máximo)o 4 (para el mínimo).



VI

prob Probabilidad para el cálculo de los cuantiles (si se usa la opción out="q" . El valor pordefecto es 0.5, equivalente al cálculo de la mediana, como con la opción out="mdn" ).

eshcol Columnas del fichero de estaciones homogeneizadas "*.esh" a incluir en el ficherode salida. El valor por defecto es 4, indicando que sólo el código de la estación (la cuarta

columna) precederá a los valores estadísticos de la tabla.

Los ficheros de salida tendrán el mismo nombre base que los de entrada, pero con una extensiónigual a la opción out escogida, con la excepción de los cuantiles, cuya extensión será qPP, dondePP se substituirá por la probabilidad seleccionada con la opción prob (pero en%).

Por consiguiente, si queremos obtener los valores normales mensuales (y su media anual) de lastemperaturas mínimas previamente homogeneizadas, podemos dar la siguiente orden:

dahstat("Tmin", 1956, 2005, 1971, 2000)

Pero si deseamos calcular las tendencias para todo el periodo de estudio 1956-2005, e incluir lascoordenadas de las estaciones (columnas 1 and 2 del fichero de salida Tmin_1956-2005.eshtras los códigos de las estaciones, deberíamos dar:

dahstat("Tmin", 1956, 2005, out="tnd", vala=1, eshcol=c(4,1,2)) 3

y de este modo obtendríamos la lista de las tendencias en un fichero de texto denominadoTmin_1956-2005.tnd que, al incluir las coordenadas de las estaciones, podría utilizarse paracartografiar las tendencias (tanto desde R como importando el fichero con un SIG).

(Fin de la introducción rápida)

3Nótese el uso de la función de concatenación de R, c, para designar un vector numérico.



VII

Guía ampliada

Índice

1. Introducción 1

2. Metodología 1

2.1. Regresión tipo II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2.2. Estimación de los datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.3. Detección y correción de datos anómalos y saltos bruscos en la media . . . . . 4

3. Aplicación 53.1. Preparación de los datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.2. Homogeneización de las series . . . . . . . . . . . . . . . . . . . . . . . . . . 5

4. Salidas 8

4.1. El fichero *.txt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

4.2. Archivo *.pdf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4.3. Ficheros *.esh

y *.dah

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5. Discusión y sugerencias 21

6. Explotación de las series homogeneizadas 23

7. Y si los datos son diarios o sub-diarios? 25

8. Otras funciones de climatol 28

8.1. Rosas de los vientos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 288.2. Climogramas de Walter y Lieth . . . . . . . . . . . . . . . . . . . . . . . . . . 29

9. Bibliografía 31

10. Anexo: Valores umbrales para las pruebas SNHT 32



1

1. Introducción

Como el lector debe saber, las estaciones meteorológicas no sólo registran las variaciones climá-ticas locales, sino que sus medidas también están afectadas por cambios en la instrumentación,

métodos de observación, la ubicación del observatorio, o su ambiente circundante (ej.: creci-miento urbano o cambios en los usos del suelo). Para filtrar las series de esas perturbacionesno deseadas es preciso aplicar un procedimiento estadístico, denominado homogeneización, demodo que al final sus datos reflejen únicamente las variaciones del clima.

Esta problemática se conoce desde hace muchos decenios. Algunos métodos antiguos se basanen pruebas estadísticas para comprobar la no estacionariedad de una única serie climatológica.Estos métodos absolutos deben evitarse, puesto que presuponen una estabilidad climática quese ha visto que no es realista. La alternativa es usar métodos de homogeneidad relativa, enlos que las pruebas de estacionariedad se aplican a series de razones o diferencias entre laestación problema y una o más series bien correlacionadas de estaciones vecinas. Peterson et

al. (1998) y Aguilar et al. (2003) pasan revista a las diferentes aproximaciones desarrolladas porlos climatólogos hasta ahora, mientras que en los próximos apartados explicaremos la estrategiseguida en este paquete.

2. Metodología

2.1. Regresión tipo II

Como en muchos otros métodos, las pruebas de homogeneidad se aplican aquí a series de di-ferencias entre la estación problema y una serie de referencia construida mediante una media(ponderada o no) de las series de las estaciones de las proximidades. Pero, a diferencia de lamayoría de ellos, la selección de estas estaciones se basa únicamente en la proximidad, y no enla correlación, para permitir el uso de las estaciones más próximas incluso cuando el periodocomún de observación es demasiado pequeño (o nulo) para poder calcular correlaciones fiables.Por tanto, mientras que el uso de las correlaciones se suele aplicar a una selección de serieslargas, nosotros podemos emplear la mayor parte de la información de nuestra red climatoló-gica. Esto implica, sin embargo, que la región estudiada debe ser climáticamente homogénea4,puesto que la presencia de abruptas fronteras geográficas pueden conducir a usar estacionespróximas pero pobremente correlacionadas para calcular las series de referencia. En este caso,

la regió debe subdividirse, para aplicar una homogeneización independiente en cada una de lassubregiones.

Este procedimiento se inspira en el método usado por Paulhus y Kohler (1952) para rellenardatos ausentes de precipitación diaria, que consiste en la interpolación espacial de precipita-ciones relativas (divididas por la precipitación normal) de estaciones vecinas. Este método delas proporciones se extiende en el paquete climatol con opciones para usar diferencias y es-tandarizaciones propiamente dichas para normalizar los datos. Las proporciones respecto a losvalores normales climatológicos son apropiados para la precipitación y otras variables que nopueden tener valores negativos y que suelen tener una distribución de probabilidad en forma de

4O, al menos, que las variaciones climáticas sean suaves.



2

L, mientras que las diferencias (o las estandarizaciones, si estas diferencias se dividen por lasdesviaciones típicas) son más apropiadas para la temperatura y otras variables que se distribuyennormalmente (o se aproximan a la normal).

Desde un punto de vista estadístico, esto equivale a aplicar una regresión lineal tipo II, en

lugar de el tipo I habitual. Éste se suele ajustar por un procedimiento de mínimos cuadrados,que minimiza las desviaciones entre los puntos (observaciones) y la recta de regresión en ladirección del eje Y (verticalmente, como en la figura 1-izquierda). Se supone con ello que lavariable independiente X está controlada por el investigador o medida con errores despreciablesrespecto a los de la variable dependiente (Sokal y Rohlf, 1969). Pero este no es el caso cuandose ajustan líneas de regresión a pares de series de una red climatológica, donde los errores sona priory similares en todas las estaciones. En este caso, las desviaciones a minimizar deberíanser las perpendiculares a la recta de regresión, como en la figura 1-derecha.

−1 0 1 2

− 1 .

0

− 0 .

5

0 .

0

0 .

5

1 .

0

1 .

5

2 .

0

2

. 5

x

y

−1 0 1 2

− 1 .

0

− 0 .

5

0 .

0

0 .

5

1 .

0

1 .

5

2 .

0

2

. 5

x

y

Figura 1: Deviaciones minimizadas por mínimos cuadrados (regresión tipo I, izquierda) y re-gresión ortogonal (regresión tipo II, derecha).

Aunque existe una expresión analítica para el ajuste de esta línea de regresión ortogonal tipoII (Daget, 1979), hay algunas alternativas que proporcionan una buena aproximación. La mássencilla es la denominada eje mayor reducido que, llamando x e y a las versiones estandarizadasde las variables ( x = ( X − m X )/s X y y = (Y − mY )/sY , donde m y s representan la media y la

desviación típica respectivamente), tiene la forma:ˆ y = x

(O bien ˆ y = − x cuando la relación es inversa, que no es el caso cuando tratamos con la mismavariable en una región climáticamente homogénea).

Una característica de esta regresión tipo II es que la varianza de la variable estimada es la mis-ma que la de la original, puesto que esta línea no tiende a la horizontal cuando el coeficiente dedeterminacion (r 2, igual a la fracción de varianza explicada) tiende a cero. Se puede argumentarque, cuando esta fracción es menor que la unidad, la varianza extra proporcionada por la regre-sión de tipo II respecto a la de mínimos cuadrados (tipo I) es espuria, pero cabe esperar altosvalores de r 2 si la red de observación es suficientemente densa, y por otra parte evitaremos el



3

problema que supone una reducción de la varianza cuando el objetivo final de nuestro estudiosea establecer la variabilidad de las series. Además, esta metodología permite corregir no sólolos cambios en la media de las series, sino también posibles cambios en su varianza 5 .

2.2. Estimación de los datos

Una vez normalizados los datos originales, se procede a estimar cada término de cada seriecomo una media ponderada de un determinado número de los datos más próximos disponiblesen cada caso. Los pesos a aplicar a los datos de referencia pueden ser todos iguales (mediasimple) o calcularse como una función inversa de la distancia d entre los sitios de observación.La función escogida para ello se formuló originalmente como 1/(1+ d 2/a), donde el parámetroa permite al investigador modular el peso relativo de las estaciones más cercanas respecto delas más alejadas, pero es más conveniente la expresión 1/(1 + d 2/h2), puesto que de este modoel nuevo parámetro h resulta ser la distancia a la que el peso se reduce a la mitad del que tendría

una estación situada en la misma posición que la de los datos a estimar6. En la figura 2 puedeverse esta función dibujada para diferentes valores de h. (Este parámetro se llama wd, por weight

distance, en la lista de argumentos de la función de homogeneización de este paquete).

0 100 200 300 400

0 . 0

0 . 2

0 . 4

0 . 6

0 . 8

1 . 0

Distancia (km)

P e s o

h (km)

1

20

50

100

200

400

Figura 2: Formas adoptadas por la función de peso según el semi-peso h (parámetro wd de lafunción homogen).

5Aunque los cambios en la varianza no se buscan explícitamente en este paquete.6Gracias a Victor Venema por esta sugerencia.



4

Pero el primer problema que debemos afrontar es que, a menos que la series estén completas, nopodemos calcular sus medias y desviaciones típicas par todo el periodo de estudio. De modo quedebemos comenzar por calcular estos parámetros únicamente a partir de los datos disponibles,usar las series estimadas (tras deshacer la normalización de los datos) para rellenar las lagunasde datos, recalcular las medias y desviaciones típicas, renormalizar los datos, y obtener nuevas

estimas de las series. Este proceso se repite hasta que el cambio máximo de cualquier datorespecto de la iteración anterior sea menor que un determinado umbral (0.005 unidades pordefecto).

2.3. Detección y correción de datos anómalos y saltos bruscos en la media

Después de haber estimado todos los datos, para cada serie original podemos calcular las seriesde anomalías (diferencias entre los datos observados y los estimados), y aplicar sobre ellaspruebas para detectar:

1. Datos anómalos: La serie de anomalías se estandariza, y las anomalías mayores de 5 (pordefecto) desviaciones típicas se borran de los datos originales.

2. Saltos en la media: A la serie de anomalías se le aplica la prueba SNHT Standard Normal

Homogeneity Test , SNHT, por Alexandersson, 1986) en dos etapas:

a) Sobre ventanas de 120 términos que se van moviendo en saltos de 60 términos (va-lores por defecto).

b) Sobre la serie completa.

Los máximos valores de SNHT (llamados tV en este paquete) y sus posiciones en cada serie seguardan en memoria, y las series con los valores más altos, si superan el umbral establecido,se cortan en la posición en que se encontró ese máximo valor de inhomogeneidad, de formaque a partir de esa posición se transfieren todos los valores a una nueva serie (con las mismascoordenadas) y se borran de la original.

Lo ideal sería repetir todo el proceso después de cortar la serie más inhomogénea, puesto queesta inhomogeneidad puede haber influido sobre la valoración de la homogeneidad de las esta-ciones vecinas. Pero esto haría el proceso muy largo si estamos tratado con un elevado númerode estaciones con muchas inhomogeneidades, de modo que se proporciona un factor de toleran-cia para permitir el corte de varias estaciones en cada pasada.

Una vez que todas las inhomogeneidades superiores al umbral se han cortado con la pruebaSNHT aplicada sobre ventanas solapadas, se repite todo el proceso aplicándo esta prueba sobrelas series completas, con lo que se pueden generar mas cortes en las series.

La prueba sobre ventanas móviles se ha implementado para evitar la existencia de múltiplessaltos en la media pueda subestimar los valores del SNHT, mientras que su aplicación a lasseries completas es más sensible y permite detectar saltos más pequeños que en la prueba sobreventanas (con menores tamaños muestrales). De todos modos, el valor por defecto del umbralfijado para la prueba sobre las series completas se ha puesto más alto que en la prueba sobre

ventanas, para evitar que se corten series debido a la presencia de tendencias locales y no de



5

saltos en la media (aunque si las tendencias locales son lo suficientemente fuertes se detectarány tratarán como si fueran saltos).

Después de haber eliminado todas las inhomogeneidades superiores a los umbrales establecidos,se realiza una nueva pasada dedicada únicamente a recalcular todos los datos ausentes (inclu-

yendo los eliminados en la detección de saltos y datos anómalos). Esto se aplica a todas lasseries, tanto si son originales (series no cortadas, o sus primeros fragmentos en caso contrario)como si se trata de las nuevas series creadas tras los cortes efectuados. En este caso, la recons-trucción de las series se efectúa únicamente con los datos de los otros fragmentos, cualquieraque sea el número de datos de referencia fijado. (Salvo cuando no existan datos originales, encuyo caso la estima se realiza según el método general).

3. Aplicación

3.1. Preparación de los datos

Las coordenadas de las estaciones y los datos climatológicos deben suministrarse como se ex-plica en la introducción rápida para que la función de homogeneización pueda leerlos correcta-mente. Otra posibilidad es que el usuario los lea de ficheros estructurados de diferente maneramediante sus propios procedimientos, pudiendo aprovechar las funciones de R para acceder abases de datos relacionales. La única precaución es que los datos deben alojarse en la memoriade R en estos dos objetos:

dat Matriz numérica que contiene los datos, de dimensiones nd, ne (donde nd y ne repre-

sentan el número de datos por estación y el número de estaciones, respectivamente). Losdatos ausentes deben especificarse como NA (el estándar de R).

est.c Tabla de datos con cinco columnas X Y Z Código Nombre, conteniendo las coorde-nadas (X e Y pueden expresarse en in grados o en km, y Z en m), códigos y nombres de lasestaciones. Estas líneas deben disponerse en el mismo orden en que aparecen los datos decada estación en el objeto dat.

3.2. Homogeneización de las series

La función de homogeneización de este paquete se llama homogen, y al llamarla deben sumins-trarse, al menos, estos tres parámetros:

varcli Acrónimo del nombre de la variable climática tratada.

anyi Año inicial del periodo de estudio.

anyf Año final del periodo de estudio.

Estos tres parámetros no tienen asignados valores por defecto, y la función los usará para de-terminar el nombre base de los ficheros de entrada y salida, como se explica en la introducciónrápida. Los demás parámetros (opcionales) que acepta la función son los siguientes:



6

nm Número de datos por año en cada estación (12 por defecto: datos mensuales. Poner nm=1para analizar datos anuales, nm=4 para los estacionales, etc).

nref Número máximo de datos de referencia a emplear para las estimas de los datos. Como seexplica en el apartado de la metodología, todos los datos se estiman como si no existieran

(para calcular las anomalías), como una media ponderada de los datos más próximos7.Este parámetro fija el número máximo de datos a usar en caso de que hubiera muchosdisponibles. (10 por defecto).

dz.max Umbral de tolerancia para los datos anómalos. Por defecto, las anomalías superio-res a 5 desviaciones típicas (de las propias series de anomalías) serán rechazadas (valorconservador).

wd Distancia (en km) a la que los datos pesarán la mitad que los de una estación localizadaen el mismo sitio de la serie a estimar. El valor por defecto es 0 para las dos primerasfases (lo que indica que todos los datos tendrán el mismo peso), y 100 para la última fase

de cálculo final de todos los datos ausentes. Se puede modificar suministrando un vectorde tres valores, como wd=c(0, 200, 50). Cualquier valor adicional será ignorado, y siel vector tuviera menos de tres elementos se repetirá el último valor las veces que seanecesario.

tVt Valor umbral para la prueba SNHT sobre ventanas escalonadas (25 por defecto).

tVf Factor de tolerancia para poder fragmentar varias series en una misma pasada. Por defectovale 0.02, lo que permite un 2% de tolerancia en cada dato de referencia. (Ej.: Si elmáximo valor de la prueba SNHT en una serie vale 30 y se han usado 10 referenciaspara el cálculo de las anomalías, la serie se cortará si el máximo valor de cualquiera de

las series de referencia es menor que 30*(1+0.02*10)=36. (Poner tVf=0 para inhabilitarla fragmentación si cualquiera de las referencias ya ha sido fragmentada en la mismapasada).

swa Tamaño del desfase a aplicar a las ventanas para la aplicación de la prueba SNHT. Elvalor por defecto es 60, lo que significa que la prueba se aplicará a los primeros 2*60términos disponibles, y luego esta ventana de 120 términos se desplazará 60 términoshacia adelante para repetir la prueba, y así sucesivamente hasta alcanzar el final de laserie. Este valor por defecto resulta adecuado para valores mensuales, pero es demasiadogrande para los anuales, y posiblemente demasiado pequeño para datos diarios.

snhtt Valor umbral para la prueba SNHT aplicada a las series completas. Por defecto tieneun valor de 50 (bastante conservador), y puede cambiarse a 0 para inhabilitar esta prueba.

mxdif Máxima diferencia de datos en iteraciones consecutivas. El cálculo iterativo de lasmedias (y, opcionalmente, las desviaciones típicas) de las series se detendrá cuando lamáxima diferencia de cualquier datos respecto a su valor en la iteración anterior sea comomáximo igual a este valor, fijado por defecto en 0.05.

7Nótese que hablamos de los datos más próximos y no de las estaciones más próximas, puesto que la disponi-bilidad de datos irá cambiando probablemente a lo largo del periodo de estudio.



7

force Parámetro lógico para forzar la fragmentación de las series incluso cuando sólo hayauna referencia disponible. Por defecto vale FALSE, evitando las fragmentaciones con unasola referencia.

a Constante a añadir a los datos tras leerlos del fichero de entrada. En combinación con el si-

guiente parámetro b, permite aplicar una transformación lineal a los datos si, por ejemplo,los datos originales vienen expresados en unidades diferentes a las deseadas. (Por defectovale 0).

b Factor a aplicar a los datos. (1 por defecto).

wz Factor a aplicar a las altitudes de la estación antes de calcular la matriz de distancias eu-clídeas. Por defecto vale 0.001, para dar a la coordenada vertical (dada en m) el mismopeso que a las horizontales (que se expresan en km).

deg Ponerlo como TRUE (verdadero) si las coordenadas geográficas se dan en grados, o dejar-

lo en su valor por defecto FALSE (falso) si se dan en km (la unidad de distancia usadainternamente en este paquete).

rtrans Transformación raíz a aplicar a los datos: 2 para raíz cuadrada, 3 para cúbica, etc. (Sepermiten números no enteros; útil si la distribución de frecuencia de la variable se alejade la normal, como suceden con la velocidad del viento, o la precipitación de regionesáridas).

std Tipo de normalización. Por defecto (3), los datos se estandarizarán restándoles la media ydividiéndolos por la desviación típica, pero si la variable estudiada tiene un cero natural(como sucede con la precipitación), puede ser más conveniente establecer std=2 (los

datos se normalizarán únicamente dividiéndolos por la media). Otra posible opción esstd=1, para restarles la media únicamente.

ndec Número de decimales de los datos de salida homogeneizados. (1 por defecto).

mndat Mínimo número de datos para que un fragmento se convierta en una nueva serie. Si sedeja en su valor por defecto (0), se fijará en la mitad del valor del parámetro swa cuando seaplique a datos diarios, y se igualará al valor de nm en caso contrario, con un valor mínimoabsoluto de 5. (Si se da un valor demasiado bajo, las medias y desviaciones típicas de lasseries no serán fiables, y lo mismo sucederá con la reconstrucción de las series).

gp Parámetro gráfico. Darle un valor:0, para no generar ninguna salida gráfica.

1, para obtener únicamente los gráficos descriptivos de los datos de entrada. (No serealizará ninguna homogeneización).

2, para obtener también los gráficos de anomalías.

3 (valor por defecto), para obtener también los gráficos de medias móviles anualesy correcciones aplicadas.

4: como con 3, pero en lugar de medias móviles anuales se representarán sumasmóviles. (Preferible cuando trabajemos con datos de precipitación).



8

leer Dar FALSE si no hay que leer los datos porque ya se ha hecho con otros procedimientosde R.

na.strings Cadena de caracteres que representa los valores ausentes. Su valor por defectoes el estándar de R, "NA" , pero se puede dar cualquier otra cadena, como por ejemplo

na.strings="-999.0" , o incluso un vector de cadenas, como en na.strings=c("-999","-999.0", "-999.0").

nclust Número máximo de estaciones para el análisis de agrupamiento. Por defecto, si elnúmero de series de entrada es mayor que 100, los gráficos iniciales descriptivos de losdatos de entrada se realizarán sobre una muestra aleatoria de nclust series.

maxite Número máximo de iteraciones para el cálculo de las medias de las series. 50 pordefecto, para evitar un tiempo de proceso demasiado largo cuando la convergencia seamuy lenta.

ini Fecha inicial. Vacía por defecto, si se fija (con formato ’AAAA-MM-DD’) se supondrá quelas series contienen datos diarios. (Ver el apartado 7 para una discusión sobre las limita-ciones de la aplicación de la función a este tipo de datos).

vmin Valor mínimo posible (límite inferior) de la variable estudiada. Aunque no tiene nin-gún valor por defecto, se usará vmin=0 si se da el valor 2 al parámetro std. (Ej.: parahomogeneizar precipitaciones o velocidades del viento).

vmax Valor máximo posible (límite superior) de la variable estudiada. Por defecto no se es-tablece ninguno pero, por ejemplo, puede ser útil usar vmax=100 y vmin=0 para datosexpresados como porcentajes, como la humedad o la insolación relativas.

verb Verbosidad. TRUE por defecto, se puede establecer a FALSE para evitar la larga salida detexto en la consola. (En cualquier caso, esa salida se grabará en el fichero de bitácora,como se explica en el apartado siguiente).

Tal como se dice en la introducción rápida, el ejemplo más simple para efectuar una homoge-neización de series con esta función es:

homogen("Tmin", 1956, 2005)

Este ejemplo se puede reproducir si se copian los correspondientes ficheros de datos y es-

taciones en el directorio de trabajo de R. Estos ficheros, llamados Tmin_1956-2005.dat yTmin_1956-2005.est, se pueden encontrar en el archivo comprimido climatol-dat.zip,disponible en http://www.climatol.eu/

Las salidas de este ejemplo se explicarán a continuación.

4. Salidas

La orden de ejemplo homogen("Tmin", 1956, 2005) genera cuatro ficheros de salida, alma-cenados en el directorio de trabajo:



9

Tmin_1956-2005.txt Un fichero de texto que guarda la información del proceso tal comosale por la consola.

Tmin_1956-2005.pdf Un fichero PDF con una colección de gráficos de diagnóstico.

Tmin_1956-2005.dah Un fichero de texto que contiene los datos homogeneizados (con losdatos ausentes rellenados). Tiene la misma estructura que el archivo de entrada Tmin_1956-2005.dat.

Tmin_1956-2005.esh Un fichero de texto con las coordenadas, nombres e informaciónadicional de las estaciones de los datos homogeneizados.

4.1. El fichero *.txt

El fichero de bitácora, en texto claro, comienza por informar de todos los parámetros de lallamada a la función (tanto explícitos como implícitos), para constancia en posibles revisiones

del proceso.Luego sigue el proceso iterativo de relleno de datos ausentes, reflejando la máxima diferenciade los datos al compararlos con la iteración anterior, identificando la estación responsable consu código entre paréntesis. Si se han rechazado datos anómalos durante este proceso, apareceránen líneas como la siguiente:

S63(7) 1966 7: 21.1 -> 14.3 (stan=6.42)

Estas líneas comienzan con el código de la estación y su número de orden (entre paréntesis)en el fichero de entrada. Luego siguen el año y el mes del dato anómalo, su valor, y una flechaseñalándo por qué valor sería sustituido, indicando entre paréntesis el valor de la anomalía

estandarizada). Nótese que el valor indicado por la flecha es sólo una estima aproximada, puestoque el relleno de lagunas de datos definitivo se realizará en la última fase del proceso.

Después del cálculo iterativo de los promedios de las series (y sus desviaciones típicas, si no seha cambiado el valor implícito std=3), se presentan los resultados de las pruebas de detecciónde cambios bruscos en la media de las series. Para cada una de ellas, identificada por su númerode orden, se da el máximo valor tV de la prueba SNHT sobre ventanas escalonadas. Al terminarde analizarlas todas, la (o las) que haya dado el valor más alto será fragmentada en dos partes,y estos cortes quedarán registrados en líneas como, por ejemplo:

M56(10) se corta en 1976 7 (95.1)

Comienzan, al igual que en las líneas de datos anómalos rechazados, con el código y el númeroordinal de la estación, indicando a continuación el año y mes del primer dato después del cortey, entre paréntesis, el valor del test (tV) en ese punto. Desde el término indicado hasta el final dela serie, los datos se borran de la serie original y se trasladan a una nueva serie, con las mismascoordenadas y añadiendo un número ordinal como sufijo del código y el nombre originales dela estación.

Estos bloques de cálculo iterativo de medias (con posible borrado de datos anómalos) y análisisde saltos se repite varias veces según el proceso va pasando por los niveles 1 (pruebas SNHTaplicadas sobre ventanas escalonadas) y 2 (aplicación clásica de SNHT sobre las series comple-tas), y después tiene lugar un nivel 3 final para el cálculo definitivo de los datos ausentes (estavez sin análisis de saltos).



10

El archivo de bitácora termina con los resultados de los cálculos finales de: The log file endswith a set of final computations, including:

ACmx Máximas autocorrelaciones absolutas. La función de autocorrelación de R, acm, se apli-

ca a las series de anomalías, guardando el máximo valor absoluto obtenido para todoslos desfases en cada serie. Elevados valores de autocorrelación pueden indicar falta dealeatoriedad en las anomalías, debiendo revisarse las series correspondientes.

SNHT Valor de la prueba SNHT de las series finales de anomalías, para evaluar la inhomoge-neidad remanente en las mismas.

RMSE Error típico (raíz cuadrada del error cuadrático medio) de los datos estimados. Se calculaa partir de las diferencias entre los datos observados y los calculados, cuando se disponede ambos. Sirve para evaluar los errores que pueden cometerse en el relleno de laguas,y puede ayudar a seleccionar los mejores valores de algunos parámetros de la función

homogen cuando se prueban varios de ellos. Por otra parte, valores elevados del errortípico pueden indicar tanto una mala calidad de la serie original como una singularidaden la ubicación de la estación (que podría estar situada en un lugar con un microclimaespecial).

PD Porcentaje de los datos originales. Cuando una serie se corta en dos o más fragmentos, estevalor sirve para identificar cuál de ellos retiene el mayor número de datos originales.

Primeramente se presentan los resúmenes estadísticos de estas magnitudes, y después se listanlos valores individuales para cada estación (original o derivada del proceso de fragmentación).

4.2. Archivo *.pdf

Otra de las salidas es una serie potencialmente larga (según el valor del parámetro gp) de gráfi-cos de diagnóstico. Las primeras figuras describen los datos de entrada: número total de datosdisponibles en cada paso temporal (figura 3), diagramas de caja (figura 4), y un histograma detodos los datos (figura 5). La inspección de estos gráficos puede revelar la existencia de da-tos muy anómalos u otro tipo de problemas en los datos de entrada, que pueden aconsejar unaacción correctora antes de repetir el proceso de homogeneización.

La siguiente figura es un gráfico de coeficientes de correlación en función de la distancia (figura6). Estos valores de correlación se calculan a partir de las series diferenciadas, para evitar el po-sible impacto de las inhomogeneidades, y se usan todos los pares de observaciones disponibles.Los coeficientes de correlación iguales a 1 o -1 se eliminan previamente, puesto que probable-mente se han originado a partir de sólo dos pares de datos, pero hay que tener en cuenta quealgunos valores de correlación pueden provenir de tres o pocos más. Aunque estos coeficientesno van a tener relevancia para el proceso de homogeneización, este gráfico puede servir paracomprobar que no haya barreras geográficas que provoquen cambios abruptos en las caracte-rísticas climáticas de la zona de estudio. En el ejemplo de la figura 6 se observan tanto valoresaltos como bajos a distancias relativamente pequeñas, indicando el impacto de las diferentescondiciones topográficas de los observatorios en las temperaturas mínimas durante las nochesdespejadas y con viento en calma.



11

1960 1970 1980 1990 2000

0

2

4

6

8

1 0

Nr. de datos de Tmin en todas las estaciones

Años

N r . d e

d a t o s

Figura 3: Número de datos disponibles.

A continuación se realiza un análisis de agrupamiento basado en la matriz de correlaciones, que

da lugar a dos nuevas figuras: un dendrograma, donde pueden verse las estaciones agrupadaspor la similaridad de las variaciones de sus datos, y un mapa de la ubicación de las mismas,identificadas por su número de orden, y con un color distinto según el grupo al que pertenecen.El objeto de este análisis es proveer una primera aproximación a una clasificación climática delas estaciones, aunque el número de grupos, escogido automáticamente por la línea de trazosroja del dendrograma, es probable que no sea el óptimo. Si los grupos son muy diferentes (estánconectados por elevadas distancias en el dendrograma) y su localización geográfica muestraáreas claramente delimitadas, el área de estudio puede incluir discontinuidades climáticas, yel investigador debe considerar la conveniencia de efectuar homogeneizaciones independientespara cada subárea climática.



12

1 2 3 4 5 6 7 8 9 10

− 8

− 6

− 4

− 2

0

2

4

Valores de Tmin (Ene)

Estaciones

V a l o r e s

Figura 4: Ejemplo de diagramas de caja de los datos.

Histograma de todos los datos

Tmin

F r e c u e n c i a

−10 −5 0 5 10 15 20

0

1 0 0

2 0 0

3 0 0

4 0 0

5 0 0

Figura 5: Histograma de todos los datos.





14

20 30 40 50 60 70

1 0

2 0

3 0

4 0

5 0

6 0

Situación de las estaciones de Tmin (2 grupos)

X (km)

Y ( k m )

1

2

34

5

6

7

8

9

10

Figura 8: Mapa de la situación de las estaciones, coloreadas por grupos.

Tras los gráficos descriptivos, encontramos los que describen el análisis de las series de anoma-lías, como en la figura 9, donde las anomalías figuran dibujadas como trazos verticales de colorazul. Cuando el valor máximo de la prueba de saltos en la media supera el umbral establecido,la posición donde se va a cortar la serie se marca con una línea vertical de trazos rojos, rotuladaen su parte superior con el valor de la prueba (redondeado por defecto). En la parte inferior delgráfico se dibuja en verde la distancia al dato más próximo en cada paso temporal, en km (conescala logarítmica).

Todas las series cortadas se muestran en gráficos similares, para permitir inspeccionar el procesode homogeneización de forma subjetiva. Los primeros cortes serán probablemente muy claros(como en la figura 9), mientras que los últimos podrían ser discutibles, especialmente si elumbral de la prueba, tVt, se fijó en un valor relativamente bajo. En este caso puede resultaraconsejable repetir el proceso con un umbral más alto.

Después de los gráficos de anomalías de las series cortadas en la primera fase siguen unosgráficos resumen, que muestran los máximos valores de las pruebas de salto de las series rema-nentes (figura 10, con barras cuyo color varía de verde hacia rojo al aumentar su valor), y un

histograma de frecuencias de los valores de todas las series (figura 11). Ambas figuras muestranla distribución de los máximos valores de las pruebas de salto, lo que permite juzgar si los va-lores más altos corresponden a series con destacadas inhomogeneidades o si más bien se sitúansimplemente en la cola derecha de la distribución de frecuencias de las pruebas de salto en lamedia.

Este bloque de pruebas de salto en la media de las series de anomalías y cortes se repite en lasegunda fase, donde la prueba de SNHT se aplica a las series completas, seguido de los corres-pondientes gráficos de barras e histograma de los máximos valores de la prueba en las seriesresultantes. A continuación aparecen otros dos gráficos referidos al número de cortes sufridopor las estaciones: un histograma del número de cortes por estación (figura 12), y un gráfico debarras indicando el número de cortes por año (figura 13). Una acumulación de muchos cortes



15

en el mismo año puede indicar un cambio en el instrumental o en el método de observación enuna parte significativa de la red8.

1960 1970 1980 1990 2000

− 4

− 2

0

2

4

Tmin at M56(10), Buena Vista

Años

A n o m a l í a s e s t a n d a r i z a d a s ( o b s e r v a c i o n e s − e s t i m a s )

1

10

100

min.d.

(km)

95

Figura 9: Análisis de las anomalías, señalando el punto de corte más significativo.8Estos cambios nunca deben aplicarse simultáneamente a toda una red, puesto que no quedarían series de

observación que sirvieran de referencia para juzgar el efecto de los mismos.



16

5 10 15 20 25

0

2 0

4 0

6 0

8 0

tV máximo por estaciones

Estaciones

t V m á x i m o

Figura 10: Máximos valores de las pruebas de salto tras el proceso de fragmentación. (Algunasseries no muestran ningún valor porque son demasiado cortas para poder aplicar la prueba enventanas escalonadas).

Histograma de los SNHT máximos

SNHT

F r e c u

e n c i a

0 5 10 15 20 25 30

0 . 0

0 . 5

1 . 0

1

. 5

2 . 0

2 . 5

3 . 0



17

Figura 11: Histograma de los máximos valores de las pruebas de salto residuales.

Número de cortes por estación

Número de cortes

N ú m e r o d e e s t a c i o n e s

0 2 4 6 8

0 . 0

0 . 5

1 . 0

1 . 5

2 . 0

2 . 5

3 . 0

Figura 12: Histograma del número de cortes por estación.

1960 1970 1980 1990 2000

0

2

4

6

8

1 0

Número de cortes por año

Años

N ú m e

r o d e c o r t e s



18

Figura 13: Número de cortes por año efectuados en el proceso.

Como se mencionó anteriormente, la tercera fase del proceso de homogeneización se dedicaal relleno de todos los datos ausentes, tanto los que ya faltaban en las series originales como

los derivados del borrado de datos anómalos y del proceso de fragmentación de las series. Estaúltima fase genera otros dos bloques de gráficos: de anomalías, similares a los de las dos fasesanteriores, y de series homogeneizadas y correcciones aplicadas.

La figura 14 muestra un ejemplo de los gráficos de anomalías finales, en los que líneas verticalesa trazos indican la localización de los máximos valores de las pruebas SNHT (en verde paraventanas escalonadas, siempre que haya un mínimo de 2*swa datos, y en negro sobre toda laserie). También se dibuja una recta de tendencia si es significativa al nivel α = 0,05.

Tras los gráficos de anomalías, se presenta un gráfico por cada serie original que muestra, enla parte superior, las medias anuales móviles (o sumas móviles, si se especificó gp=4), y en laparte inferior, las correcciones aplicadas en cada reconstrucción (véase el ejemplo de la figura15).

Las últimas figuras consisten en histogramas de anomalías normalizadas (colorenado de rojolas frecuencias de las que exceden el umbral de corrección de datos anómalos), y de valoresmáximos de las pruebas de salto (tVx y SNHT). Adviértase que éstos pueden ser mayores quesus umbrales si, como en la aplicación por defecto, el peso de las observaciones vecinas esinferior en la última etapa de cálculo de todos los valores ausentes que en las fases anterioresde detección y corrección.

El último gráfico del fichero PDF generado representa los valores SNHT frente a los erroresRMSE (figura 16), de forma que pueda visualizarse de forma rápida la calidad (o singularidad)

de cada una de las series reconstruidas.



19

1960 1970 1980 1990 2000

−

4

− 2

0

2

4

Tmin at S33(5), Pastores

Años

A n o m a l í a s e s t a n d a r i z a d a s ( o b s e r v a c i o n e s − e s t i m a s )

1

10

100

min.d.

(km)

17

25

Figura 14: Anomalías de las series finales, con las localizaciones de los máximos valores deSNHT y, si es significativa, la recta de la tendencia general.

4

5

6

7

8

9

10

Tmin at S11(3), Miraflores

M e d i a s a n u a l e s m ó v i l e s

1960 1970 1980 1990 2000

−6

−4

−2

0

2

Años

T é r m i n o s c o r r e c t o r e s

Figura 15: Serie original (en negro) y reconstrucciones de los valores anuales móviles (arriba),y correcciones aplicadas a cada fragmento (abajo).



20

0.0 0.2 0.4 0.6 0.8 1.0 1.2

0

1 0

2 0

3 0

4 0

5 0

Calidad/singularidad de las estaciones

RMSE

S N H T

1

2

3

4

5

6

7

8

9 10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

Figura 16: Gráfico de los valores de SNHT y RMSE de cada serie final (original ofragmentada).

4.3. Ficheros *.esh y *.dah

Los ficheros *.esh y *.dah son equivalentes a los ficheros de entrada *.est y *.dat, perocontienen los resultados de la homogeneización. Sin embargo, el archivo de estaciones homo-geneizadas *.esh presenta información adicional, como podemos ver en las primeras líneas delfichero Tmin_1956-2005.esh generado en el ejercicio de ejemplo:

27 53.9 456 "S03" "La Perla" 79 1 0 1231.8 26.5 123 "S08" "El Palmeral" 11 2 0 8.449.2 30 154 "S11" "Miraflores" 31 3 0 5.1

En cada línea aparecen los datos siguientes (los cinco primeros son los mismos que en el ficherode entrada Tmin_1956-2005.est):

1 Longitud, X.

2 Latitud, Y.

3 Altitud, Z.

4 Código de la estación, Cód.



21

5 Nombre de la estación, Nombre.

6 Porcentaje de datos originales, PD.

7 Índice de la estación original en el fichero de entrada, io.

8 Marca binaria que indica si la estación estaba funcionando al final del periodo de estudio (1)o no (0), op.

9 Máximo valor SNHT, SNHT.

X y Y se expresan en las mismas unidades (km o grados) que en el fichero de entrada. En cuantoal índice de la estación original (io), su propósito es identificar qué fragmentos pertenecen ala misma serie original. Así, la octava estación de nuestro ejemplo (Esmeraldas), se ha cortadodos veces, y por tanto en el fichero Tmin_1956-2005.esh aparecen tres fragmentos, y se hareconstruido una serie completa para cada uno de ellos (disponible en Tmin_1956-2005.dah):

31.6 56.2 498 "S40" "Esmeraldas" 48 8 0 21.131.6 56.2 498 "S40-2" "Esmeraldas-2" 7 8 0 5.531.6 56.2 498 "S40-3" "Esmeraldas-3" 8 8 0 3.1

Por estas líneas (que no aparecen consecutivas en el fichero) podemos ver que todas pertenecen ala misma serie original, dado que: a) Tienen las mismas coordenadas; b) sus códigos y nombresson iguales, excepto por el sufijo numérico que se ha añadido para diferenciarlas; y c) su índicede estación original io es el mismo (8). Pero hay que tener en cuenta que los sufijos numéricosno tienen porqué seguir el orden cronológico de los fragmentos en la serie original, puesto

que se crean por orden de importancia del salto en la media. En nuestro ejemplo, si buscamoslas palabras S40 y se corta en el fichero de bitácora Tmin_1956-2005.txt, encontramos lassiguientes dos líneas, que indican que el primer corte (que da lugar a la serie S40-2) tiene lugaren marzo de 2000, mientras que el segundo corte tiene lugar en un punto anterior (marzo de1996), aunque cree la serie S40-3:

S40(8) se corta en 2000 3 (47.2)S40(8) se corta en 1996 3 (28.5)

5. Discusión y sugerencias

Si se necesitan con rapidez datos homogeneizados para un proyecto determinado, el investiga-dor puede verse tentado a usar esta función de homogeneización como una caja negra, pero esaconsejable que revise los ficheros de salida para comprobar si los parámetros usados, tantofijados en la llamada a la función como los establecidos por defecto, son apropiados para la redclimática objeto de estudio. Hay que tener en cuenta que los valores óptimos de los parámen-tros variarán según el elemento climático de que se trate, su variabilidad espacial, y la densidadtemporal y espacial de las observaciones, y por consiguiente no pueden proveerse valores deaplicación universal.



22

Es más, los parámetros escogidos pueden ser óptimos o no dependiendo del objeto final delanálisis de las series. Ejemplo: si lo que se desea es obtener normales climáticas, los ajustes dela varianza no tendrán importancia, en tanto que éstos serán cruciales si el objetivo es calcularperiodos de retorno de valores extremos. En este último caso, se puede limitar la disminución devarianza asociada a las estimas ponderadas fijando una distancia de peso pequeña en la tercera

fase, tal como wd=c(0,200,30)), o evitar las ponderaciones totalmente especificando en estafase final de cálculo de todos los valores ausentes sólo se usará una estación de referencia(nref=c(10, 10, 1)).

Por consiguiente, debería prestarse atención a los gráficos de diagnóstico y ver si quedan inho-mogeneidades pendientes de corregir, en cuyo caso habría que bajar los umbrales de tVt y/osnhtt, o si, por el contrario, valores demasiado bajos de estos parámetros han producido unaexcesiva fragmentación de las series. Algunos autores han publicado valores críticos (ej.: Kha-liq and Ouarda, 2007), y en el anexo de esta guía se comentan unas pruebas con series aleatoriasrealizadas expresamente durante el desarrollo de Climatol.

De igual modo, dependiendo de la curtosis de la variable estudiada, puede que se hayan borradodemasiados (o demasiado pocos) datos anómalos. El valor por defecto, 5 desviaciones típicas,es bastante conservador. Puede ajustarse a las necesidades de cada caso, e incluso fijar distintosvalores para cada fase del proceso. Por ejemplo, dz.max=c(6, 3.5, 9) sólo eliminará lasmayores anomalías en la primera fase, y será más drástico en la segunda, para no eliminar yaningún otro dato en la última fase (salvo que su anomalía superara las 9 desviaciones típicas,caso improbable tras las dos primeras fases correctivas).

No olvidar fijar deg=TRUE si las coordenadas de las estaciones están expresadas en grados (com-probando que se da primero la longitud, y luego la latitud), y también escoger el tipo de nor-malización apropiado, prefiriendo std=2 para las variables con un cero natural (como la preci-

pitación o la velocidad del viento) y aplicando una transformación raíz si el histograma de losdatos muestra una distribución con clara forma de L. Nótese que std=1 aplicará correccionesconstantes a los datos, y por tanto no se tendrán en cuenta posibles diferencias estacionales enlas inhomogeneidades, ni se ajustará ningún cambio en la varianza.

Si se va a homogeneizar un pequeño número de series, es recomendable fijar tVf=0 para evitarcortar muchas series a la vez. En estos casos puede darse la situación de que, en algún pasotemporal del periodo de estudio (normalmente al principio, que es cuando suele haber menosobservatorios en funcionamiento), sólo se disponga de datos en una o dos series. En cualquierpaso de tiempo debe de existir al menos un dato válido en alguna de las series para que elprograma pueda funcionar, pero en este caso los datos de ese paso de tiempo de todas lasdemás series se rellenarán tomando ese dato como única referencia, cuya calidad no se podrácomprobar. Si en lugar de un único dato se dispone de dos, y sus anomalías o pruebas de saltoresultan demasiado grandes, el problema será decidir cuál de los dos datos es el incorrecto. Porconsiguiente, cuando sólo existan dos datos en un determinado paso temporal, no se borraráninguno ni se efectuará ningún corte en las series en ese punto, limitándose el programa aadvertir si se han sobrepasado los umbrales correspondientes en el fichero de bitácora, como enlas líneas siguientes:



23

Para datos anó-malos:

... Sólo tiene 1 referencia! (No se elimina)

Para saltos en lamedia:

... podría cortarse en ..., pero sólo tiene una referencia

(Los puntos serán reemplazados por la información sobre la estación y el dato que resultensospechosos).

En estos casos, si la prueba de salto en la media resulta significativa, el único modo de decidircuál de los dos es el sospechoso es basarse en los metadatos. La historia de las estaciones puedecontener información sobre cuál de las dos estaciones sufrió un cambio de emplazamiento u otracircunstancia que pueda haber afectado a las observaciones, en cuyo caso puede procederse acortar manualmente la serie y volver a lanzar el proceso de homogeneización.

Nunca se insistirá suficientemente en lo importante que es guardar registrados todos los cambiosque afecten a un observatorio o sus alrededores, y varios métodos de homogeneización hacenuso de ellos (Aguilar et al., 2003). Desgraciadamente los metadados suelen ser muy incompletos

o incluso no existir en absoluto, y es por esto por lo que este paquete funciona sin ellos, aunquees totalmente recomendable que el usuario los tenga en cuenta para comprobar si los resultadosde la homogeneización son consistentes con ellos.

Otra posibilidad cuando tratamos de homogeneizar unas pocas series y sólo tenemos dos dispo-nibles en algún subperiodo es complementarlas con otras derivadas de productos de reanálisisen puntos de rejilla próximos, aunque estos reanálisis difícilmente se extenderán hacia atrás másallá de mediados del siglo XX.

Resumiendo, resulta recomendable realizar homogeneizaciones con diferentes parámetros ycomprobar cuál de ellas resulta más satisfactoria. Para evitar la reescritura de los archivos de sa-

lida de cada proceso, pueden renombrarse con la función outrename, que añadirá un sufijo a sunombre base. Ejemplo: si queremos conservar las salidas anteriores como Tmin_1956-2005-old.*,usaremos la orden:

outrename("Tmin", 1956, 2005, "old")

Una vez optimizados los parámetros para la aplicación de homogen a una base de datos concreta,pueden conservarse para futuras homogeneizaciones de la misma, como sería el caso cuandoesa base se está actualizando con nuevos datos con el paso del tiempo, siendo recomendableentonces rehomogeneizarla una vez al año, dado que los datos añadidos pueden servir paraconfirmar o rechazar inhomogeneidades localizadas en la parte final de las series.

6. Explotación de las series homogeneizadas

Una vez obtenido un conjunto de series satisfactoriamente homogeneizadas, el investigador eslibre de aplicar sus propios análisis a las mismas y obtener valores estadísticos y gráficos quemuestren la variabilidad espacial y temporal del elemento climático objeto de estudio. Para faci-litar algunos de los cálculos más frecuentes, este paquete incluye la función dahstat, que puedeinvocarse con los siguientes parámetros, de los cuales sólo los tres primeros han de asignarseexplícitamente (los demás adoptan por defecto los valores que aparecen entre paréntesis):



24

varcli Acrónimo del nombre de la variable climática estudiada.

anyi Primer año del periodo de estudio.

anyf Último año del periodo de estudio.

anyip Primer año del periodo de cálculo (anyi).

anyfp Último año del periodo de cálculo (anyf).

nm Número de datos por año en cada estación (12).

ndec Número de decimales de los valores calculados (1).

vala Valor anual a calcular (2). Pueden dársele los valores 0 (no calcular ningún valor anual),1 (sumar los nm valores del año), 2 (calcular su media, que es el cálculo por defecto), 3 (elvalor anual será el máximo de los nm valores), o 4 (el mínimo de los valores).

mnpd Filtrar las series que no posean este porcentaje mínimo de datos originales ( 0).

mnsh Filtrar las series cuyo SNHT sea superior a éste (0).

out Parámetro estadístico a calcular (el nombre del fichero de salida llevará esta extensión):"med" Medias (parámetro por defecto)."mdn" Medianas."max" Máximos."min" Mínimos."std" Desviaciones típicas.

"q" Cuantiles (ver el parámetro prob)."tnd" Tendencias.Cualquier otra opción evitará realizar cálculo alguno, pero leerá las series homogeneiza-das, permitiendo al usuario la aplicación de sus propios análisis.

prob Probabilidad para el cálculo de los cuantiles (si se ha establecido la opción out="q" .0.5 por defecto, lo que produce los mismos resultados que out="mdn" ).

func Poner func=TRUE para filtrar las series que no estuvieran en funcionamiento al final delperiodo de estudio. (FALSE por defecto, con lo que se usarán todas las series).

pernum Número de años sobre los que expresar los valores de las tendencias (100).

eshcol Columnas del fichero de estaciones homogeneizadas (*.esh) a incluir en el ficherode salida (4 por defecto, para identificar los datos calculados sólo por el código de cadaestación).

sep Cadena de caracteres que debe usarse para separa los datos de salida (" " por defecto).

eol Estilo de finalización de línea. (Código de nueva línea, "\n", por defecto).



25

Si no se especifica otra cosa, los cálculos deseados se aplicarán a todas las series, dado que losvalores por defecto de los parámetros mnpd, mnsh y func no filtrarán ninguna serie.

En el fichero de salida, los valores calculados estarán separados por un espacio en blanco, peroeste comportamiento puede cambiarse con el parámetro sep. Por ejemplo, se pueden obtener

valor separados por punto y coma estableciendo sep=’;’. (Recordar que en R las cadenas decaracteres pueden especificarse tanto con dobles comillas como con apóstrofres).

Los ficheros de salida tendrán el mismo nombre base que los demás, y su extensión será lacorrespondiente a la opción out elegida, con la excepción de los cuantiles, cuya extensión seráqPP, siendo PP la probabilidad escogida con el parámetro prob expresada en%.

Por tanto, para obtener las normales mensuales del periodo 1971-2000 a partir de las tempera-turas mínimas previamente homogeneizadas, haríamos:

dahstat("Tmin", 1956, 2005, 1971, 2000)

Pero si lo que queremos es calcular las tendencias para todo el periodo de estudio 1956-2005,expresadas ◦C/década (en lugar de por siglo) con dos decimales, e incluir las coordenadas delas estaciones (columnas 1 y 2 del fichero Tmin_1956-2005.esh) tras los códigos de estación,ordenaríamos:

dahstat("Tmin",1956,2005,out="tnd",pernum=10,vala=1,ndec=2,eshcol=c(4,1,2))9

De este modo obtendríamos la lista de tendencias en un fichero llamado Tmin_1956-2005.tnd,que podría ser importado por un SIG para producir un mapa de tendencias. Alternativamente,podríamos generar ese mapa sin abandonar R, con la ayuda de otros paquetes (Bivand et al.,2008). Ejemplo del comienzo de ese fichero de salida:

"Cód." "X" "Y" "Ene" "Feb" "Mar" "Abr" "May" "Jun" "Jul" "Ago" "Sep" "Oct" "Nov" "Dic" "Anual""S03" 27 53.9 0.04 0.01 0.05 -0.07 0.13 0.19 0.08 0.17 -0.21 0.1 -0.1 0.12 0.53"S08" 31.8 26.5 -0.02 -0.01 0.15 0.06 0.17 0.21 0.05 0.13 -0.15 0.12 -0.04 0.08 0.75"S11" 49.2 30 -0.01 0.05 0.08 0.05 0.24 0.23 0.11 0.12 -0.19 0.1 -0.09 0.07 0.76...

7. Y si los datos son diarios o sub-diarios?

En los últimos años ha aumentado el interés por la homogeneización de datos diarios. Esta es

una tarea bastante difícil, puesto que la detección de saltos en la media de las series es básica-mente un problema de relación señal/ruido, y los valores diarios son generalmente demasiadoruidosos para permitir esa detección. Por tanto, mientras prosiguen los trabajos en busca de lastécnicas apropiadas para abordar este reto, este paquete no debería aplicarse a la detección decambios en las medias en series de datos diarios, a no ser que esos cambios sean suficientementegrandes o la variabilidad de los datos muy pequeña.

Un ejemplo real de baja variabilidad se encontró al investigar un cambio en la media de unaestación termométrica situada en un aeropuerto, gracias a que la cercana ubicación de los sen-sores de temperatura en pista proporcionaron una referencia muy próxima para estudiar la serie

9

La función de concatenación c sirve en R para suministrar un vector de valores.



26

problemática, y eso permitió detectar que una operación de mantenimiento defectuosa produ- jo un cambio en la media de 0.9◦C. La proximidad de los registros de referencia hizo posiblerealizar esa detección en las diferencias de los datos diezminutales, cosa que hubiera sido muyproblemática si esa referencia no hubiera existido, incluso con series de datos diarios.

El principal problema con los datos subdiarios cuando no hay referencias muy próximas resideen la falta de sincronía entre las medidas, puesto que el paso de perturbaciones frontales océlulas convectivas posiblemente tormentosas tendrá lugar a distintas horas en las estaciones dela red de observación. Esto puede pasar incluso con series diarias de precipitación, cuando unaguacero tiene lugar alrededor de la hora de la observación y es asignado a un día o al siguientesegún la hora en la que alcanza a cada pluviómetro.

De todas formas, aunque generalmente sea desaconsejable la detección de saltos en la mediaen series de datos diarios, es en ellas donde debería realizarse el control de calidad, antes queen las series agregadas de datos mensuales. Nótese que un error de 10 ◦C al leer o transcribir latemperatura máxima de un día disminuye a sólo 0.17◦C en la temperatura media mensual (si

se calcula como promedio de las máximas y las mínimas). Por consiguiente, siempre que seaposible la detección y correción de datos anómalos debería realizarse sobre las medidas origina-les, tanto diarias (lo más frecuente en el caso de estaciones de aficionados colaboradores) comoa instervalos más cortos (normalmente en estaciones meteorológicas automáticas), aunque eneste último caso la ya mencionada falta de sincronización puede hacer muy difícil disponer deestaciones de referencia útiles.

Para evitar conflictos en los nombres de los ficheros de datos diarios y mensuales, aquéllos sedistinguen en este paquete añadiendo el sufijo -d al acrónimo de la variable. Por ejemplo, sihemos estado trabajando con los valores mensuales del fichero Tmin_1956-2005.dat, el quecontuviera los valores diarios se llamaría Tmin-d_1956-2005.dat (y el fichero de estaciones

Tmin-d_1956-2005.est sería una mera copia del Tmin_1956-2005.est). La llamada a lafunción homogen sería en este caso:

homogen("Tmin", 1956, 2005, nm=0, tVt=0, ini="1956-01-01")

Los gráficos generados en Tmin-d_1956-2005.pdf pueden revelar algún salto importante enla media que podría valer la pena corregir. En este caso, puede hacerse una nueva aplicaciónde homogen estableciendo valores apropiados para tVt y snhtt) en lugar del valor cero sumi-nistrado antes para evitar el análisis de saltos. También será importante usar una ventana swagrande para la prueba escalonada de SNHT, puesto que la persistencia de un determinado tipode circulación atmosférica puede inducir algún periodo de anomalías diarias altamente autoco-

rrelacionadas, cosa más difícil de ver en las series de datos mensuales.Después de haber obtenido las series de datos diarios con los datos ausentes rellenados y losdatos anómalos corregidos mediante la función homogen, el usuario puede desear obtener lascorrespondientes series mensuales. Para ello puede hacer uso de la función dd2m, que compartemuchos parámetros con dahstat: varcli, anyi, anyf, anyip, anyfp y ndec. Los otros pará-metros de dd2m también resultarán familiares, por ser iguales o similares a otros que ya hemosvisto con anterioridad:

ini Fecha inicial, sin valor por defecto. Debe suministrarse con el formato ’AAAA-MM-DD’(AAAA=Año, MM=mes y DD=día), para permitir una correcta asignación de cada datodiario al mes que le corresponda, puesto que los datos diarios no tienen porqué empezar



27

el 1 de enero.

valm Valor mensual a calcular: 1 (suma), 2 (media, el valor por defecto), 3 (máximo), o 4(mínimo).

nmin Número mínimo de datos diarios disponibles en un mes para calcular el valor mensual(15 por defecto).

na.strings Código de ausencia de dato en el archivo de datos diarios ( "NA" por defecto,el estándar de R).

Por tanto, aplicaríamos esta función a los datos diarios “homogeneizados” (aunque lo más se-guro es que no hayamos corregido saltos en la media) de nuestro ejemplo del siguiente modo:

dd2m("Tmin", 1956, 2005, ini=’1956-01-01’)

El fichero de salida se llamará Tmin-m_1956-2005.dah, con el nuevo sufijo -m que nos indicaque contiene datos mensuales calculados a partir del fichero de datos diarios, evitando así so-breescribir un posible fichero Tmin_1956-2005.dah ya existente. Esta será la única salida dedd2m, pudiendo el usuario incluir estas series manualmente en una base de datos más amplia.Si se necesitan, las coordenadas y nombres de las estaciones se pueden tomar directamente delfichero Tmin-d_1956-2005.est, puesto que los datos mensuales conservarán el mismo ordende estaciones.



28

8. Otras funciones de climatol

Este paquete incluye dos funciones adicionales que no guardan relación con la homogeneiza-ción, sino que sirven para generar gráficos de rosas de los vientos y diagramas de Walter y Lieth.

En los siguientes apartados se pueden ver ejemplos de aplicación de las mismas.

8.1. Rosas de los vientos

La función para su generación se llama rosavent10, y acepta los siguientes parámetros:

frec Tabla de datos con las frecuencias del viento.

fnum Número de referencias circulares a dibujar (4 circunferencias por defecto).

fint Incrementos (en%) de las referencias circulares (5 por defecto).

flab Parámetro para indicar qué circunferencias deben rotularse: 1 (sólo la más externa), 2(todas, el comportamiento por defecto), o cualquier otro valor (no se rotulará ninguna).

ang Ángulo donde situar las etiquetas de las circunferencias (3*pi/16).

col Colores para rellenar los polígonos de frecuencias (rainbow(10,.5,.92,start=.33,end=.2)).

margen Vector de márgenes para el gráfico (para ser pasado a la función par, ver la ayuda deparámetros gráficos de R. Por defecto vale c(0, 0, 4, 0)).

key Fijarlo a FALSE si no se desea la leyenda que aparecería si se dan más de una fila (intervalosde velocidad) de frecuencias.

uni Unidades del viento para encabezar la leyenda (’m/s’).

... Cualquier otro parámetro gráfico que quiera establecerse (como el título de la figura, etc).

Ejemplo: Supongamos que tenemos las siguientes frecuencias en una tabla de datos llamadafrecvto (que podemos haber leído de un fichero o calculado por otros medios):

N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW

0-3 59 48 75 90 71 15 10 11 14 20 22 22 24 15 19 333-6 3 6 29 42 11 3 4 3 9 50 67 28 14 13 15 56-9 1 3 16 17 2 0 0 0 2 16 33 17 6 5 9 2

Entonces, la orden siguiente generaría el gráfico de la figura 17:

rosavent(frecvto, 4, 4, ang=-3*pi/16, main="Rosa anual del viento")

No hay ninguna restricción en cuanto al número de columnas de la tabla de datos, con tal de quela primera de ellas corresponda a las frecuencias del viento de dirección norte. (Esta función notiene en cuenta la cabecera de las columnas).

10Contracción del catalán rosa dels vents.



29

Rosa anual del viento

N

S

EW

4 %

8 %

12 %

16 %

0−33−6

6−9

> 9

m/s

Figura 17: Ejemplo de una rosa de los vientos obtenida con la función rosavent.

8.2. Climogramas de Walter y Lieth

Los diagramas climáticos se han usado desde hace mucho como un medio de sintetizar el climade un lugar, y a los botánicos Bagnouls y Gaussen se les ocurrió trazar la línea de precipitacionesmensuales a una escala doble que la de las temperaturas, para poder distinguir de una manerasencilla los meses húmedos de los secos (según que la línea de las precipitaciones se sitúepor encima o por debajo de la de las temperaturas), y lo llamaron diagrama ombrotérmico

(Bagnouls y Gaussen, 1957). No mucho más tarde, Walter y Lieth mejoraron ese diagramaañadiendo información climática suplementaria, señalando los meses con heladas probables oseguras, y encogiendo la escala de las precipitaciones cuando sobrepasa los 100 mm mensualespara permitir su aplicación a todo el mundo, incluso en las zonas más lluviosas (Walter y Lieth,1960).

La función diagwl nos permite generar este tipode diagrama climático a partir de una tabla dedatos que contenga las medias mensuales de precipitación total y temperaturas máximas diarias,mínimas diarias y mínimas mensuales. Los parámetros que admite son (entre paréntesis, losvalores por defecto):

dat Datos climáticos mensuales para generar el diagrama.



30

est Nombre de la estación climatológica ("" ).

alt Altitud de la estación climatológica (NA).

per Periodo de cálculo de los datos mensuales ("" ).

margen Márgenes del gráfico (c(4,4,5,4)).

mlab Iniciales de los meses para rotular el eje X: "es" en español, "en" en inglés; concualquier otro valor rotularán con números del 1 al 12 ("" ).

pcol Color del trazo de precipitaciones ("#005ac8" ).

tcol Color del trazo de temperaturas ("#e81800" ).

pfcol Color de relleno para heladas probables ("#79e6e8" ).

sfcol Color de relleno para heladas seguras ("#09a0d1" ).

shem Fijar a TRUE si la estación está ubicada en el hemisferio sur (FALSE).

p3line Fijar a TRUE para dibujar una línea suplementaria de precipitaciones a escala triple dela temperatura11 (FALSE).

... Otros parámetros gráficos que se desee establecer.

Como ejemplo, supongamos que ya hemos leído nuestras medias climáticas mensuales comouna tabla de datos llamada datcli, y que son las siguientes:

Ene Feb Mar Abr May Jun Jul Ago Sep Oct Nov Dic

Prec. 97.4 69.3 85.5 71.1 48.9 25.1 8.1 37.2 81.6 144.8 110.6 126.5

Max.t. 15.4 16.1 17.2 19.7 23.9 27.9 31.3 31.4 26.5 22.9 18.2 15.8

Min.t. -0.1 -0.4 1.9 4.9 8.3 11.9 14.8 15.5 13.4 9.7 4.6 2.2

Ab.m.t. -5.1 -7.0 -3.5 -1.7 3.4 8.2 11.6 12.2 9.0 3.0 -1.7 -3.6

Para generar el climograma de la figura 18, llamaríamos a la función de este modo:

diagwl(datcli,est="Estación de ejemplo",alt=100,per="1961-90",mlab="es")

Puede verse el gráfico de las medias mensuales de precipitación y temperatura, anotado con las

medias anuales de ambos elementos (en la parte superior) y las temperaturas máximas diariasmedias del mes más cálido y mínimas diarias medias del mes más frío (en el margen izquierdo).La probabilidad de helada se muestra mediante rectángulos achatados adyacentes al eje de 0◦C.Los meses en que el promedio de temperatura mínima diaria es igual o inferior a cero podemosestar seguros de que habrá heladas, y el rectángulo se rellena (por defecto) con un azul másoscuro que si sólo es igual o inferior a cero la temperatura mínima absoluta del mes, en cuyocaso consideramos que las heladas pueden aparecer o no. La trama de líneas azules verticalesindica los meses húmedos, mientras que la trama de puntos rojos señala los áridos, pudiendohacernos una idea de la intensidad de la aridez o el exceso hídrico apreciando el área cubiertapor cada tipo de trama.

11Sugerencia de Bogdan Rosca



31

0

10

20

30

40

50

0

20

40

60

80

100

300

°C mm

Estación de ejemplo (100 m)

1961−90 14.7°C 906 mm

31.4

−0.4

E F M A M J J A S O N D

Figura 18: Ejemplo de un diagrama de Walter y Lieth obtenido con la función diagwl.

9. Bibliografía

Aguilar E, Auer I, Brunet M, Peterson TC, Wieringa J (2003): Guidelines on climate metadata

and homogenization. WCDMP-No. 53, WMO-TD No. 1186. World Meteorological Organiza-tion, Geneve.

Alexandersson H (1986): A homogeneity test applied to precipitation data. Jour. of Climatol.,6:661-675.

Bagnouls F, Gaussen H (1957): Les climats biologiques et leurs classifications. Ann. de Geogr.,355:193-220.

Bivand RS, Pebesma EJ, Gómez-Rubio V (2008): Applied Spatial Data Analysis with R. Sprin-ger, 376 pp.

Daget J (1979): Les modèles mathematiques en écologie. Collection d’Écologie 8, 172 pp, Mas-son, Paris.

Khaliq MN, Ouarda TBMJ (2007): On the critical values of the standard normal homogeneitytest (SNHT). Int. J. Climatol., 27:681687.

Paulhus JLH, Kohler MA (1952): Interpolation of missing precipitation records. Month. Weath.

Rev., 80:129-133.



32

Peterson TC, Easterling DR, Karl TR, Groisman P, Nicholls N, Plummer N, Torok S, Auer I,Böhm R, Gullett D, Vincent L, Heino R, Tuomenvirta H, Mestre O, Szentimrey T, Salinger J,Førland E, Hanssen-Bauer I, Alexandersson H, Jones P, Parker D (1998): Homogeneity Adjust-ments of ’In Situ’ Atmospheric Climate Data: A Review. Int. J. Climatol., 18:1493-1518.

Sokal RR, Rohlf PJ (1969): Introduction to Biostatistics. 2nd

edition, 363 pp, W.H. Freeman,New York.

Walter H, Lieth H (1960): Klimadiagramm Weltatlas. G. Fischer, Jena.

10. Anexo: Valores umbrales para las pruebas SNHT

Si bien ya se han publicado valores críticos para la prueba SNHT de Alexandersson, se realiza-ron simulaciones tipo Monte Carlo adaptadas específicamente a la manera en que se aplica estaprueba en el paquete climatol:

Se generaron 2000 series de ruido blanco de 600 términos con la función de R rnorm parasimular series de 50 años de anomalías mensuales de una estación homogénea con series dereferencia también homogéneas. A cada una de estas series se aplicaron saltos de 0,0 (ningún

salto), 0,5, 1,0, 1,5 y 2,0 desviaciones típicas justo en mitad de la serie (a partir del término301), y se realizó la prueba SNHT en ventanas de 2*swa términos escalonados en pasos de swa= 6, 12, 24, 48, 60, 90, 120, 180, 240 y 300 términos. De cada una de las pruebas se guardó elmáximo valor del estadístico T y la posición donde se encontró. Por tanto, el número total deresultados obtenidos fue: 2000 series * 5 saltos * 10 tamaños de semiventana = 100000.

En primer lugar analizaremos los resultados de las series homogéneas (aquéllas en las que no seintrodujo salto alguno). Estudiando la cola derecha de la distribución acumulada empírica de losvalores máximos del estadístico T de la prueba SNHT, podemos obtener los valores umbralespara evitar falsas detecciones de saltos en la media con niveles de confianza del 90%, 95%,99%, 99.5% y 99.9%. La figura 19 muestra esos umbrales, para estos niveles de confianza y

para los 10 desfases de swa términos de una ventana de tamaño 2*swa. La irregularidad de losgráficos debe atribuirse al azar, pero el aspecto general no debe diferir mucho si se realizara unnúmero mucho mayor de simulaciones. Resulta curioso el máximo que presentan los gráficoscuando se usan ventanas de tamaño medio, puesto que los valores críticos publicados hastaahora muestran un incremento constante (aunque asintótico) al aumentar el tamaño muestral.Pero aquí la prueba sólo se aplica una vez para cada serie cuando la ventana la contiene porcompleto (cuando swa=300), mientras que con ventanas pequeñas la prueba se aplica variasveces sobre la misma serie, permitiendo al estadístico T alcanzar valores más altos.



33

1 0

1 5

2 0

tV's para ruido blanco

swa (muestras de 2*swa términos)

t V

6 12 24 48 60 90 120 180 240 300

Prob.

0.9990.995

0.99

0.95

0.9

6

12

24

48

60

90

120

180

240

300

Figura 19: Valores umbrales (tV) de diez pruebas SNHT aplicadas a ventanas de 2*swa términosescalonadas swa términos hacia adelante sobre series de ruido blanco, para cinco probabilidadesde evitar falsas detecciones de saltos en la media.

Lo siguiente que queremos saber es qué tan buenos son estos valores umbrales tV a la hora dedetectar correctamente los saltos en las medias. Para averiguarlo se contaron los valores de tVsuperiores a los umbrales, y se calificaron como correctos si la localización del salto tenía unerror inferior a 12 términos12 y como erróneos en caso contrario. Ambos se contabilizaron porseparado para cada uno de los 10 valores de swa, las 4 magnitudes de salto y los 5 niveles deconfianza.

La figura 20 resume los resultados de las proporciones de aciertos, mostrando que los saltos

de 0,5 desviaciones típicas son bastante difíciles de detectar, incluso con los mayores tamañosmuestrales, para los que el índice de aciertos es de alrededor del 63% para los cinco niveles deconfianza de evitar falsas detecciones. Los saltos de 1 desviación típica se detectan con mayorfiabilidad: el 95% cuando la prueba se aplica a las series completas (swa=300), y más del 90%incluso para muestras de 120 términos (swa=60), con tal de tolerar una probabilidad de falsospositivos del 10% (nivel de confianza de 0,90). Los saltos mayores (de 1,5 y 2 desviacionestípicas) se detectan casi totalmente con ventanas escalonadas de unos 100 términos o más (desdeswa=48 en adelante).

En cuanto a los falsos saltos (figura 21), con los mayores tamaños muestrales la probabilidadde detectar los de 0,5 desviaciones típicas en una posición errónea alcanza hasta un 35%, mien-

12Un año, si tratamos con series mensuales



34

tras que con las muestras más pequeñas la probabilidad cae a menos del 1% si el umbral dedetección (tVt) se establece suficientemente alto (niveles de confianza superiores a 0,99). Enlas simulaciones de 1 desviación típica esta probabilidad de localización errónea es de un 5 a6% en la mayoría de casos, y para magnitudes de salto mayores es prácticamente desprecia-ble (excepto cuando los tamaños muestrales más pequeños se combinan con bajos niveles de

confianza).

0

1 0 0

2 0 0

3 0 0

4 0 0

5 0 0

6 0 0

0.5 std. deviation shifts

Confidence level of no false breaks

H i t

r a t e

( p e r t h o u s a n d )

0.9 0.95 0.99 0.995 0.999

0

2 0 0

4 0 0

6 0 0

8 0 0



H i t

r a t e


0.9 0.95 0.99 0.995 0.999

0

2 0 0

4 0 0

6 0 0

8 0 0

1 0 0 0



H i t r a t e


0.9 0.95 0.99 0.995 0.999

0

2 0 0

4 0 0

6 0 0

8 0 0

1 0 0 0



H i t r a t e


0.9 0.95 0.99 0.995 0.999

6

12

24

48

60

90

120

180

240

300

Figura 20: Proporciones de acierto para diferentes magnitudes de salto, niveles de confianza deevitar falsas detecciones, y tamaños muestrales de 2*swa términos escalonados.

Por consiguiente, las mayores dificultades se encuentrarán cuando tratemos de detectar saltos





36

0

1 0 0

2 0 0

3 0 0



F a l s e

h i t r a t e


0.9 0.95 0.99 0.995 0.999

0

2 0

4 0

6 0

8 0



F a l s e

h i t r a t e


0.9 0.95 0.99 0.995 0.999

0

2 0

4 0

6 0

8 0



F a l s e

h i t r a t e


0.9 0.95 0.99 0.995 0.999

0

2

0

4 0

6 0

8 0



F a l s e

h i t r a t e


0.9 0.95 0.99 0.995 0.999

6

12

24

48

60

90

120

180

240

300

Figura 21: Proporción de falsa detección para diferentes magnitudes de salto, niveles de confian-za para evitar cortes falsos, y tamaños de ventana de 2*swa términos (escalonadas swa términos,y mostradas con trazos diferentes).

El valor por defecto de swa en la función homogen se ha fijado en 60 (línea azul claro en lasfiguras), como un compromiso entre una buena probabilidad de detección de saltos y un altopoder de discriminación cuando en la serie se presentan más de un salto en la media (situaciónbastante frecuente). Y en la segunda fase del proceso de homogeneización, la aplicación dela prueba SNHT a las series completas permitirá detectar los saltos menores que no se hayancorregido con las pruebas sobre ventanas escalonadas.

Importante nota final: Estos umbrales de tVt se han obtenido a partir de series sintéticas de



37

ruido blanco, pero en el mundo real, las series de anomalías mostrarán inevitablemente algúngrado de autocorrelación y tendencias generales o locales, dependiendo de la variable climática,su variabilidad espacial, la densidad de la red de observación, y el tipo de dato (anual, estacional,mensual, diario, ...). Es por ello por lo que los valores por defecto de tVt y snhtt se hanfijado en la función homogen notablemente más altos que los obtenidos en las simulaciones de

Monte Carlo, y por tanto es aconsejable ajustarlos empíricamente, con ayuda de los gráficos dediagnóstico de una primera aplicación exploratoria, para adaptar los resultados a las necesidadesde cada caso particular.

0 . 0

0 . 2

0 . 4

0 . 6

0 . 8

1 . 0

Goodness index for 0.5 sd shifts


G o o d n e s s i n d e x

0.9 0.95 0.99 0.995 0.999

swa

6

12

24

48

60

90

120

180

240

300

0 . 0

0 . 2

0 . 4

0 . 6

0 . 8

1 . 0

Goodness index for 1.0 sd shifts


G o o d n e s s i n d e x

0.9 0.95 0.99 0.995 0.999

swa

6

12

24

48

60

90

120

180

240

300

Figura 22: Índice de bondad para saltos de 0 5 (arriba) y 1 0 (abajo) desviaciones típicas

Documents

Guia Para Usar Programas r y Clima