Handbook for hydrological data time series management

Librerías usadas

>>> import pandas as pd 

>>> import numpy as np 

>>> from matplotlib import pyplot as plt

Lectura de datos

Los datos utilizados para el manejo de datos pluviométricos fueron descargados del visor de estaciones de la Autoridad Nacional del Agua. El proceso de lectura se realizará para formato .csv o formato .txt. En este caso utilizará un formato .txt. Dado que el archivo cuenta con datos de días y las horas en columnas separadas, se sobreescribirá la columna de fecha para unir ambas variables. Por último, esta se convertirá a formato datetime.

Code 

>>> data = pd.read_csv("precip.txt", sep="\t")

>>> data.FECHA = data.FECHA + " " + data.HORA

>>> data.FECHA = pd.to_datetime(data.FECHA, dayfirst=True)

>>> data.head()

Output

                FECHA   HORA  VALOR
0 1971-04-03 07:00:00  07:00    0.0
1 1971-04-03 19:00:00  19:00    0.0
2 1971-04-04 07:00:00  07:00    0.0
3 1971-04-04 19:00:00  19:00    0.0
4 1971-04-05 07:00:00  07:00    0.0

Exploración de datos

Cuando se tiene un conjunto de datos, se debe realizar una exploración para identificar posibles valores atípicos o cierto comportamiento en la información. La función .describe() permite obtener la  principal información estadísticas de una variable numérica; y para este caso es útil en la variable VALOR, que representa la precipitación.

Code 

>>> data.describe()

Output

              VALOR

count  32035.000000

mean       0.253164

std        3.190159

min        0.000000

25%        0.000000

50%        0.000000

75%        0.000000

max      165.100000

Una de las ventajas de haber convertido la variables data.FECHA a datatime es que se puede hacer operaciones con ellas, por ejemplo, identificar la fecha mínima, máxima o la diferencia entre ellas. 

Code & Output

>>> print ('Primera fecha: ', data.FECHA.min())

Primera fecha:  1971-04-03 07:00:00

>>> print('Última fecha: ', data.FECHA.max())

Última fecha:  2017-12-31 19:00:00

>>> period = data.FECHA.max() - data.FECHA.min()

>>> print('Período: ', period)

Período:  17074 days 12:00:00

La exploración de datos también involucra verificar la existencia de datos faltantes, para ello se ejecuta lo siguiente.

Code & Output

>>> print('Cantidad de datos faltantes : ', sum(data.VALOR.isnull()))

Cantidad de datos faltantes :  0

Evidentemente no hay datos faltantes; sin embargo, puede ser que en los datos haya saltos temporales que no se muestran en la data procesada por ANA. Entonces para ello se va a comparar los días presentes en data.FECHA con la cantidad de valores que deberían haber si el tiempo estuviese completo. Esta operación requiere conocer la fecha mínima, máxima, que fueron calculadas anteriormente, y el periodo, que en este caso son 12 horas. 

Code

>>> range_data = pd.date_range(str(data.FECHA.min()), str(data.FECHA.max()), freq='12H')

>>> range_data

Output

DatetimeIndex(['1971-04-03 07:00:00', '1971-04-03 19:00:00',

               '1971-04-04 07:00:00', '1971-04-04 19:00:00',

               '1971-04-05 07:00:00', '1971-04-05 19:00:00',

               '1971-04-06 07:00:00', '1971-04-06 19:00:00',

               '1971-04-07 07:00:00', '1971-04-07 19:00:00',

               ...

               '2017-12-27 07:00:00', '2017-12-27 19:00:00',

               '2017-12-28 07:00:00', '2017-12-28 19:00:00',

               '2017-12-29 07:00:00', '2017-12-29 19:00:00',

               '2017-12-30 07:00:00', '2017-12-30 19:00:00',

               '2017-12-31 07:00:00', '2017-12-31 19:00:00'],

              dtype='datetime64[ns]', length=34150, freq='12H')

Code & Output

>>> print('Debería haber ', len(range_data),' datos, pero hay', len(data), '\nHay ', len(range_data)-len(data), ' datos faltantes.')

Debería haber  34150  datos, pero hay 32035

Hay  2115  datos faltantes.

>>> print('Porcentaje de completitud de datos : ', round(len(data)/len(range_data)*100,1), '%')

Porcentaje de completitud de datos :  93.8 %

Entonces, se evidencia que aunque no había datos de precipitación faltantes en la data, el parámetro FECHA no contiene la información completa. Se observa un 93.8% de completitud de datos de fecha, mismo que debería someterse a diferentes pruebas para garantizar que sea fiable procesar la información. Ahora la pregunta es, ¿cómo se pueden conocer los datos de fecha faltantes? El proceso se realizará verificando si la información contenida en range_data está contenida en data.FECHA. En el caso que no, se almacenará en el array m_dates que representa las fechas faltantes (missing dates).

Code

>>> m_dates = np.array([])

>>> for i in range (len(range_data)):

...     if sum(range_data[i]==data.FECHA) == 0:

...         m_dates = np.append(m_dates, range_data[i])

... 

>>> print('Datos faltantes:\n', m_dates)

Output

Datos faltantes:

 [Timestamp('1975-06-01 07:00:00', freq='12H')

 Timestamp('1975-06-01 19:00:00', freq='12H')

 Timestamp('1975-06-02 07:00:00', freq='12H') ...

 Timestamp('2017-06-30 07:00:00', freq='12H')

 Timestamp('2017-06-30 19:00:00', freq='12H')

 Timestamp('2017-07-01 07:00:00', freq='12H')]

Para un mejor manejo de los datos puede crearse un data frame que almacene la información de los datos faltantes, para ello se propone el siguiente código, similar al anterior. 

Code

>>> m_Value = np.ones(len(range_data))

