Programa para realizar el control de datos de precipitación (3 observaciones)¶
| Fecha de Actualización | Descripción de la Actualización | Responsables | |------------------------|---------------------------------------------------------------|----------------------| | 2025-05-13 | Propuesta de control de datos de pricipitación | Bruce Tumbaco | | 2025-05-13 | Revision de la propuesta de control de datos de precipitación | Sandra Torres | | |
El siguiente programa sirve para realizar el control de datos de precipitación Porpuesta realizada en el proyecto ENANDES+
Dirección de estudios, investigación y desarrollo hidrometeorologico.
Copyright. INAMHI <www.inamhi.gob.ec>. Todos los derechos reservados.
Version 1.0¶
Paso 1: Se debe conectar a los datos del Banco Nacional de información Hidrometeorológico¶
In [1]:
from banadih import con_db
In [2]:
# Definir la consulta SQL para obtener los datos de precipitacion
query_precipitacion = """
SELECT *
FROM convencionales."_017140801h" h
WHERE h.id_estacion = 63759
AND h.fecha_toma_dato BETWEEN '1950-01-01' AND '2025-12-31';
"""
# Ejecutar la consulta y almacenar los resultados en un DataFrame
data= con_db.query_to_dataframe(query_precipitacion)
# Mostrar los primeros registros del DataFrame
print(data)
id id_estacion id_usuario fecha_toma_dato \
0 49306142 63759 999 2008-04-30 12:00:00
1 49285777 63759 999 1979-11-01 12:00:00
2 49285778 63759 999 1979-11-02 12:00:00
3 49285779 63759 999 1979-11-03 12:00:00
4 49285780 63759 999 1979-11-04 12:00:00
... ... ... ... ...
41914 49378140 63759 999 2011-08-25 00:00:00
41915 49378142 63759 999 2011-08-26 00:00:00
41916 49378144 63759 999 2011-08-27 00:00:00
41917 49378145 63759 999 2011-08-28 00:00:00
41918 49378147 63759 999 2011-08-29 00:00:00
fecha_ingreso valor validado id_estado_proceso acumulada \
0 2025-07-21 14:51:05 0.3 False 32 False
1 2025-07-21 14:51:05 0.0 False 32 False
2 2025-07-21 14:51:05 0.0 False 32 False
3 2025-07-21 14:51:05 0.0 False 32 False
4 2025-07-21 14:51:05 0.0 False 32 False
... ... ... ... ... ...
41914 2025-07-21 14:51:05 0.0 False 32 False
41915 2025-07-21 14:51:05 0.0 False 32 False
41916 2025-07-21 14:51:05 -888.8 False 32 False
41917 2025-07-21 14:51:05 0.1 False 32 False
41918 2025-07-21 14:51:05 0.1 False 32 False
id_estado_migracion id_origen_dato
0 None 6
1 None 6
2 None 6
3 None 6
4 None 6
... ... ...
41914 None 6
41915 None 6
41916 None 6
41917 None 6
41918 None 6
[41919 rows x 11 columns]
In [3]:
# Información necesaria de la estación
Nombre= "Estación Querochaca- Código: M0258"
Codigo="M0258"
Longitud= -78.605539
Latitud= -1.3671
Altitud= "2865 m s.n.m"
Demarcacion= "Pastaza"
Provincia= "Tungurahua"
Canton= "Cevallos"
In [4]:
pip install folium
Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: folium in /home/jupyter-btumbaco/.local/lib/python3.10/site-packages (0.19.5) Requirement already satisfied: branca>=0.6.0 in /home/jupyter-btumbaco/.local/lib/python3.10/site-packages (from folium) (0.8.1) Requirement already satisfied: jinja2>=2.9 in /opt/tljh/user/lib/python3.10/site-packages (from folium) (3.1.4) Requirement already satisfied: numpy in /home/jupyter-btumbaco/.local/lib/python3.10/site-packages (from folium) (1.26.4) Requirement already satisfied: requests in /opt/tljh/user/lib/python3.10/site-packages (from folium) (2.32.3) Requirement already satisfied: xyzservices in /opt/tljh/user/lib/python3.10/site-packages (from folium) (2024.9.0) Requirement already satisfied: MarkupSafe>=2.0 in /opt/tljh/user/lib/python3.10/site-packages (from jinja2>=2.9->folium) (2.1.5) Requirement already satisfied: charset-normalizer<4,>=2 in /opt/tljh/user/lib/python3.10/site-packages (from requests->folium) (3.1.0) Requirement already satisfied: idna<4,>=2.5 in /opt/tljh/user/lib/python3.10/site-packages (from requests->folium) (3.4) Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/tljh/user/lib/python3.10/site-packages (from requests->folium) (1.26.15) Requirement already satisfied: certifi>=2017.4.17 in /opt/tljh/user/lib/python3.10/site-packages (from requests->folium) (2022.12.7) [notice] A new release of pip is available: 24.3.1 -> 25.2 [notice] To update, run: pip install --upgrade pip Note: you may need to restart the kernel to use updated packages.
In [5]:
import folium
from IPython.display import display
from folium.plugins import FloatImage
# Datos de la estación meteorológica
nombre_estacion = Nombre
latitud = Latitud
longitud = Longitud
# Crear el mapa centrado en la estación
mapa = folium.Map(location=[latitud, longitud], zoom_start=12)
# Agregar un marcador en la estación meteorológica
folium.Marker(
location=[latitud, longitud],
popup=f"{nombre_estacion}\nLongitud: {longitud} Latitud: {latitud}",
tooltip="Haz clic para más info",
icon=folium.Icon(color="red", icon="cloud")
).add_to(mapa)
# Agregar un título al mapa (superior)
titulo_html = f"""
<div style="position: fixed;
top: 10px; left: 50px; width: 400px;
background-color: white; z-index:9999; font-size:16px;
padding: 10px; border-radius: 5px; box-shadow: 2px 2px 5px rgba(0,0,0,0.3);">
<b>{nombre_estacion}</b><br>
Longitud: {longitud} | Latitud: {latitud}
</div>
"""
mapa.get_root().html.add_child(folium.Element(titulo_html))
# Mostrar el mapa en Jupyter Notebook
display(mapa)
Make this Notebook Trusted to load map: File -> Trust Notebook
Paso 2: Se debe tomar los datos necesarios "fecha y el valor"¶
In [6]:
# Seleccionar solo las columnas necesarias
data_filtrada = data[['fecha_toma_dato', 'valor']]
# Mostrar las primeras filas para verificar
print(data_filtrada.head())
fecha_toma_dato valor 0 2008-04-30 12:00:00 0.3 1 1979-11-01 12:00:00 0.0 2 1979-11-02 12:00:00 0.0 3 1979-11-03 12:00:00 0.0 4 1979-11-04 12:00:00 0.0
Paso 3: Se grafica un boxplot para observar la distribucion de los valores de precipitación¶
In [7]:
import pandas as pd
import matplotlib.pyplot as plt
# Asegurar que los valores sean tipo float64 (evitar Decimals)
data_filtrada['valor'] = pd.to_numeric(data_filtrada['valor'], errors='coerce')
# Eliminar valores nulos antes de graficar
data_filtrada_limpia = data_filtrada.dropna(subset=['valor'])
# Crear figura y boxplot
plt.figure(figsize=(8, 6))
plt.boxplot(data_filtrada_limpia['valor'].astype(float), vert=True)
plt.ylabel("Precipitación (mm)", fontsize=16)
plt.grid(True)
# Eliminar etiquetas del eje X
plt.xticks([])
# Crear título multilínea centrado
titulo = (
f"{Nombre}\n"
f"Demarcación: {Demarcacion}, Provincia: {Provincia}, Cantón: {Canton}\n"
f"Latitud: {Latitud}°, Longitud: {Longitud}°, Altitud: {Altitud}"
)
# Posición del título (modificable con x e y)
plt.suptitle(titulo, x=0.54, y=0.96, fontsize=16) # Tamaño del título ajustado
# Ajustar tamaño de los ticks del eje Y
plt.tick_params(axis='y', labelsize=16)
plt.savefig(f"boxplot_precipitacion_{Codigo}.png", dpi=200)
plt.tight_layout()
plt.show()
/tmp/ipykernel_175955/3674985530.py:5: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy data_filtrada['valor'] = pd.to_numeric(data_filtrada['valor'], errors='coerce')
In [8]:
# Filtrar los valores menores a -800 mm en la columna 'valor'
valores_erroneos = data_filtrada[data_filtrada['valor'] < -800]
# Contar cuántos son
cantidad_erroneos = valores_erroneos.shape[0]
# Imprimir el resultado
print(f"Cantidad de valores de precipitación menores a -800 mm: {cantidad_erroneos}")
Cantidad de valores de precipitación menores a -800 mm: 675
Paso 4: Se reemplaza los valores -888.8 por cero (trazas)¶
In [9]:
# Reemplazar los valores menores a -800 mm por 0 en la columna 'valor'
data_filtrada.loc[data_filtrada['valor'] < -800, 'valor'] = 0
Paso 5: Se establece un rango desde 0 hasta 300 para precipitación¶
In [10]:
import numpy as np
# Reemplazar los valores fuera del rango 0–250 mm por NaN
data_filtrada.loc[(data_filtrada['valor'] < 0) | (data_filtrada['valor'] > 300), 'valor'] = np.nan
# Mostrar cuántos valores fueron marcados como NaN
valores_fuera_de_rango = data_filtrada['valor'].isna().sum()
print(f"Se marcaron {valores_fuera_de_rango} valores como NaN por estar fuera del rango 0–300 mm.")
Se marcaron 865 valores como NaN por estar fuera del rango 0–300 mm.
Paso 6: Distribución gama por percentiles¶
In [11]:
import pandas as pd
import numpy as np
from scipy.stats import gamma
import matplotlib.pyplot as plt
# Asegurarse que 'data_filtrada' sea una copia independiente
data_filtrada = data_filtrada.copy()
# Convertir fecha si no está en datetime (esto previene errores)
data_filtrada['fecha_toma_dato'] = pd.to_datetime(data_filtrada['fecha_toma_dato'], errors='coerce')
# Extraer la hora de forma segura
data_filtrada.loc[:, 'hora'] = data_filtrada['fecha_toma_dato'].dt.hour
In [12]:
from scipy.stats import gamma
import numpy as np
def plot_gamma_por_hora(data, hora, ax):
# Filtrar solo valores válidos (no nulos y ≥ 0.1 mm)
data_hora = data[
(data['hora'] == hora) &
(data['valor'] >= 0.1) &
(~data['valor'].isna())
]
# Ajustar distribución Gamma
valores = data_hora['valor'].values
shape, loc, scale = gamma.fit(valores, floc=0) # loc = 0 es habitual en precipitación
# Percentiles de interés
percentiles = [0.95, 0.975, 0.99, 0.995]
percentiles_gamma = gamma.ppf(percentiles, shape, loc=loc, scale=scale)
# Histograma normalizado (barras azul #005ac8) + curva Gamma
x = np.linspace(0, valores.max(), 200)
ax.hist(
valores,
bins=30,
density=True,
alpha=0.4,
color='#005ac8', # tono azul solicitado
edgecolor='black'
)
ax.plot(
x,
gamma.pdf(x, shape, loc=loc, scale=scale),
color='orange',
label='Gamma ajustada'
)
# Líneas verticales para los percentiles
for p, val in zip(percentiles, percentiles_gamma):
ax.axvline(
val,
linestyle='--',
color='gray',
label=f'{int(p*100)}%'
)
# Estética
ax.set_title(f'Distribución Gamma - {hora}h')
ax.set_xlabel('Precipitación (mm)')
ax.set_ylabel('Densidad')
ax.legend()
In [13]:
def plot_gamma_por_hora(data, hora, ax):
# Filtrar solo valores válidos (no nulos y mayores o iguales a 0.1 mm)
data_hora = data[(data['hora'] == hora) & (data['valor'] >= 0.1) & (~data['valor'].isna())]
# Ajustar distribución gamma a los datos
valores = data_hora['valor'].values
shape, loc, scale = gamma.fit(valores, floc=0) # loc=0 es común en precipitaciones
# Calcular percentiles deseados
percentiles = [0.95, 0.975, 0.99, 0.995]
percentiles_gamma = gamma.ppf(percentiles, shape, loc=loc, scale=scale)
# Histograma normalizado + curva gamma
x = np.linspace(0, valores.max(), 200)
ax.hist(valores, bins=30, density=True, alpha=0.4, color='black', edgecolor='black')
ax.plot(x, gamma.pdf(x, shape, loc=loc, scale=scale), color='orange', label='Gamma ajustada')
# Líneas verticales para los percentiles
for p, val in zip(percentiles, percentiles_gamma):
ax.axvline(val, linestyle='--', color='gray', label=f'{int(p*100)}%')
# Estética
ax.set_title(f'Distribución Gamma - {hora}h')
ax.set_xlabel('Precipitación (mm)')
ax.set_ylabel('Densidad')
ax.legend()
In [14]:
fig, axes = plt.subplots(3, 1, figsize=(8, 14), sharey=True)
for i, hora in enumerate([0, 12, 18]):
plot_gamma_por_hora(data_filtrada, hora, axes[i])
axes[i].tick_params(axis='both', labelsize=16) # Tamaño de números en ejes
# Solo la última figura lleva etiqueta en el eje X
if hora == 18:
axes[i].set_xlabel("Precipitación (mm)", fontsize=16)
else:
axes[i].set_xlabel("")
axes[i].set_ylabel("Frecuencia", fontsize=16)
# Crear título multilínea centrado
titulo = (
f"{Nombre}\n"
f"Demarcación: {Demarcacion}, Provincia: {Provincia}, Cantón: {Canton}\n"
f"Latitud: {Latitud}°, Longitud: {Longitud}°, Altitud: {Altitud}"
)
plt.suptitle(titulo, x=0.52, y=0.96, fontsize=16)
# Ajustar espacio para suptitle
plt.tight_layout(rect=[0, 0, 1, 0.96])
# Guardar figura con nombre personalizado
plt.savefig(f"distribucion_gamma_por_hora_{Codigo}.png", dpi=200)
plt.show()
In [15]:
from scipy.stats import gamma
# Diccionario para guardar resultados
percentiles_resultado = {}
# Definir percentiles de interés
percentiles = [0.95, 0.975, 0.99, 0.995]
# Calcular y almacenar por hora
for hora in [0, 12, 18]:
data_hora = data_filtrada[(data_filtrada['hora'] == hora) &
(data_filtrada['valor'] >= 0.1) &
(~data_filtrada['valor'].isna())]
valores = data_hora['valor'].values
shape, loc, scale = gamma.fit(valores, floc=0)
valores_percentiles = gamma.ppf(percentiles, shape, loc=loc, scale=scale)
print(f"\nHora {hora}h:")
for p, val in zip(percentiles, valores_percentiles):
print(f" Percentil {int(p*100)}%: {val:.2f} mm")
Hora 0h: Percentil 95%: 4.72 mm Percentil 97%: 5.98 mm Percentil 99%: 7.66 mm Percentil 99%: 8.95 mm Hora 12h: Percentil 95%: 8.88 mm Percentil 97%: 11.24 mm Percentil 99%: 14.41 mm Percentil 99%: 16.83 mm Hora 18h: Percentil 95%: 4.15 mm Percentil 97%: 5.21 mm Percentil 99%: 6.63 mm Percentil 99%: 7.70 mm
In [16]:
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import gamma
# Calcular los máximos secos mensuales por hora
fig, axes = plt.subplots(1, 3, figsize=(18, 5), sharey=True)
for idx, hora in enumerate([0, 12, 18]):
data_hora = data_filtrada[(data_filtrada['hora'] == hora) &
(data_filtrada['valor'] >= 0.1) &
(~data_filtrada['valor'].isna())].copy()
# Agregar columna mes
data_hora['mes'] = pd.to_datetime(data_hora['fecha_toma_dato']).dt.month
# Ajustar gamma y calcular percentil 99.5% (0.995) por mes
percentil_995 = []
for mes in range(1, 13):
data_mes = data_hora[data_hora['mes'] == mes]['valor']
if len(data_mes) > 10:
shape, loc, scale = gamma.fit(data_mes, floc=0)
p_995 = gamma.ppf(0.995, shape, loc=loc, scale=scale)
else:
p_995 = None
percentil_995.append(p_995)
# Boxplot
sns.boxplot(ax=axes[idx], data=data_hora, x='mes', y='valor', color='lightyellow')
axes[idx].plot(range(12), percentil_995, color='gray', marker='o', linewidth=2)
axes[idx].set_title(f'Largo secuencia seca - {hora}h')
axes[idx].set_xlabel('Mes')
axes[idx].set_ylabel('Precipitación (mm)')
titulo = (
f"{Nombre}\n"
f"Demarcación: {Demarcacion}, Provincia: {Provincia}, Cantón: {Canton}\n"
f"Latitud: {Latitud}°, Longitud: {Longitud}°, Altitud: {Altitud}"
)
# Posición del título (modificable con x e y)
plt.suptitle(titulo, x=0.5, y=0.96, fontsize=12) # Centrado y un poco más arriba
plt.tight_layout()
plt.show()
In [17]:
# Crear boxplots por hora con línea del percentil 99.5% por hora
fig, ax = plt.subplots(figsize=(10, 6))
# Lista de horas y sus nombres
horas = [0, 12, 18]
labels = ['0h', '12h', '18h']
valores_plot = []
# Calcular valores y percentiles para cada hora
percentiles_995 = []
for hora in horas:
data_hora = data_filtrada[(data_filtrada['hora'] == hora) &
(data_filtrada['valor'] >= 0.1) &
(~data_filtrada['valor'].isna())]
valores = data_hora['valor'].values
valores_plot.append(valores)
# Ajustar distribución gamma y obtener percentil 99.5%
shape, loc, scale = gamma.fit(valores, floc=0)
p_995 = gamma.ppf(0.995, shape, loc=loc, scale=scale)
percentiles_995.append(p_995)
# Crear el boxplot
ax.boxplot(valores_plot, labels=labels, patch_artist=True, boxprops=dict(facecolor='lightyellow'))
# Agregar línea horizontal por cada percentil 99.5%
for i, p in enumerate(percentiles_995, start=1):
ax.hlines(p, i - 0.3, i + 0.3, colors='gray', linestyles='--', label=f'P99.5%' if i == 1 else None)
#ax.set_title('Distribución de precipitación por hora de observación')
ax.set_ylabel('Precipitación (mm)')
ax.set_xlabel('Hora')
ax.legend()
plt.tight_layout()
# Guardar la figura con 200 DPI
plt.savefig("Boxplot_preci", dpi=200)
titulo = (
f"{Nombre}\n"
f"Demarcación: {Demarcacion}, Provincia: {Provincia}, Cantón: {Canton}\n"
f"Latitud: {Latitud}°, Longitud: {Longitud}°, Altitud: {Altitud}"
)
# Posición del título (modificable con x e y)
plt.suptitle(titulo, x=0.54, y=0.96, fontsize=12) # Centrado y un poco más arriba
plt.tight_layout()
plt.show()
/tmp/ipykernel_175955/1599002823.py:27: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. ax.boxplot(valores_plot, labels=labels, patch_artist=True, boxprops=dict(facecolor='lightyellow'))
In [18]:
print(data_filtrada)
fecha_toma_dato valor hora 0 2008-04-30 12:00:00 0.3 12 1 1979-11-01 12:00:00 0.0 12 2 1979-11-02 12:00:00 0.0 12 3 1979-11-03 12:00:00 0.0 12 4 1979-11-04 12:00:00 0.0 12 ... ... ... ... 41914 2011-08-25 00:00:00 0.0 0 41915 2011-08-26 00:00:00 0.0 0 41916 2011-08-27 00:00:00 0.0 0 41917 2011-08-28 00:00:00 0.1 0 41918 2011-08-29 00:00:00 0.1 0 [41919 rows x 3 columns]
Paso 6: Z score modificado¶
In [19]:
import numpy as np
# Asegurarse de que 'fecha_toma_dato' sea datetime
data_filtrada['fecha_toma_dato'] = pd.to_datetime(data_filtrada['fecha_toma_dato'])
# Inicializar columnas
data_filtrada['z_score'] = np.nan
data_filtrada['sospechoso_z'] = False
# Evaluar por hora
for hora in [7, 13, 19]:
df_hora = data_filtrada[data_filtrada['hora'] == hora].copy()
df_hora = df_hora.sort_values('fecha_toma_dato').reset_index()
z_scores = []
for i, row in df_hora.iterrows():
fecha = row['fecha_toma_dato']
valor = row['valor']
ventana = df_hora[(df_hora['fecha_toma_dato'] >= fecha - pd.Timedelta(days=2)) &
(df_hora['fecha_toma_dato'] <= fecha + pd.Timedelta(days=2))]['valor'].dropna()
if len(ventana) < 5 or np.isnan(valor):
z_scores.append(np.nan)
continue
q1 = ventana.quantile(0.25)
q3 = ventana.quantile(0.75)
ssd = (q3 - q1) / 1.349
mediana = ventana.median()
z = (valor - mediana) / ssd if ssd > 0 else np.nan
z_scores.append(z)
# Guardar valores Z en la columna original
data_filtrada.loc[df_hora['index'], 'z_score'] = z_scores
# Calcular percentil 99.5% global de Z
umbral_995 = data_filtrada['z_score'].abs().quantile(0.995)
# Etiquetar como sospechoso si Z > P99.5
data_filtrada['sospechoso_z'] = data_filtrada['z_score'].abs() > umbral_995
print(f"Umbral Z (percentil 99.5%) calculado: {umbral_995:.2f}")
Umbral Z (percentil 99.5%) calculado: nan
In [20]:
import matplotlib.pyplot as plt
import numpy as np
fig, axes = plt.subplots(3, 1, figsize=(15, 10), sharex=True)
for idx, hora in enumerate([0, 12, 18]):
df_hora = data_filtrada[data_filtrada['hora'] == hora].copy()
df_hora = df_hora.sort_values('fecha_toma_dato')
# Cálculo rolling
rol = df_hora['valor'].rolling(window=5, center=True)
q1 = rol.quantile(0.25)
q3 = rol.quantile(0.75)
mediana = rol.median()
ssd = (q3 - q1) / 1.349
z = (df_hora['valor'] - mediana) / ssd
z[np.abs(z) > 10] = np.nan # Eliminar valores anómalos
# Umbral del percentil 99.5 %
z_percentil_995 = np.nanquantile(np.abs(z), 0.995)
# Serie de valores Z
axes[idx].plot(df_hora['fecha_toma_dato'], z,
label='Valor Z', color='blue')
# Líneas horizontales de referencia
for lim in [1.5, 2.0, 2.5]:
axes[idx].axhline(y=lim, color='gray', linestyle='--',
linewidth=1, label=f'lim:{lim}' if idx == 0 else "")
axes[idx].axhline(y=3.0, color='gray', linestyle='--',
linewidth=1.2, label='lim:3' if idx == 0 else "")
axes[idx].axhline(y=z_percentil_995, color='black', linestyle='-',
linewidth=2, label='lim:P99.5%' if idx == 0 else "")
axes[idx].set_ylim(0, 10)
axes[idx].set_title(f'Estadístico Z a lo largo del tiempo – {hora} h', fontsize=16)
axes[idx].set_ylabel('Z', fontsize=16)
axes[idx].tick_params(axis='both', labelsize=16)
axes[2].set_xlabel('Fecha', fontsize=16)
# Leyenda fuera del área de gráficos
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels,
loc='center left',
bbox_to_anchor=(0.85, 0.5),
frameon=False,
fontsize=16)
# Título principal
titulo = (
f"{Nombre}\n"
f"Demarcación: {Demarcacion}, Provincia: {Provincia}, Cantón: {Canton}\n"
f"Latitud: {Latitud}°, Longitud: {Longitud}°, Altitud: {Altitud}"
)
plt.suptitle(titulo, x=0.44, y=0.96, fontsize=16)
# Ajuste de layout para respetar el suptitle y la leyenda
fig.tight_layout(rect=[0, 0, 0.85, 0.96])
# Guardar con 200 DPI
plt.savefig(f"Z_score_preci_{Codigo}.png", dpi=200, bbox_inches='tight')
plt.show()
In [21]:
import matplotlib.pyplot as plt
import numpy as np
# Crear figura y ejes para 3 subgráficos (uno por hora)
fig, axes = plt.subplots(3, 1, figsize=(15, 10), sharex=True)
for idx, hora in enumerate([0, 12, 18]):
df_hora = data_filtrada[data_filtrada['hora'] == hora].copy()
df_hora = df_hora.sort_values('fecha_toma_dato')
# Calcular umbral del percentil 99.5% de precipitación para esta hora
umbral_p995 = df_hora['valor'].dropna().quantile(0.995)
# Graficar la serie de precipitación
axes[idx].plot(df_hora['fecha_toma_dato'], df_hora['valor'],
label='Precipitación (mm)', color='blue')
# Agregar línea horizontal del percentil 99.5%
axes[idx].axhline(y=umbral_p995, color='black', linestyle='-',
linewidth=2, label='Umbral P99.5%')
# Estética
axes[idx].set_ylabel('Precipitación (mm)', fontsize=16)
axes[idx].set_title(f'Precipitación a lo largo del tiempo – {hora} h', fontsize=16)
axes[idx].tick_params(axis='both', labelsize=16)
axes[idx].set_ylim(bottom=0)
# Solo la primera figura lleva leyenda
if idx == 0:
axes[idx].legend(fontsize=16)
# Etiqueta en eje X solo en el último panel
axes[2].set_xlabel('Fecha', fontsize=16)
# Título general con metadatos
titulo = (
f"{Nombre}\n"
f"Demarcación: {Demarcacion}, Provincia: {Provincia}, Cantón: {Canton}\n"
f"Latitud: {Latitud}°, Longitud: {Longitud}°, Altitud: {Altitud}"
)
plt.suptitle(titulo, x=0.54, y=0.96, fontsize=16)
# Ajustar el layout dejando espacio para el título
fig.tight_layout(rect=[0, 0, 1, 0.94])
# Guardar figura
plt.savefig(f"precipitacion_umbral_p995_{Codigo}.png", dpi=200, bbox_inches='tight')
# Mostrar figura
plt.show()
In [22]:
# Aplicar Isolation Forest con configuración personalizada y sin paralelización
from sklearn.ensemble import IsolationForest
data_filtrada['sospechoso_isolation'] = False
for hora in [0, 12, 18]:
df_hora = data_filtrada[(data_filtrada['hora'] == hora) & (~data_filtrada['valor'].isna())]
valores = df_hora[['valor']]
if len(valores) >= 20:
modelo = IsolationForest(
n_estimators=1000,
max_samples='auto',
contamination=0.005, # 99.5% normal
n_jobs=1,
random_state=123
)
modelo.fit(valores)
predicciones = modelo.predict(valores)
indices_sospechosos = df_hora.index[predicciones == -1]
data_filtrada.loc[indices_sospechosos, 'sospechoso_isolation'] = True
# Mostrar los primeros registros actualizados
print(data_filtrada.head())
fecha_toma_dato valor hora z_score sospechoso_z \ 0 2008-04-30 12:00:00 0.3 12 NaN False 1 1979-11-01 12:00:00 0.0 12 NaN False 2 1979-11-02 12:00:00 0.0 12 NaN False 3 1979-11-03 12:00:00 0.0 12 NaN False 4 1979-11-04 12:00:00 0.0 12 NaN False sospechoso_isolation 0 False 1 False 2 False 3 False 4 False
In [23]:
# Eliminar la columna z_score si existe
if 'z_score' in data_filtrada.columns:
data_filtrada = data_filtrada.drop(columns=['z_score'])
In [24]:
# Reaplicar Isolation Forest sin paralelización (n_jobs=1)
from sklearn.ensemble import IsolationForest
import matplotlib.pyplot as plt
data_filtrada['sospechoso_isolation'] = False
for hora in [0, 12, 18]:
df_hora = data_filtrada[(data_filtrada['hora'] == hora) & (~data_filtrada['valor'].isna())]
valores = df_hora[['valor']]
if len(valores) >= 20:
modelo = IsolationForest(
n_estimators=1000,
max_samples='auto',
contamination=0.005,
n_jobs=1,
random_state=123
)
modelo.fit(valores)
predicciones = modelo.predict(valores)
indices_sospechosos = df_hora.index[predicciones == -1]
data_filtrada.loc[indices_sospechosos, 'sospechoso_isolation'] = True
# Graficar puntos: normales en azul, sospechosos en rojo
fig, axes = plt.subplots(3, 1, figsize=(15, 10), sharex=True)
for idx, hora in enumerate([0, 12, 18]):
df_hora = data_filtrada[data_filtrada['hora'] == hora].copy()
df_hora = df_hora.sort_values('fecha_toma_dato')
normales = df_hora[df_hora['sospechoso_isolation'] == False]
sospechosos = df_hora[df_hora['sospechoso_isolation'] == True]
axes[idx].scatter(normales['fecha_toma_dato'], normales['valor'],
color='blue', s=10, label='Normal')
axes[idx].scatter(sospechosos['fecha_toma_dato'], sospechosos['valor'],
color='red', s=10, label='Sospechoso')
axes[idx].set_title(f'Isolation Forest - {hora}h', fontsize=16)
axes[idx].set_ylabel('Precipitación (mm)', fontsize=16)
axes[idx].tick_params(axis='both', labelsize=16)
# Mostrar leyenda solo en el primer subgráfico
if idx == 0:
axes[idx].legend(fontsize=16)
# Etiqueta solo en el último eje
axes[2].set_xlabel('Fecha', fontsize=16)
# Título general con metadatos
titulo = (
f"{Nombre}\n"
f"Demarcación: {Demarcacion}, Provincia: {Provincia}, Cantón: {Canton}\n"
f"Latitud: {Latitud}°, Longitud: {Longitud}°, Altitud: {Altitud}"
)
plt.suptitle(titulo, x=0.54, y=0.96, fontsize=16)
# Ajustar layout para no solapar título ni recortar
plt.tight_layout(rect=[0, 0, 1, 0.96])
# Guardar figura sin recortar
plt.savefig(f"Isolation_forest_preci_{Codigo}.png", dpi=200, bbox_inches='tight')
# Mostrar
plt.show()
In [25]:
# Paso 1: Detectar ≥ 3 días consecutivos con el mismo valor (≠ 0)
data_filtrada['repetidos_consecutivos'] = False
for hora in [0, 12, 18]:
df_hora = data_filtrada[data_filtrada['hora'] == hora].copy().sort_values('fecha_toma_dato')
valores = df_hora['valor'].values
indices = df_hora.index
for i in range(2, len(valores)):
if (
valores[i] == valores[i - 1] == valores[i - 2] and
valores[i] != 0 and
not np.isnan(valores[i])
):
data_filtrada.loc[indices[i], 'repetidos_consecutivos'] = True
# Paso 2: Detectar ≥ 3 días consecutivos con precipitación igual a 0
data_filtrada['rachas_cero'] = False
for hora in [0, 12, 18]:
df_hora = data_filtrada[data_filtrada['hora'] == hora].copy().sort_values('fecha_toma_dato')
valores = df_hora['valor'].values
indices = df_hora.index
for i in range(2, len(valores)):
if (
valores[i] == valores[i - 1] == valores[i - 2] == 0 and
not np.isnan(valores[i])
):
data_filtrada.loc[indices[i], 'rachas_cero'] = True
# Paso 3: Detectar gaps de más de 2 días sin datos
data_filtrada['lagunas_datos'] = False
for hora in [0, 12, 18]:
df_hora = data_filtrada[(data_filtrada['hora'] == hora)].copy().sort_values('fecha_toma_dato')
fechas = df_hora['fecha_toma_dato']
gaps = fechas.diff().dt.days.fillna(1)
indices_con_gap = df_hora.index[gaps > 2]
data_filtrada.loc[indices_con_gap, 'lagunas_datos'] = True
# Paso 4: Crear columna vacía para validación manual
data_filtrada['validacion_libretas'] = ''
print(data_filtrada)
fecha_toma_dato valor hora sospechoso_z sospechoso_isolation \
0 2008-04-30 12:00:00 0.3 12 False False
1 1979-11-01 12:00:00 0.0 12 False False
2 1979-11-02 12:00:00 0.0 12 False False
3 1979-11-03 12:00:00 0.0 12 False False
4 1979-11-04 12:00:00 0.0 12 False False
... ... ... ... ... ...
41914 2011-08-25 00:00:00 0.0 0 False False
41915 2011-08-26 00:00:00 0.0 0 False False
41916 2011-08-27 00:00:00 0.0 0 False False
41917 2011-08-28 00:00:00 0.1 0 False False
41918 2011-08-29 00:00:00 0.1 0 False False
repetidos_consecutivos rachas_cero lagunas_datos validacion_libretas
0 False False False
1 False False False
2 False False False
3 False True False
4 False True False
... ... ... ... ...
41914 False True False
41915 False True False
41916 False True False
41917 False False False
41918 False False False
[41919 rows x 9 columns]
In [26]:
# Generar nombre de archivo con base en el código de estación
nombre_archivo = f"control_precipitacion_observaciones_2025_{Codigo}.csv"
ruta_salida = f"{nombre_archivo}"
# Guardar el DataFrame con el nombre personalizado
data_filtrada.to_csv(ruta_salida, index=False)
ruta_salida
Out[26]:
'control_precipitacion_observaciones_2025_M0258.csv'