Cuaderno 3: Visualización y comparación con 1979

Cuaderno 3: Visualización y comparación con 1979¶

Mapas, series de tiempo, control de calidad de las observaciones, comparación con 1979 y animación. Requiere los resultados del cuaderno 2 (output/resultados.npz). Runtime: CPU, 1–2 minutos.

Qué es y qué no es la comparación. Las ocho observaciones de la tabla data/observaciones_1979.csv están sin verificar contra las fuentes primarias, y varias coordenadas son aproximadas (dos quedan tierra adentro). El modelo es lineal y sin inundación, de modo que su η en el mar no es un run-up. Se muestran las diferencias para ubicar el orden de magnitud y los límites del modelo, no para validarlo.

import sys, subprocess, warnings
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

IN_COLAB = 'google.colab' in sys.modules
if IN_COLAB:
    subprocess.run([sys.executable, '-m', 'pip', 'install', '-q', 'rasterio', 'pyproj'], check=True)
import rasterio
from rasterio.crs import CRS
from pyproj import Transformer

WORK = Path('/content') if IN_COLAB else Path('.')
OUT = WORK / 'output'
OUT.mkdir(exist_ok=True)
RAW = 'https://raw.githubusercontent.com/scr-ungrd/cuadernos-reproducibles-crc003/main/data/'

def datos(nombre):
    """Ruta de un archivo de data/: local (../data) o descargado del repositorio (Colab)."""
    for base in (Path('../data'), Path('data'), WORK / 'data'):
        if (base / nombre).exists():
            return base / nombre
    import requests
    destino = WORK / 'data'
    destino.mkdir(exist_ok=True)
    r = requests.get(RAW + nombre, timeout=120)
    r.raise_for_status()
    (destino / nombre).write_bytes(r.content)
    return destino / nombre

# Epicentro del 12-dic-1979 (1,598° N, 79,359° O, según el artículo; por verificar con la fuente primaria)
EPI = (-79.359, 1.598)                       # (lon, lat)
CRS_UTM = CRS.from_epsg(32617)               # UTM 17N, meridiano central 81° O (error de escala < 0,2 % en el dominio)
TR = Transformer.from_crs(4326, 32617, always_xy=True)
G = 9.81

def ejes_lonlat(ax, paso=1):
    """Marcas de longitud y latitud sobre un mapa en coordenadas UTM (km); oeste a la izquierda, norte arriba."""
    lons = np.arange(-81, -76.5, paso); lats = np.arange(0, 4.5, paso)
    xl, yl = ax.get_xlim(), ax.get_ylim()
    ax.set_xticks([TR.transform(l, 2.0)[0] / 1e3 for l in lons]); ax.set_xticklabels([f'{abs(l):g}°O' for l in lons])
    ax.set_yticks([TR.transform(-79.0, l)[1] / 1e3 for l in lats]); ax.set_yticklabels([f'{l:g}°N' for l in lats])
    ax.set_xlim(xl); ax.set_ylim(yl)
    ax.set_xlabel('Longitud'); ax.set_ylabel('Latitud')

from matplotlib.animation import FuncAnimation, PillowWriter
if not (OUT / 'resultados.npz').exists():
    raise FileNotFoundError('Falta output/resultados.npz: ejecute primero el cuaderno 2.')
R = np.load(OUT / 'resultados.npz')
H = R['H'].astype('float64'); maxeta, tarr = R['maxeta'], R['tarr']
snap_t, snaps, rec, dt = R['snap_t'], R['snaps'], R['rec'], float(R['dt'])
x0, y0, dx = float(R['x0']), float(R['y0']), float(R['dx'])
ny, nx = H.shape
ext = [x0 / 1e3, (x0 + dx * nx) / 1e3, (y0 - dx * ny) / 1e3, y0 / 1e3]   # oeste a la izquierda, norte arriba
est = pd.read_csv(OUT / 'estaciones.csv'); obs = pd.read_csv(datos('observaciones_1979.csv'))
ex, ey = TR.transform(*EPI)
tt = np.arange(rec.shape[0]) * dt
print(f'Resultados: malla {nx}×{ny}, {len(snap_t)} instantáneas (0–{snap_t[-1]/60:.0f} min), {rec.shape[1]} estaciones en el mar')

