17.1. Mapa de Calor de Sismicidad

Hacer mapa de calor (heatmap) de las zonas con mayor sismicidad del sur de Sudamérica y del Mar de Scotia.

Introducción:
En el mapa de sismicidad se observa que los sismos cubren gran parte de la región cordillerana y del Mar de Scotia. Sin embargo, debido a la gran cantidad no es posible definir cuales son las zonas de mayor y menor sismicidad. Una forma de observar esto es mediante un mapa de calor (heatmap) que permita ver cuantos sismos ocurrieron por áreas. Para esto definimos bloques de 10 x 10 km (100 km2), contamos la cantidad de sismos ocurridos en cada uno y creamos una grilla que de dibuja sobre un mapa base.

Datos:

• Topografía: GMRTv3_3_Arg.grd. Grilla GMRT (ver anexo 4).
• Sismos: Distintos archivos de texto plano separados por comas (query*.csv link archivos comprimidos) descargados del USGS (ver Base de datos). Los archivos tienen 1 linea de encabezado (-h1). Los datos de ubicación (long, lat) están en las columnas 3 y 2 respectivamente (-i2,1).

Script

0. Sismicidad: El mapa se graficó con los siguientes comandos:

gmt makecpt > "temp_seis.cpt" -Crainbow -T0/700 -Z -I
gmt psxy -R -J -O -K "E:\Facultad\Datos_Geofisicos\Sismicidad\USGS\query_*.csv" -h1 -i2,1,3   >> %OUT% -Sc0.05 -C"temp_seis"

1. Mapa Base: Los siguientes comandos crean el mapa base con límites de países, provincias y línea de costa. Ver 5.2 y 5.5 para más detalles.

gmt makecpt -T-11000,0,9000 -Cdodgerblue2,white > %color%
gmt grdgradient %CUT% -A270 -G%SHADOW% -Ne0.5
gmt grdimage -R -J -O -K %CUT% -C%color% -I%SHADOW% >> %OUT%
gmt pscoast -R -J -O -K >> %OUT% -Df -N1/0.30
gmt pscoast -R -J -O -K >> %OUT% -Df -N2/0.25,-.
gmt pscoast -R -J -O -K -Df -W1/faint >> %OUT%
gmt psbasemap -R -J -O -K >> %OUT% -Ln0.11/0.075+c-32:00+w800k+f+l -F+gwhite+p+i+s
gmt psbasemap -R -J -O -K >> %OUT% -Ba8f4/a4f2

2. Resolución:  Se define la variable %res% como 10 km (res=10k). Se utiliza para crear los bloques (punto 3) y definir la resolución de la grilla (punto 4). Las dimensiones de los bloques será diferente si se utilizan unidades planas (km, m) o coordenadas (grados, minutos, segundos).

SET    res=10k

3. Procesar Sismos: blockmean lee los datos de sismos (query*.csv), define los bloques y cuenta cuantos sismos ocurrieron en cada uno. Se utilizan los siguientes argumentos:

• -I%res%: define el tamaño de los bloques, según la variable %res%. En este caso es de 10 x 10 km (100 km2). 
• -C: establece que los datos del archivo de salida estén en el centro del bloque (en lugar de la posición promedio).
• -h y -i: respectivamente definen cuantas filas conforman el encabezado (o header; -h) y en que columnas están los datos (-i).
• -Sn: cuenta la cantidad de sismos ocurridos en cada bloque.

gmt blockmean -R "E:\Facultad\Datos_Geofisicos\Sismicidad\USGS\query_*.csv" -Sn -C -h1 -i2,1 > "temp_Heatmap.xyz" -I%res%

El archivo de salida ("temp_Heatmap.xyx") contiene una primera línea (comentada) con el comando utilizado para crear el archivo y luego 3 columnas con la longitud y latitud del centro de cada bloque y la cantidad de sismos que en cada uno.

4. Crear Grilla: xyz2grd crea una grilla a partir de "temp_Heatmap.xyz".

• -I%res%: define la resolución. Utilizamos la variable %res% (10 km).
• -G%CUT2%: nombre del archivo de salida.

gmt xyz2grd -R -G%CUT2% "temp_Heatmap.xyz" -I%res% 


5. Analizar datos: Debido a la distribución asimétrica de los datos es necesario analizar su distribución (para luego crear una paleta de color para el mapa de color; CPT). Hay 2 métodos para realizarlo:

I. Histograma: pshistogram analiza los datos de "temp_Heatmap.xyz",

• -W5: define el ancho de clase.
• -Io: muestra los datos en la pantalla. 
• -i2: selecciona la columna con los datos (la 3a), 
• -Z0: cuenta la cantidad de eventos por intervalo de clase

gmt pshistogram "temp_Heatmap.xyz" -W5 -Io -i2 -Z0


II. grdinfo: -T muestra el intervalo de los datos se extienden entre 0 y 485 (para un intervalo de 5).

• -T+a: permite obtener el rango de valores los datos para el 99% (+a1), 97,5 (+a2.5) y 95% (+a5).

grdinfo %CUT2% -T5
grdinfo %CUT2% -T5+a1
grdinfo %CUT2% -T5+a2.5
grdinfo %CUT2% -T5+a5

6. Paleta de Color: grd2cpt crea una paleta de color que se utilizara para pintar el mapa de calor. Se guarda como %color2%.

• -Z: paleta continua. 
• -C: nombre de la CPT maestra (hot).
• -L: intervalo de la CPT de 0 a 45 (97,5 % de los datos).
• -Di: utiliza el color para el valor máximo (45) para los valores mayores.

gmt grd2cpt %CUT2% > %color2% -Z -Di -Chot -L0/45


7. Dibujar Mapa de Calor: grdimage dibuja el mapa de calor.

• -C%color2%: usa la CPT. 
• -Q: bloques con valor 0 no son pintados. 
• -t25: establece la transparencia de la capa para las áreas pintadas (25%).

gmt grdimage -R -J -O -K %CUT2% -C%color2% >> %OUT% -Q -t25

7. Escala de colores: psscale dibuja la escala con los siguientes argumentos:

• -C%color2%: Utiliza la paleta %color2%.
• -DjRT+o1.7c/0.7+w9/0.618c+ef: Define la posición y dimensiones de la escala.
• -Baf+l: Dibuja las divisiones de la escala y la leyenda (texto entre "").
• -F: Dibuja un recuadro con fondo blanco (+g), borde principal (+p) e interno (+i) y con sombra (+s).

gmt psscale -R -J -O -K -C%color2% >> %OUT% -DjRT+o1.7c/0.7+w9/0.618c+ef -Baf+l"Cantidad de Sismos cada 100 km@+2@+"  -F+gwhite+p+i+s