Análisis de Sismicidad de Colombia

import sys

# ── Detección automática del entorno de ejecución ────────────────────────────
_ES_JUPYTERLITE = 'pyodide' in sys.modules
_ES_COLAB       = 'google.colab' in sys.modules

if _ES_JUPYTERLITE:
    import micropip
    await micropip.install(['pyodide-http', 'folium', 'scipy', 'seaborn'])
    import pyodide_http
    pyodide_http.patch_all()
    print('✅ Entorno: JupyterLite (navegador)')

else:
    import subprocess
    paquetes = [
        'pandas<3.0.0',   # <3.0.0 por compatibilidad con bqplot y otros paquetes del entorno
        'numpy',
        'matplotlib',
        'seaborn',
        'folium',
        'requests',
        'scipy',
    ]
    print('📦 Verificando e instalando dependencias...')
    subprocess.run(
        [sys.executable, '-m', 'pip', 'install', '-q'] + paquetes,
        check=True
    )
    entorno = 'Google Colab' if _ES_COLAB else 'local / Binder'
    print(f'✅ Entorno: {entorno}')

import requests
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
import folium
from folium.plugins import MarkerCluster, HeatMap
from scipy import stats
from datetime import datetime
import warnings
warnings.filterwarnings('ignore')

# Configuración de estilos
plt.rcParams.update({
    'figure.dpi': 120,
    'font.family': 'sans-serif',
    'axes.spines.top': False,
    'axes.spines.right': False,
})
UNGRD_BLUE = '#0154a5'
UNGRD_ACCENT = '#2563eb'

print('✅ Dependencias cargadas correctamente')

# Parámetros de consulta — puedes modificar estos valores
PARAMS = {
    'format':     'csv',
    'starttime':  '2000-01-01',
    'endtime':    '2024-12-31',
    'minlatitude':  -5.0,   # sur de Colombia
    'maxlatitude':  13.0,   # norte (incluye Caribe)
    'minlongitude': -82.0,  # oeste (Pacífico)
    'maxlongitude': -65.0,  # este (Orinoquía)
    'minmagnitude': 3.0,
    'orderby': 'time',
}

API_URL = 'https://earthquake.usgs.gov/fdsnws/event/1/query'

print('🌐 Descargando catálogo sísmico del USGS...')
print(f"   Período: {PARAMS['starttime']} → {PARAMS['endtime']}")
print(f"   Área: Colombia y área de influencia")
print(f"   Magnitud mínima: M ≥ {PARAMS['minmagnitude']}")

response = requests.get(API_URL, params=PARAMS, timeout=120)
response.raise_for_status()

from io import StringIO
df_raw = pd.read_csv(StringIO(response.text))
print(f"\n✅ Descarga completa: {len(df_raw):,} eventos sísmicos")

# Seleccionar y renombrar columnas relevantes
df = df_raw[['time', 'latitude', 'longitude', 'depth', 'mag', 'magType', 'place']].copy()
df.columns = ['tiempo', 'latitud', 'longitud', 'profundidad_km', 'magnitud', 'tipo_mag', 'lugar']

# Convertir tiempo a datetime
df['tiempo'] = pd.to_datetime(df['tiempo'], utc=True)
df['año']    = df['tiempo'].dt.year
df['mes']    = df['tiempo'].dt.month

# Eliminar registros sin magnitud o profundidad
df = df.dropna(subset=['magnitud', 'profundidad_km', 'latitud', 'longitud'])

# Clasificación de profundidad focal
def clasificar_profundidad(d):
    if d <= 70:   return 'Superficial (≤ 70 km)'
    elif d <= 300: return 'Intermedia (70–300 km)'
    else:          return 'Profunda (> 300 km)'

df['clase_profundidad'] = df['profundidad_km'].apply(clasificar_profundidad)

print(f"Registros válidos: {len(df):,}")
print(f"\nEstadísticas descriptivas:")
df[['magnitud', 'profundidad_km']].describe().round(2)