# Prueba de orientación y de fuente: el máximo de η inicial debe estar cerca del epicentro (decenas de km), no a cientos
eta_ini = snaps[0].astype('float64')
i, j = np.unravel_index(np.argmax(eta_ini), eta_ini.shape)
d_max = np.hypot(x0 + dx * (j + .5) - ex, y0 - dx * (i + .5) - ey) / 1e3
assert d_max < 150, f'la fuente está a {d_max:.0f} km del epicentro'
print(f'✓ El máximo de η inicial ({eta_ini.max():.1f} m) está a {d_max:.0f} km del epicentro')

Resultados: malla 1065×1084, 31 instantáneas (0–60 min), 4 estaciones en el mar
✓ El máximo de η inicial (3.1 m) está a 75 km del epicentro

1. Control de calidad de las observaciones¶

Antes de comparar, se revisa que las observaciones sean físicamente coherentes con una onda larga $c=\sqrt{gH}$ : con la distancia de círculo máximo $d$ al epicentro y el arribo $t$ , la profundidad equivalente es $H_{eq}=(d/t)^2/g$ . En el océano frente a Colombia las profundidades son de 3 000–4 000 m (hasta ~5 000 m en el dominio); una profundidad equivalente mayor sobre una ruta que cruza plataforma y talud indica un dato sospechoso. Es el mismo control del cuaderno 318 del libro Teoría y métodos de análisis de riesgos.

