513372 - Linux, Scripts y GMT

CLASE 10 - más símbolos

back

10.0

El script usado en esta clase es ejemplo_mapa5.sh. Es una modificación del script anterior.

10.1 símbolos desde un archivo de texto: volcanes

Podemos usar psxy para graficar símbolos desde un archivo de datos, como southern_Chile_and_Argentina.txt. En la misma manera que antes, el comando psxy en GMT toma coordenadas de longitud, latitud y grafica un símbolo en esta posición, también podemos mandar varias coordenadas al comando al mismo tiempo con un pipe:

more southern_Chile_and_Argentina.txt  | awk '{print $6, $5}' | grep -v lon | psxy -J${projection} -R${bounds} $verbose} -St0.5c -G255/0/0 -W2 -O -K >> ${psfile}

Aquí tomamos el archivo de datos, aislamos la longitud y latitud de los volcanes usando una combinación de awk y grep, y los graficamos con triángulos de 0.5 cm (-St0.5c), color rojo (-G255/0/0), con un borde negro de dos puntos (-W2). Esta línea de codigo va dentro de un script, como antes.

10.2 símbolos con diferentes colores y tamaños

Aquí trabajamos con el archivo de datos global_seismicity_feb27-apr19_2010.txt, de la sismicidad global de las próximas semanas después del terremoto en Chile 2010. 
Podemos aislar las columnas de longitud, latitud, profundidad y magnitud, si existe una magnitud, usando el siguiente comando:

    more global_seismicity_feb27-apr19_2010.txt | sed 's/,//g' | awk '{if ( $10 > 0 ) print $5, $4, $6, $10}'

Si esto es combinado con psxy, la tercera columna da el color del símbolo, y la cuarta su tamaño. Por eso, necesitamos una paleta de colores para los terremotos a diferentes profundidades. Podemos generar una , o podemos modificar una. Por ejemplo, empezamos con GMT_copper.cpt, que da diferentes tonos de cobre entre valores de cero y uno (revisa el .cpt), y generamos una nueva paleta con el comando

    makecpt -C/home/matt/GMT/GMT4.5.2/share/cpt/GMT_copper.cpt -T0/100/10 -I > depths.cpt


El -T especifica la escala del .cpt (0 a 100), y el ancho de sus bandas (10). El -I especifica que tomamos los colores de GMT_copper.cpt al revés. Pruebe el comando sin el -I para ver la diferencia. Esta paleta de colores no es perfecta, pero se nota que la paleta original (GMT_copper.cpt) es muy simple, y es fácil generar otras paletas entre 0 y 1 con diferentes colores, y después modificar el rango con la opción -T en makecpt.

Ahora estamos listos para graficar los terremotos, con

more global_seismicity_feb27-apr19_2010.txt | sed 's/,//g' | awk '{if ( $10 > 0 ) print $5, $4, $6, ($10*$10*$10/500)}' | sort -k3 -r | psxy -J${projection} -R${bounds} ${verbose} -Sc -Cdepths.cpt -W2 -O -K >> ${psfile}
• Note que uso  sort -k3 -r para ordenar los terremotos por su profundidad. Con eso grafico los mas profundos primeros.
• Note que este catalogo era preliminar, entonces las profundidades de los terremotos no son muy precisas.
• El color del símbolo depende de su profundidad (-Cdepths.cpt). Puedo graficar la escala con:
psscale  -D16c/2.5c/5c/0.5c -Cdepths.cpt -Ba50f10/a50f10:"Depth (km)": ${verbose} -O -K >> ${psfile}

