Visualización Geospacial Avanzada con PyGMT

Desarrollo del proyecto final

En este notebook construyes, paso a paso, el proyecto final del curso: un análisis visual de la sismicidad del sur de México en septiembre de 2017, el mes del sismo de Chiapas (magnitud 8.1, el día 8) y del sismo de Puebla-Morelos (magnitud 7.1, el día 19). Combinarás un mapa de relieve con los sismos, un mapa de ubicación, un perfil que muestra cómo se hunde la placa de Cocos bajo México y un resumen numérico. Los datos son reales: el catálogo usgs_quakes de PyGMT y el relieve SRTM15 del servidor de GMT.

Paso 1: preparar los datos

Cargas el catálogo, filtras la región de estudio y calculas los números que vas a citar en las conclusiones. Separar los datos de la visualización permite rehacer un mapa sin repetir el filtrado:

import pygmt

region = [-106, -88, 12, 22]
sismos = pygmt.datasets.load_sample_data("usgs_quakes")
sur = sismos[sismos.longitude.between(region[0], region[1])
             & sismos.latitude.between(region[2], region[3])].copy()
print(f"Sismos en la región: {len(sur)}")
print(f"Profundidad media: {sur.depth.mean():.1f} km")
print(sur.nlargest(3, "mag")[["time", "mag", "depth", "place"]])

Paso 2: el mapa principal

El mapa principal reúne las técnicas del curso: relieve con sombreado, sismos con color por profundidad y tamaño por magnitud, y una escala de color con unidades. Usas dos paletas: una para el relieve y otra para la profundidad de los sismos. Por eso la segunda makecpt va justo antes de dibujar los sismos.

grid = pygmt.datasets.load_earth_relief(resolution="05m", region=region)
sombra = pygmt.grdgradient(grid=grid, azimuth=315, normalize="e0.4")
fig = pygmt.Figure()
pygmt.makecpt(cmap="gray", series=[-6000, 5000])
fig.grdimage(grid=grid, shading=sombra, projection="M16c", cmap=True,
             frame=["WSne+tSismicidad del sur de México, septiembre de 2017", "af"])
fig.coast(shorelines="0.4p,black", borders="1/0.5p,gray20")
pygmt.makecpt(cmap="lajolla", series=[0, 150])
fig.plot(x=sur.longitude, y=sur.latitude, size=0.02 * 2 ** (sur.mag - 3),
         fill=sur.depth, cmap=True, style="cc", pen="0.2p,black")
fig.colorbar(frame="af+lProfundidad (km)")
fig.savefig("mapa_sismicidad.png")

El relieve va en grises a propósito: así el color queda reservado para la variable que importa, la profundidad. Después agregas el mapa de ubicación con fig.inset y etiquetas para los dos sismos principales con fig.text.

Paso 3: el perfil de subducción

Un mapa muestra dónde ocurren los sismos; un perfil muestra a qué profundidad. Trazas una línea perpendicular a la costa de Chiapas, que empieza en el océano (punto A) y termina tierra adentro (punto B), y proyectas sobre ella los sismos cercanos con pygmt.project: cada sismo recibe su distancia a lo largo del perfil. Con unit=True las distancias salen en kilómetros width limita los sismos a una franja de 150 km a cada lado de la línea y length="w" descarta los que quedan más allá de A o de B (sin él, project usa la línea infinita):

A, B = [-95.5, 13.5], [-92.5, 18.0]
datos = sur[["longitude", "latitude", "depth", "mag"]]
perfil = pygmt.project(data=datos, center=A, endpoint=B, convention="pz",
                       unit=True, width=[-150, 150], length="w")
perfil.columns = ["distancia", "profundidad", "magnitud"]
print(len(perfil), "sismos en la franja del perfil")

fig = pygmt.Figure()
fig.plot(x=perfil.distancia, y=perfil.profundidad, region=[0, 600, 0, 150],
         projection="X15c/-6c", style="c0.15c", fill="red3", pen="0.2p",
         frame=["WSne", "xaf+lDistancia a lo largo del perfil (km)", "yaf+lProfundidad (km)"])
fig.savefig("perfil_subduccion.png")

El eje vertical va invertido (-6c) para que la profundidad crezca hacia abajo, como en un corte geológico. Arriba del perfil puedes agregar la topografía de la misma línea: pygmt.project(..., generate=5) crea puntos cada 5 km y pygmt.grdtrack lee la elevación de la grilla en cada uno.

Paso 4: resumen y conclusiones

Cierra el notebook con una tabla corta y tres o cuatro frases. La tabla cuenta los sismos por rango de profundidad (someros, de menos de 70 km; intermedios, de 70 a 300 km) y da la magnitud máxima de cada grupo. Las conclusiones responden preguntas concretas: dónde se concentran los sismos, cómo cambia la profundidad al alejarse de la costa y qué relación tiene eso con la placa que se hunde.

import pandas as pd

sur["tipo"] = pd.cut(sur.depth, bins=[0, 70, 300],
                     labels=["somero", "intermedio"], include_lowest=True)
resumen = sur.groupby("tipo", observed=True).agg(
    sismos=("mag", "size"), mag_max=("mag", "max"), prof_media=("depth", "mean"))
print(resumen.round(1))

Trampas comunes

• Usar la misma paleta para relieve y sismos. El lector no sabe qué color pertenece a qué variable.
• Un perfil en dirección arbitraria. Para ver la placa que se hunde, el perfil debe cruzar la trinchera de forma perpendicular.
• Una franja demasiado ancha en project. Mezcla sismos de zonas distintas y borra la forma de la placa.
• Olvidar length="w". Sin él, entran sismos alineados con el perfil pero a miles de kilómetros.
• Confiar en la región «recordada». PyGMT reutiliza la última region de la sesión; pásala explícitamente al dibujar una grilla nueva.
• Conclusiones sin números. Cada afirmación debe apoyarse en un conteo, una media o una figura del notebook.