def dist_km(lat1, lon1, lat2, lon2):
    p1, p2 = np.radians(lat1), np.radians(lat2)
    a = np.sin((p2 - p1) / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(np.radians(lon2 - lon1) / 2) ** 2
    return 2 * 6371.0 * np.arcsin(np.sqrt(a))

qc = obs.copy()
qc['dist_km'] = dist_km(EPI[1], EPI[0], qc.lat, qc.lon).round(0)
qc['v_media_ms'] = (qc.dist_km * 1e3 / (qc.tiempo_arr_min * 60)).round(0)
qc['H_equiv_m'] = (qc.v_media_ms ** 2 / G).round(-1)
qc = qc.merge(est[['estacion', 'estado']], on='estacion')
dominio = (qc.lon.between(-81.2, -76.9)) & (qc.lat.between(-0.2, 4.2))
qc.loc[~dominio, 'estado'] = 'fuera del dominio'
qc['sospechoso'] = np.where(qc.H_equiv_m > 4000, 'H equiv. > 4 000 m', '')
qc[['estacion', 'dist_km', 'runup_m', 'tiempo_arr_min', 'v_media_ms', 'H_equiv_m', 'sospechoso', 'estado', 'fuente']]

2. Elevación máxima de la ola¶

def fondo(ax):
    ax.imshow(np.where(H > 0, 0.97, 0.86), cmap='gray', vmin=0, vmax=1, extent=ext)
def costa(ax):
    ax.contour(np.linspace(ext[0], ext[1], nx), np.linspace(ext[3], ext[2], ny), H, levels=[0.5], colors='k', linewidths=.5)

fig, axs = plt.subplots(1, 2, figsize=(15, 7.5))
ax = axs[0]; fondo(ax)
ax.contour(np.linspace(ext[0], ext[1], nx), np.linspace(ext[3], ext[2], ny), H, levels=[200, 1000, 3000], colors='#8c99a6', linewidths=.4)
im = ax.imshow(np.where(H >= 10, np.maximum(maxeta, 0), np.nan), cmap='YlOrRd', vmin=0, vmax=5, extent=ext, alpha=.9)
costa(ax); plt.colorbar(im, ax=ax, shrink=.75, label='η máxima en 60 min, profundidad ≥ 10 m (m)')
ax.plot(ex / 1e3, ey / 1e3, '*', c='red', mec='k', ms=18, label='Epicentro 1979'); ax.legend(loc='upper left'); ejes_lonlat(ax)
ax.set_title('Elevación máxima de la superficie del mar (LSWE, 0–60 min)')

ax = axs[1]; fondo(ax)
zona = (H >= 10) & (H < 200)
im2 = ax.imshow(np.where(zona, np.maximum(maxeta, 0), np.nan), cmap='YlOrRd', vmin=0, vmax=8, extent=ext)
costa(ax); plt.colorbar(im2, ax=ax, shrink=.75, label='η máxima, profundidad 10–200 m (m)')
tx, ty = TR.transform(-78.80, 1.85)
ax.set_xlim(tx / 1e3 - 55, tx / 1e3 + 55); ax.set_ylim(ty / 1e3 - 50, ty / 1e3 + 50)
for _, r in qc.iterrows():
    x, y = TR.transform(r.lon, r.lat)
    if abs(x - tx) < 55e3 and abs(y - ty) < 50e3:
        ax.plot(x / 1e3, y / 1e3, 'ko', ms=4); ax.annotate(r.estacion, (x / 1e3, y / 1e3), textcoords='offset points', xytext=(4, 4), fontsize=8)
ejes_lonlat(ax, paso=.25); ax.set_title('Zona costera de Tumaco (η máx. en plataforma)')
plt.tight_layout(); plt.savefig(WORK / 'max_eta_map.png', dpi=130, bbox_inches='tight'); plt.show()
print(f'η máxima (profundidad ≥ 10 m): {maxeta[H>=10].max():.2f} m | plataforma 10–200 m: {maxeta[zona].max():.2f} m')
print(f'(en celdas de < 10 m el modelo lineal llega a {maxeta[H>0].max():.0f} m: amplificación sin fricción ni rotura, no se interpreta)')

η máxima (profundidad ≥ 10 m): 8.93 m | plataforma 10–200 m: 8.93 m
(en celdas de < 10 m el modelo lineal llega a 25 m: amplificación sin fricción ni rotura, no se interpreta)

3. Series de tiempo y arribos¶

El arribo simulado es el primer instante con η > 0,1 m (primer ascenso del nivel; en Tumaco la deformación inicial es un descenso y el primer ascenso es la llegada de la ola). La lectura se toma en la celda de mar descrita en el cuaderno 2.

ok = est[est.estado == 'ok'].reset_index(drop=True)
sim = []
for k, r in ok.iterrows():
    s = rec[:, k]; ta = tt[np.argmax(s > 0.1)] / 60 if (s > 0.1).any() else np.nan
    sim.append(dict(estacion=r.estacion, H_celda_m=r.H_celda_m, eta_max_m=s.max(), eta_min_m=s.min(),
                    t_arribo_sim_min=ta, t_cresta_min=tt[np.argmax(s)] / 60))
sim = pd.DataFrame(sim).merge(obs[['estacion', 'runup_m', 'tiempo_arr_min']], on='estacion')
sim['razon_amplitud'] = (sim.eta_max_m / sim.runup_m).round(2)
sim['dif_arribo_min'] = (sim.t_arribo_sim_min - sim.tiempo_arr_min).round(1)

fig, axs = plt.subplots(2, 2, figsize=(14, 7.5), sharex=True)
for ax, (_, r), k in zip(axs.flat, sim.iterrows(), range(len(sim))):
    s = rec[:, k]
    ax.plot(tt / 60, s, c='#0154a5', lw=1.6)
    ax.axhline(0, c='k', lw=.4)
    ax.axvline(r.tiempo_arr_min, c='#c0392b', ls='--', lw=1.2, label=f'Arribo observado ({r.tiempo_arr_min:.0f} min)*')
    ax.axvline(r.t_arribo_sim_min, c='#0154a5', ls=':', lw=1.2, label=f'Arribo simulado ({r.t_arribo_sim_min:.0f} min)')
    ax.axhline(r.runup_m, c='#c0392b', ls='-.', lw=.8, alpha=.6, label=f'Run-up observado ({r.runup_m} m)*')
    ax.set_title(f'{r.estacion} (celda de {r.H_celda_m:.0f} m de profundidad)', fontsize=10); ax.set_ylabel('η (m)'); ax.legend(fontsize=7, loc='upper left'); ax.grid(alpha=.3)
for ax in axs[1]: ax.set_xlabel('Tiempo desde el sismo (min)')
fig.text(.01, .005, '* Observaciones de data/observaciones_1979.csv, sin verificar contra la fuente primaria.', fontsize=8)
plt.tight_layout(); plt.savefig(WORK / 'time_series.png', dpi=130, bbox_inches='tight'); plt.show()
sim.round(2)

4. Comparación indicativa con las observaciones¶

fig, axs = plt.subplots(1, 2, figsize=(13, 5))
ax = axs[0]
ax.plot([0, 40], [0, 40], 'k--', lw=1, label='1:1')
ax.scatter(sim.tiempo_arr_min, sim.t_arribo_sim_min, s=70, c='#0154a5', zorder=3)
for _, r in sim.iterrows(): ax.annotate(r.estacion, (r.tiempo_arr_min, r.t_arribo_sim_min), textcoords='offset points', xytext=(5, 4), fontsize=8)
ax.set(xlim=(5, 30), ylim=(5, 30), xlabel='Arribo observado* (min)', ylabel='Arribo simulado (min)', title='Tiempos de arribo'); ax.legend(); ax.grid(alpha=.3)
ax = axs[1]; xp = np.arange(len(sim)); w = .38
ax.bar(xp - w / 2, sim.runup_m, w, label='Run-up observado*', color='#c0392b', alpha=.85)
ax.bar(xp + w / 2, sim.eta_max_m, w, label='η máxima simulada (mar, sin inundación)', color='#0154a5', alpha=.85)
ax.set_xticks(xp); ax.set_xticklabels(sim.estacion, rotation=20, ha='right', fontsize=8); ax.set_ylabel('m'); ax.legend(fontsize=8); ax.grid(axis='y', alpha=.3)
ax.set_title('Amplitudes (no son magnitudes equivalentes)')
plt.tight_layout(); plt.savefig(WORK / 'validation.png', dpi=130, bbox_inches='tight'); plt.show()

print(f'Estaciones comparables: {len(sim)} de {len(obs)}')
print('Excluidas:'); print(qc[qc.estado != 'ok'][['estacion', 'estado']].to_string(index=False))
print(f'\nArribo: diferencia simulado − observado entre {sim.dif_arribo_min.min():+.0f} y {sim.dif_arribo_min.max():+.0f} min')
print(f'Amplitud: η simulada / run-up observado entre {sim.razon_amplitud.min():.1f} y {sim.razon_amplitud.max():.1f}')

Estaciones comparables: 4 de 8
Excluidas:
         estacion                                                 estado
         San Juan sin mar de ≥ 10 m a < 10 km (el más cercano, a 101 km)
        El Charco  sin mar de ≥ 10 m a < 10 km (el más cercano, a 21 km)
     Buenaventura  sin mar de ≥ 10 m a < 10 km (el más cercano, a 19 km)
Salinas (Ecuador)                                      fuera del dominio

Arribo: diferencia simulado − observado entre -7 y +8 min
Amplitud: η simulada / run-up observado entre 1.2 y 2.0

Contraste con resultados publicados¶

Dos trabajos de la literatura simularon el escenario de 1979 en Tumaco con modelos no lineales con cálculo de inundación (Ávila-Guzmán et al. usan TUNAMI-N2) y distribuciones de deslizamiento de Poveda y Pulido (2019). Se citan aquí como referencia de orden de magnitud; su fuente, su malla y sus definiciones difieren de las de este cuaderno.

tum = sim.estacion.str.startswith("Tumaco")
pub = pd.DataFrame([
    ['Altura máxima del tsunami, isla de Tumaco, escenario 1979 (modelo de déficit de deslizamiento)', '2,30 m', 'Escobar et al. (2020), tabla 3'],
    ['Altura máxima del tsunami, isla de Tumaco, escenario 1979 (modelo de banda ancha)', '1,88 m', 'Escobar et al. (2020), tabla 3'],
    ['Llegada de la primera ola a Tumaco', '22 min', 'Ávila-Guzmán et al. (2026)'],
    ['Este cuaderno: η máxima en celdas de mar de 10-25 m, Tumaco', f'{sim.eta_max_m[tum].min():.1f}-{sim.eta_max_m[tum].max():.1f} m', 'LSWE lineal, 450 m'],
    ['Este cuaderno: primer ascenso > 0,1 m en Tumaco', f'{sim.t_arribo_sim_min[tum].min():.0f}-{sim.t_arribo_sim_min[tum].max():.0f} min', 'LSWE lineal, 450 m']],
    columns=['Magnitud', 'Valor', 'Fuente'])
pub

5. Sensibilidad (resumen del cuaderno 2)¶

sp = OUT / 'sensibilidad.csv'
if sp.exists():
    sens = pd.read_csv(sp, index_col=0)
    cols = [c for c in sens.columns if 'η máx' in c]
    display(sens[cols].describe().loc[['min', 'max']].T.round(2).rename(columns={'min': 'mínimo', 'max': 'máximo'}))
    display(sens[[c for c in sens.columns if 'arribo' in c]].describe().loc[['min', 'max']].T.round(1).rename(columns={'min': 'mínimo', 'max': 'máximo'}))
else:
    print('Ejecute la sección 6 del cuaderno 2 para generar output/sensibilidad.csv')

6. Animación de la propagación¶

Mapa con el oeste a la izquierda y el norte arriba; la estrella marca el epicentro y la falla se muestra con su proyección. η en metros (|η| < 0,1 m se deja transparente).

cf_t = None
fig, ax = plt.subplots(figsize=(8.2, 8))
fondo(ax)
ax.contour(np.linspace(ext[0], ext[1], nx), np.linspace(ext[3], ext[2], ny), H, levels=[200, 1000, 3000], colors='#8c99a6', linewidths=.4)
costa(ax)
agua = ax.imshow(np.full(H.shape, np.nan), cmap='RdBu_r', vmin=-3, vmax=3, extent=ext, alpha=.9)
plt.colorbar(agua, ax=ax, shrink=.7, label='η (m)')
ax.plot(ex / 1e3, ey / 1e3, '*', c='red', mec='k', ms=18, zorder=10, label='Epicentro 1979')
for nombre, (lon, lat) in {'Tumaco': (-78.76, 1.82), 'Esmeraldas': (-79.65, 0.97), 'Buenaventura': (-77.02, 3.88)}.items():
    x, y = TR.transform(lon, lat); ax.plot(x / 1e3, y / 1e3, 'ko', ms=4, zorder=9)
    ax.annotate(nombre, (x / 1e3, y / 1e3), textcoords='offset points', xytext=(5, 4), fontsize=8)
ax.legend(loc='upper left'); ejes_lonlat(ax)
ax.set_xlim(ext[0] + 5, ext[1] - 5); ax.set_ylim(ext[2] + 5, ext[3] - 5)
titulo = ax.set_title('')

def cuadro(k):
    e = snaps[k].astype('float32'); agua.set_array(np.where((np.abs(e) > 0.1) & (H >= 10), e, np.nan))
    titulo.set_text(f't = {snap_t[k]/60:.0f} min — Tsunami de Tumaco 1979 (LSWE)')
    return [agua, titulo]

anim = FuncAnimation(fig, cuadro, frames=len(snap_t), interval=250, blit=False)
gif = WORK / 'tsunami_propagation.gif'
anim.save(str(gif), writer=PillowWriter(fps=4), dpi=80)
plt.close(fig)
print(f'✓ Animación: {len(snap_t)} cuadros, t = 0–{snap_t[-1]/60:.0f} min, {gif.stat().st_size/1e6:.1f} MB')

✓ Animación: 31 cuadros, t = 0–60 min, 2.3 MB

from IPython.display import Image, display
display(Image(str(gif)))

7. Resumen¶

print('=' * 62); print('RESUMEN — TSUNAMI DE TUMACO 1979 (versión 2)'); print('=' * 62)
print(f'  Modelo:              LSWE, malla escalonada, {dx:.0f} m, dt = {dt:.0f} s')
print(f'  Fuente:              Okada (1985), Mw {8.3}, deslizamiento uniforme {float(R["slip"]):.1f} m (geometría supuesta)')
print(f'  η máx. (≥ 10 m prof.): {maxeta[H>=10].max():.2f} m | plataforma 10–200 m: {maxeta[zona].max():.2f} m')
print(f'  Estaciones:          {len(sim)} comparables de {len(obs)} (resto excluido, ver sección 4)')
print(f'  Arribo sim − obs:    {sim.dif_arribo_min.min():+.0f} a {sim.dif_arribo_min.max():+.0f} min (obs. sin verificar)')
print(f'  η sim / run-up obs:  {sim.razon_amplitud.min():.1f} a {sim.razon_amplitud.max():.1f}')
for f in ['max_eta_map.png', 'time_series.png', 'validation.png', 'tsunami_propagation.gif', 'fuente_okada.png']:
    p = WORK / f; print(f'  {"✓" if p.exists() else "✗"} {f}')

==============================================================
RESUMEN — TSUNAMI DE TUMACO 1979 (versión 2)
==============================================================
  Modelo:              LSWE, malla escalonada, 450 m, dt = 1 s
  Fuente:              Okada (1985), Mw 8.3, deslizamiento uniforme 7.4 m (geometría supuesta)
  η máx. (≥ 10 m prof.): 8.93 m | plataforma 10–200 m: 8.93 m
  Estaciones:          4 comparables de 8 (resto excluido, ver sección 4)
  Arribo sim − obs:    -7 a +8 min (obs. sin verificar)
  η sim / run-up obs:  1.2 a 2.0
  ✓ max_eta_map.png
  ✓ time_series.png
  ✓ validation.png
  ✓ tsunami_propagation.gif
  ✓ fuente_okada.png