# Distribución por clase de profundidad
conteo_prof = df['clase_profundidad'].value_counts()
print('Distribución por tipo de sismo:')
for clase, n in conteo_prof.items():
    print(f'  {clase}: {n:,} ({n/len(df)*100:.1f}%)')

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# --- Panel superior: sismos por año ---
sismos_anuales = df.groupby('año').size().reset_index(name='conteo')
axes[0].bar(sismos_anuales['año'], sismos_anuales['conteo'],
            color=UNGRD_BLUE, alpha=0.85, edgecolor='white', linewidth=0.5)
axes[0].axhline(sismos_anuales['conteo'].mean(), color='tomato',
                linestyle='--', linewidth=1.5, label=f"Media anual: {sismos_anuales['conteo'].mean():.0f}")
axes[0].set_title('Sismicidad anual en Colombia (M ≥ 3.0)', fontsize=13, fontweight='bold', pad=10)
axes[0].set_ylabel('Número de eventos', fontsize=11)
axes[0].legend(fontsize=10)
axes[0].set_xlim(1999.5, 2024.5)

# --- Panel inferior: magnitud máxima por año ---
mag_anual = df.groupby('año')['magnitud'].max().reset_index()
sc = axes[1].scatter(mag_anual['año'], mag_anual['magnitud'],
                     s=mag_anual['magnitud']**2.5 * 3, c=mag_anual['magnitud'],
                     cmap='YlOrRd', alpha=0.85, edgecolors=UNGRD_BLUE, linewidths=0.5)
plt.colorbar(sc, ax=axes[1], label='Magnitud máxima')
axes[1].set_title('Magnitud máxima registrada por año', fontsize=13, fontweight='bold', pad=10)
axes[1].set_ylabel('Magnitud (Mw)', fontsize=11)
axes[1].set_xlabel('Año', fontsize=11)
axes[1].set_xlim(1999.5, 2024.5)

plt.tight_layout()
plt.savefig('sismicidad_temporal.png', bbox_inches='tight', dpi=150)
plt.show()
print('Figura guardada: sismicidad_temporal.png')

def calcular_gutenberg_richter(magnitudes, m_min=3.0, bin_size=0.1):
    """Calcula la curva Gutenberg-Richter y estima la b-value por mínimos cuadrados."""
    bins = np.arange(m_min, magnitudes.max() + bin_size, bin_size)
    counts, edges = np.histogram(magnitudes, bins=bins)
    cum_counts = np.cumsum(counts[::-1])[::-1]  # frecuencia acumulada
    mag_centers = (edges[:-1] + edges[1:]) / 2

    # Solo usar bins con al menos 1 evento
    mask = cum_counts > 0
    log_n = np.log10(cum_counts[mask])
    mags  = mag_centers[mask]

    # Regresión lineal para estimar a y b
    slope, intercept, r, p, se = stats.linregress(mags, log_n)
    b_value = -slope
    a_value = intercept

    return mags, log_n, a_value, b_value, r**2


mags_gr, log_n_gr, a, b, r2 = calcular_gutenberg_richter(df['magnitud'])

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# --- Gráfico G-R ---
axes[0].scatter(mags_gr, log_n_gr, color=UNGRD_BLUE, s=25, alpha=0.8,
                label='Datos observados', zorder=3)
m_fit = np.linspace(mags_gr.min(), mags_gr.max(), 100)
axes[0].plot(m_fit, a - b * m_fit, color='tomato', linewidth=2,
             label=f'Ajuste G-R: b = {b:.2f}  (R² = {r2:.3f})')