• Los símbolos van a ser circulos (-Sc); si no especifico un tamaño aquí, va a tomar el tamaño de la cuarta columna de los datos. En nuestro caso usamos 
$10*$10*$10/500 para que un terremoto de magnitud 8 tenga un tamaño de ~1 cm, y magnitud 4 tenga tamaño ~0.12 cm. Por supuesto, puedo poner una escala al lado del mapa con:
echo "-69.5 -34.5 35 1.024" | psxy -J${projection} -R${bounds} ${verbose} -Sc -Cdepths.cpt -N -W2 -O -K >> ${psfile}
echo "-69.2 -34.5 12 0 1 LM 8.0" | pstext -J${projection} -R${bounds} ${verbose} -G0 -N -W255 -O -K >> ${psfile}
echo "-69.5 -34.8 35 .432" | psxy -J${projection} -R${bounds} ${verbose} -Sc -Cdepths.cpt -N -W2 -O -K >> ${psfile}
echo "-69.2 -34.8 12 0 1 LM 6.0" | pstext -J${projection} -R${bounds} ${verbose} -G0 -N -W255 -O -K >> ${psfile}
echo "-69.5 -35.1 35 .128" | psxy -J${projection} -R${bounds} ${verbose} -Sc -Cdepths.cpt -N -W2 -O -K >> ${psfile}
echo "-69.2 -35.1 12 0 1 LM 4.0" | pstext -J${projection} -R${bounds} ${verbose} -G0 -N -W255 -O -K >> ${psfile}
Donde grafico símbolos de un cierto tamaño al lado del mapa, y la magnitud que coresponde. El -N en el comando  permite graficar cosas afuera del borde del mapa.

En conclusión, podemos tomar un set de datos y graficarlo con símbolos cuyos tamaños y colores dependen de los datos. En este caso, los datos son tomdos de IRIS: http://www.iris.edu/SeismiQuery/sq-events.htm. Ahora, tomando los resultados de la sismicidad global durante esta época, usando el catálogo "PREFERRED", va a tener las profundidades, y ubicaciones, de los sismos con mayor precisión.

Este ejemplo usa sismicidad, pero se  puede extender a cualquier set de datos que tengan latitud, longitud y ciertos valores que se necesitan graficar como símbolos.

10.3 mecanismos focales

El comando psmeca es para graficar mecanismos focales de terremotos. Por ejemplo, podemos buscar el mecanismo focal del terremoto del 2010 Chile en el Harvard CMT catalog, y elegir que su formato esta en
"GMT psmeca input". Para graficarlo, usamos el comando:

echo "-73.15 -35.98 23 1.04 -0.04 -1.00 0.30 -1.52 -0.12 29 X Y 201002270634A" | psmeca ${portrait} -J${projection} -R${bounds} ${verbose} -L -Sm1.0c -G255/0/0 -O -K >> ${psfile}

Note que si queremos graficar varios mecanismos focales, podemos guardar ellos en un archivo .txt y correr el comando :
more mecanismos.txt
| psmeca ${portrait} -J${projection} -R${bounds} ${verbose} -L -Sm1.0c -G255/0/0 -O -K >> ${psfile}

Esta sección es especifica a sismología, no a geofísica en general, así que voy a dejarlo hasta aquí. Los sismólogos entre ustedes pueden leer mas sobre los mecanismos focales y el manual para psmeca para más conocimiento.

10.4 indent maps

Siempre, cuando se grafica una región, es necesario poner un mapa al lado, con esta región marcada en un mapa general del país o continente. Para esto, podemos graficar la costa de una  región diferente en la esquina de nuestro mapa

bounds2=-90/-50/-50/-25
pscoast -JM3c -W1p -R${bounds2} -Y11.25 -X0.25 -Dh -G128 -S255 -O -K >> ${psfile}

Note que el -Y y -X son elegidos para ubicar este mapa pequeño en la esquina.

Para marcar el borde de la región original, podemos usar psxy de nuevo:

echo ${bounds} | sed 's/\// /g' | awk '{printf"%s %s\n %s %s\n %s %s\n %s %s\n %s %s\n", $1, $3, $2, $3, $2, $4, $1, $4, $1, $3}' | psxy -R${bounds2} -JM3c -A -W0.5p -O >> ${psfile}

Si no es obvio por que se usa un printf, pruebe los siguientes tres comandos en un terminal:

bounds="-75.5/-70/-38.5/-34"

echo ${bounds} | sed 's/\// /g'

echo ${bounds} | sed 's/\// /g' | awk '{printf"%s %s\n %s %s\n %s %s\n %s %s\n %s %s\n", $1, $3, $2, $3, $2, $4, $1, $4, $1, $3}'

10.5

Acá esta la imagen final usando el script generado en esta clase. No esta listo, pero espero que ahora ustedes pueden tomar un set de datos, que consiste de latitud, longitud y unos valores, y representarlo graficamente.



back