>>> for i in range (len(range_data)):

...     if sum(range_data[i]==data.FECHA) == 0:

...         m_Value[i] = 0

...

>>> m_data = pd.DataFrame()    

>>> m_data['Date'] = range_data

>>> m_data['Value'] = m_Value  

>>> m_data = m_data.set_index('Date')

>>> m_data.head()

Output

                     Value

Date

1971-04-03 07:00:00    1.0

1971-04-03 19:00:00    1.0

1971-04-04 07:00:00    1.0

1971-04-04 19:00:00    1.0

1971-04-05 07:00:00    1.0

Para este caso se ha asignado el valor de 1 a la fecha que sí se encuentra en la data original, y 0 a los que no se encuentran. Esto permitirá poder visualizar los datos en gráficos como el que se presenta a continuación. 

Code

>>> plt.subplot(211)

>>> m_bar = m_data.groupby(['Value']).Value.count() 

>>> m_bar.plot(kind='bar', color = ['red', 'green'])

>>> plt.xticks([0,1],['Missing', 'Present'])

>>> plt.title('Barplot Missing Dates')

>>> plt.ylabel('Amount of data')

>>>

>>> plt.subplot(212)

>>> plt.plot(m_data, color = 'black', marker = '|', linewidth = 0)

>>> plt.yticks([0,1],['Missing', 'Present'], rotation = 90)

>>> plt.title ("Time series - Missing data")

>>> plt.xlabel('Dates')

>>>

>>> plt.tight_layout()

>>> plt.show()

Output

Figura 1. Gráficos de datos faltantes.

Mediante gráficas se puede presentar información de manera resumida, y en este caso permiten observar la cantidad de datos faltantes y su disposición a lo largo del tiempo. Después de haber realizado la exploración de datos, debería someterse la data a una limpieza (data cleaning). Para esta oportunidad se va a omitir este proceso, asumiendo que los datos son representativos. 

Análisis de datos

Escalas temporales

Se puede analizar a diferentes escalas temporales según el tiempo de datos que tengamos. Para este caso se tiene una escala de dos datos por día, por lo que se podría analizar a escala diaria, mensual o anual, preferiblemente.

Primero se debe indexar los datos de FECHA para realizar el procesamiento sin errores. 

Code

>>> data = data.set_index('FECHA')

Posteriormente se realizará el remuestreo (resample) de los datos, el cual se maneja con la siguiente estructura en Python. 

file = data_frame.reample("Frecuencia").Operación()

• Frecuencia: es la frecuencia a la cual se van a reordenar los datos, por ejemplo "d" es diario, "m" es mensual, y "y" es anual. Para establecer frecuencias diferentes, por ejemplo cada cinco meses, se puede especificar antes del código, "5m". 
• Operación: aquí se establece el tipo de operación correspondiente al grupo remuestreado, por ejemplo, para precipitaciones se utiliza la suma .sum(). También se puede obtener promedios como .mean(), máximos .max() y mínimos .min().  

Escala diaria

Code & Output

>>> daily_data = data.resample("d").sum()

>>> daily_data.head()

            VALOR

FECHA

1971-04-03    0.0

1971-04-04    0.0

1971-04-05    0.0

1971-04-06    1.4

1971-04-07    0.0

Escala mensual

Code & Output

>>> monthly_data = data.resample("m").sum()

>>> monthly_data.head()

            VALOR

FECHA

1971-04-30    3.3

1971-05-31    1.4

1971-06-30    0.5

1971-07-31    0.0

1971-08-31    0.0

Escala anual

Code & Output

>>> annual_data = data.resample("y").sum()

>>> annual_data.head()

            VALOR

FECHA

1971-12-31    7.0

1972-12-31  185.7

1973-12-31  137.8

1974-12-31   12.2

1975-12-31   37.7

Visualización de las series de tiempo

Plot

El ploteo de datos se realizará para el remuestreo anual. 

Code

>>> plt.plot(annual_data, "go-", label = "Rainfall")

>>> plt.title('Annual rainfall data')

>>> annual_data_mean = np.mean(annual_data.VALOR)

>>> annual_data_median = np.median(annual_data.VALOR)

>>> plt.axhline(annual_data_mean, color = "red", label= f"Mean {round(annual_data_mean,2)} mm")

>>> plt.axhline(annual_data_median, color = "blue", label = f"Median {round(annual_data_median,2)} mm")

>>> plt.xlabel("Year")

>>> plt.ylabel("Rainfall [mm/year]")

>>> plt.legend(), plt.grid(), plt.show()

Output

Figura 2. Gráfico de precipitación anual desde 1971 hasta 2017.

Diagrama de cajas

Los diagramas de caja (boxplot) permiten conocer la distribución de los datos, así como la identificación de datos atípicos. Estos pueden utilizarse a diferentes a escalas temporales, aunque dependiendo la cantidad de datos se puede observar un diagrama de cajas por representativo. En este caso, por ejemplo cuando es a escala mensual, casi todos los datos se encuentran fuera de la caja, mientras que en escala anual, los datos se pueden visualizar y analizar de mejor manera. 

Code

>>> plt.subplot(121)

>>> plt.boxplot(monthly_data.VALOR)

>>> plt.xticks([1], ["Data"])

>>> plt.ylabel('Monthly rainfall [mm/month]')

>>> plt.title('Monthly rainfall data')

>>> plt.subplot(122)

>>> plt.boxplot(annual_data.VALOR)

>>> plt.xticks([1], ["Data"])

>>> plt.ylabel('Annual rainfall [mm/year]')

>>> plt.title('Annual rainfall data')

>>>

>>> plt.tight_layout()

>>> plt.show()

Output

Figura 3. Diagrama de cajas de precipitación mensual y anual.