axes[0].set_xlabel('Magnitud (M)', fontsize=11)
axes[0].set_ylabel('log₁₀ N (frecuencia acumulada)', fontsize=11)
axes[0].set_title('Relación de Gutenberg-Richter\nColombia 2000–2024', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# --- Histograma de magnitudes ---
axes[1].hist(df['magnitud'], bins=np.arange(3.0, df['magnitud'].max() + 0.2, 0.2),
             color=UNGRD_BLUE, alpha=0.85, edgecolor='white', linewidth=0.4)
axes[1].axvline(df['magnitud'].median(), color='tomato', linestyle='--',
                linewidth=1.5, label=f"Mediana: M {df['magnitud'].median():.1f}")
axes[1].set_xlabel('Magnitud (M)', fontsize=11)
axes[1].set_ylabel('Frecuencia', fontsize=11)
axes[1].set_title('Distribución de Magnitudes', fontsize=12, fontweight='bold')
axes[1].legend(fontsize=10)

plt.tight_layout()
plt.savefig('gutenberg_richter.png', bbox_inches='tight', dpi=150)
plt.show()

print(f'📊 Parámetros Gutenberg-Richter para Colombia:')
print(f'   b-value: {b:.3f}  (valor global típico ≈ 1.0)')
print(f'   a-value: {a:.3f}')
print(f'   R²:      {r2:.3f}')

COLORES_PROF = {
    'Superficial (≤ 70 km)':    '#e74c3c',
    'Intermedia (70–300 km)':   '#f39c12',
    'Profunda (> 300 km)':      '#2980b9',
}

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# --- Histograma de profundidades ---
for clase, color in COLORES_PROF.items():
    subset = df[df['clase_profundidad'] == clase]
    axes[0].hist(subset['profundidad_km'], bins=50, alpha=0.7, color=color,
                 label=f'{clase} (n={len(subset):,})', edgecolor='none')
axes[0].set_xlabel('Profundidad focal (km)', fontsize=11)
axes[0].set_ylabel('Frecuencia', fontsize=11)
axes[0].set_title('Distribución de Profundidades Focales', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=9)

# --- Magnitud vs. Profundidad ---
for clase, color in COLORES_PROF.items():
    subset = df[df['clase_profundidad'] == clase]
    axes[1].scatter(subset['profundidad_km'], subset['magnitud'],
                    alpha=0.15, s=8, color=color, label=clase)
axes[1].set_xlabel('Profundidad focal (km)', fontsize=11)
axes[1].set_ylabel('Magnitud (M)', fontsize=11)
axes[1].set_title('Relación Profundidad–Magnitud', fontsize=12, fontweight='bold')
axes[1].legend(fontsize=9, markerscale=3)

plt.tight_layout()
plt.savefig('profundidad_focal.png', bbox_inches='tight', dpi=150)
plt.show()

# Top 15 eventos por magnitud
top_sismos = df.nlargest(15, 'magnitud')[['tiempo', 'magnitud', 'profundidad_km',
                                           'clase_profundidad', 'latitud', 'longitud', 'lugar']]
top_sismos['tiempo'] = top_sismos['tiempo'].dt.strftime('%Y-%m-%d %H:%M UTC')
top_sismos = top_sismos.reset_index(drop=True)
top_sismos.index += 1
top_sismos.columns = ['Fecha/Hora', 'Magnitud', 'Prof. (km)', 'Clase', 'Lat', 'Lon', 'Lugar']
top_sismos

def color_por_profundidad(profundidad):
    if profundidad <= 70:    return '#e74c3c'   # rojo — superficial
    elif profundidad <= 300: return '#f39c12'   # naranja — intermedia
    else:                    return '#2980b9'   # azul — profunda

# Centrar mapa en Colombia
mapa = folium.Map(
    location=[4.5, -74.0],
    zoom_start=5,
    tiles='CartoDB positron'
)

# Usar solo M >= 4.0 para no sobrecargar el mapa
df_mapa = df[df['magnitud'] >= 4.0].copy()

# Capa de calor (heatmap) de todos los eventos M >= 3.0
heat_data = df[['latitud', 'longitud', 'magnitud']].values.tolist()
HeatMap(heat_data, name='Densidad sísmica (M ≥ 3.0)',
        radius=6, blur=8, max_zoom=8).add_to(mapa)

# Capa de marcadores para M >= 4.0
cluster = MarkerCluster(name='Epicentros M ≥ 4.0')
for _, row in df_mapa.iterrows():
    radio = max(3, row['magnitud'] ** 2.2 * 0.5)
    popup_txt = (
        f"<b>Magnitud:</b> M {row['magnitud']:.1f}<br>"
        f"<b>Profundidad:</b> {row['profundidad_km']:.0f} km ({row['clase_profundidad']})<br>"
        f"<b>Fecha:</b> {str(row['tiempo'])[:10]}<br>"
        f"<b>Lugar:</b> {row['lugar']}"
    )
    folium.CircleMarker(
        location=[row['latitud'], row['longitud']],
        radius=radio,
        color=color_por_profundidad(row['profundidad_km']),
        fill=True, fill_opacity=0.65, weight=0.5,
        popup=folium.Popup(popup_txt, max_width=280),
        tooltip=f"M {row['magnitud']:.1f}"
    ).add_to(cluster)

cluster.add_to(mapa)

# Leyenda
leyenda_html = """
<div style="position:fixed; bottom:30px; left:30px; z-index:1000;
     background:white; padding:12px 16px; border-radius:8px;
     box-shadow:0 2px 8px rgba(0,0,0,0.2); font-size:12px; font-family:sans-serif;">
  <b>Profundidad focal</b><br>
  <span style='color:#e74c3c;'>●</span> Superficial ≤ 70 km<br>
  <span style='color:#f39c12;'>●</span> Intermedia 70–300 km<br>
  <span style='color:#2980b9;'>●</span> Profunda > 300 km<br>
  <hr style='margin:6px 0; border-color:#eee;'>
  <small>Tamaño ∝ Magnitud · Catálogo USGS 2000–2024</small>
</div>
"""
mapa.get_root().html.add_child(folium.Element(leyenda_html))
folium.LayerControl().add_to(mapa)

# Guardar y mostrar
mapa.save('mapa_sismicidad_colombia.html')
print(f'✅ Mapa generado con {len(df_mapa):,} epicentros (M ≥ 4.0)')
print('   Abre mapa_sismicidad_colombia.html para visualización interactiva')
mapa

from folium.plugins import TimestampedGeoJson

# Filtrar eventos M >= 4.0 para reducir el volumen de la animación
df_anim = df[df['magnitud'] >= 4.0].copy()
df_anim = df_anim.sort_values('tiempo').reset_index(drop=True)

# Construir features GeoJSON con timestamps individuales
features = []
for _, row in df_anim.iterrows():
    mag  = row['magnitud']
    prof = row['profundidad_km']
    color = color_por_profundidad(prof)
    feature = {
        'type': 'Feature',
        'geometry': {
            'type': 'Point',
            'coordinates': [row['longitud'], row['latitud']]
        },
        'properties': {
            'time': row['tiempo'].strftime('%Y-%m-%dT%H:%M:%S'),
            'popup': (
                f"<b>M {mag:.1f}</b><br>"
                f"Profundidad: {prof:.0f} km ({row['clase_profundidad']})<br>"
                f"Fecha: {str(row['tiempo'])[:10]}<br>"
                f"{row['lugar']}"
            ),
            'style': {'weight': 0},
            'icon': 'circle',
            'iconstyle': {
                'fillColor': color,
                'fillOpacity': 0.75,
                'stroke': True,
                'color': '#222',
                'weight': 0.4,
                'radius': max(4, mag ** 2.2 * 0.55)
            }
        }
    }
    features.append(feature)

# Mapa base con fondo oscuro para mejor contraste
mapa_anim = folium.Map(
    location=[4.5, -74.0],
    zoom_start=5,
    tiles='CartoDB dark_matter'
)

# Capa animada temporal — período mensual, ventana de 3 meses visible
TimestampedGeoJson(
    {'type': 'FeatureCollection', 'features': features},
    period='P1M',
    duration='P3M',
    transition_time=400,
    loop=False,
    auto_play=False,
    add_last_point=True,
    date_options='YYYY-MM-DD'
).add_to(mapa_anim)

# Leyenda
leyenda_anim = """
<div style='position:fixed; bottom:30px; left:30px; z-index:1000;
     background:rgba(20,20,20,0.85); color:#eee; padding:12px 16px;
     border-radius:8px; box-shadow:0 2px 8px rgba(0,0,0,0.5);
     font-size:12px; font-family:sans-serif;'>
  <b>Profundidad focal</b><br>
  <span style='color:#e74c3c;'>●</span> Superficial ≤ 70 km<br>
  <span style='color:#f39c12;'>●</span> Intermedia 70–300 km<br>
  <span style='color:#2980b9;'>●</span> Profunda > 300 km<br>
  <hr style='margin:6px 0; border-color:#555;'>
  <small>Tamaño ∝ Magnitud · M ≥ 4.0 · 2000–2024</small>
</div>
"""
mapa_anim.get_root().html.add_child(folium.Element(leyenda_anim))

mapa_anim.save('mapa_animado_sismicidad.html')
print(f'✅ Animación generada: {len(df_anim):,} eventos M ≥ 4.0 · {len(features)} frames mensuales')
mapa_anim

meses_nombres = ['Ene', 'Feb', 'Mar', 'Abr', 'May', 'Jun',
                 'Jul', 'Ago', 'Sep', 'Oct', 'Nov', 'Dic']

fig, axes = plt.subplots(1, 2, figsize=(13, 4))

# --- Sismos por mes (todos los años) ---
por_mes = df.groupby('mes').size().reset_index(name='conteo')
axes[0].bar(por_mes['mes'], por_mes['conteo'],
            color=[plt.cm.Blues(0.4 + 0.5 * i / 12) for i in range(12)],
            edgecolor='white')
axes[0].axhline(por_mes['conteo'].mean(), color='tomato',
                linestyle='--', linewidth=1.5, label='Media mensual')
axes[0].set_xticks(range(1, 13))
axes[0].set_xticklabels(meses_nombres)
axes[0].set_title('Distribución mensual de sismos\n(suma 2000–2024)', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Número de eventos')
axes[0].legend()

# --- Magnitud media mensual ---
mag_mes = df.groupby('mes')['magnitud'].mean()
axes[1].plot(mag_mes.index, mag_mes.values, 'o-',
             color=UNGRD_BLUE, linewidth=2, markersize=7)
axes[1].fill_between(mag_mes.index, mag_mes.values,
                     mag_mes.values.min(), alpha=0.1, color=UNGRD_BLUE)
axes[1].set_xticks(range(1, 13))
axes[1].set_xticklabels(meses_nombres)
axes[1].set_title('Magnitud promedio mensual', fontsize=12, fontweight='bold')
axes[1].set_ylabel('Magnitud media (M)')

plt.tight_layout()
plt.savefig('variacion_estacional.png', bbox_inches='tight', dpi=150)
plt.show()

print('=' * 60)
print('  RESUMEN — SISMICIDAD DE COLOMBIA 2000–2024')
print('=' * 60)
print(f'  Total de eventos M ≥ 3.0:     {len(df):,}')
print(f'  Magnitud máxima registrada:   M {df["magnitud"].max():.1f}')
print(f'  Magnitud media:               M {df["magnitud"].mean():.2f}')
print(f'  Profundidad media:            {df["profundidad_km"].mean():.1f} km')
print(f'  Eventos superficiales (≤70):  {(df["profundidad_km"]<=70).sum():,} ({(df["profundidad_km"]<=70).mean()*100:.1f}%)')
print(f'  Período de análisis:          2000–2024 (25 años)')
print(f'  b-value Gutenberg-Richter:    {b:.3f}')
print('=' * 60)
print()
print('Archivos generados:')
print('  📊 sismicidad_temporal.png')
print('  📊 gutenberg_richter.png')
print('  📊 profundidad_focal.png')
print('  📊 variacion_estacional.png')
print('  🗺️  mapa_sismicidad_colombia.html')