513372 - Linux, Scripts y GMT

CLASE 11 - grillas

back

En esta clase, trabajamos con un set de datos que representa algo que se puede tomar en terreno; en este caso es datos tomados de la anomalía de Bouguer de gravedad:

bouguer.xyz

Contiene tres columnas: longitud (x), latitud (y), y anomalía de Bouguer en mgal (z); y representa una serie de mediciones del Altiplano. Se puede usar sort para ver los límites de latitud y longitud que contienen los datos. Este datos viene de la Universidad Libre de Berlin: http://www.fu-berlin.de/sfb/sfb267/results/data_catalogue/central_andean_data/gravitiy_data.html

Este set de datos no está contínuo, es decir que solamente mediciones están tomadas donde existen caminos o aceso al Altiplano. Podemos usar el comando surface para generar una superficie (es decir, una grilla contínua) que pasa por los puntos donde existen mediciones. Este comando requiere un factor de tensión entre 0 y 1 que impone restricciones sobre la curvatura de la superficie para evitar oscilaciones o fálsos máximos/mínimos en la superficie, y un límite de convergencia para definir cuando la superficie esta lista (es decir, que aparece bien a los datos, no hagamos más perturbaciones a la superficie). Probemos el siguiente en el terminal:

surface bouguer.xyz -Gbouguer.grd -R-71.5/-64/-29.6/-24.9 -I1m -T0.25 -C0.1 -V

Aquí tomamos el set de datos de gravidad bouguer.xyz y generamos una grilla (bouguer.grd) que pasa por sus puntos dentro de una cierta región que contiene los datos (-R), con una resolución de un minuto (es decir (1/60)°) para la grilla (-I1m), un factor de tensión de 0.25 (-T0.25), y un límite de convergencia de 0.1 mgal  con -C0.1 (es decir que cuando la iteración para generar la grilla esta cambiando los puntos un valor menor que 0.1, se para la iteración y tenemos nuestro resulto final).

Con este comando hemos generado una grilla bouguer.grd que representa los datos. Se puede revisar los contenidos de cualquier grilla con el comando grd2xyz para ver en el terminal los punto de la grilla. Revisen bouguer.grd usando este comando.

También queremos una paleta de colores con un rango suficiente para cubrir los valores de la grilla. Podemos crear uno pero en este caso esta más facil usar una paleta ya creada para el gráfico. Usamos la paleta de colores GMT_red2green.cpt (que debe estar incluido en la instalación de GMT) para mostrar eso. Si miren la paleta, se pueden ver que la paleta varía de -1 (rojo) a +1 (verde) mientras que los datos tiene valores en los cientos de mgal; entonces queremos usar el comando grd2cpt para generar una paleta apropriada para los datos.

grd2cpt bouguer.grd -C/home/matt/GMT/GMT4.5.2/share/cpt/GMT_red2green.cpt -T= > bouguer.cpt

Con el comando de arriba estoy tomando una paleta de colores y cambiando su rango para que se aplica a los valores de una cierta grilla. (Obviamente, hay que cambiar el camino at .cpt para que da la ubicación de sus paletas). El -T= está una opción para que la paleta de colores generada (bouguer.cpt) sea simétrica, entonces una anomalía de Bouguer de cero será puesta en blanco, anomalías negativas en rojo y anomalías positivas en verde. Revisen la paleta y la grilla para los datos y noten que con esta paleta podemos graficar los datos con un rango de colores aceptables.

Ahora estamos listo para generar/modificar un script para graficar los datos. Bajen el script de abajo que maneja el set de datos:

ejemplo_bouguer.sh

Este script es una modificación de los scripts anteriores. Empiece definiendo unos variables (hay que cambiar "palette" y "topography" para dar los caminos a sus archivos), después grafíca el basemap y la topografía, con iluminación, usando psbasemap, grdgradient y grdimage. Los próximos comandos son los de surface y grd2cpt para definir la superficie que representa los datos y la paleta con que vamos a graficarlo.

Después vienen los próximos dos líneas que usan el comando grdsample:

grdsample ${region}.int -Gtopo.int -R-71.5/-64/-29.6/-24.9 -I1m -V -F
grdsample bouguer.grd -Gbouguer2.grd -R-71.5/-64/-29.6/-24.9 -I1m -V -F

¿Que estoy haciendo aquí? Bueno, siempre la anomalía de gravedad está asociada con la topografía (piensen en los cursos de Geofísica Tierra Sólida o Geodesia; pero ellos que no han hecho estas cursos noten que la interpretación de esta imagen no es necesaria en el curso de GMT) entonces quiero mostrar la anomalía de gravedad y la topografía en el mismo imagen. Para esto, puedo usar el comando grdsample para forzar (-F) los puntos de las grilla de Bouguer (bouguer.grd) y la grilla de la iluminación de topografía (${region}.int) a exáctamente los mismos puntos, cada grilla en la misma región (-R) con la misma resolución de 1 minuto (-I). Estos comandos creen las grillas bouguer2.grd y topo.int respectivamente. Ahora puedo graficar la grilla de Bouguer en colores de la paleta bouguer.cpt pero con iluminación asociada con la topografía:

grdimage bouguer2.grd -Cbouguer.cpt -B${ticks} -J${projection} -Itopo.int -R${bounds} ${portrait} ${verbose} -O -K >> ${psfile}

Después vienen 5 líneas para

(i) graficar los puntos de datos (bouguer.xyz) como círculos en el mapa para ver donde la grilla está confiable.
(ii) poner contornos a la superficie.
(iii) y (iv) graficar la costa, y los lagos en un cierto azul definido por ${lakes}.
(v) poner una escala al lado del gráfico para la paleta de colores usada.

more bouguer.xyz | psxy -J${projection} -R${bounds} ${verbose} -Sc0.1c -Cbouguer.cpt -W3 -O -K >> ${psfile}
grdcontour bouguer2.grd -B${ticks} -J${projection} -C50 -A100+k0/0/0+g255/255/255 -R${bounds} ${portrait} ${verbose} -W0.5p/0/0/0 -Gd5c -O -K >> ${psfile}
pscoast -B${ticks} -J${projection} -R${bounds} ${portrait} ${verbose} -D${resolution} -W${coastline} ${borders} ${rivers}  -O -K >> ${psfile}
pscoast -B${ticks} -J${projection} -R${bounds} ${portrait} ${verbose} ${lakes} -D${resolution} -W${coastline} ${borders} ${rivers}  -O -K >> ${psfile}
psscale  -D16c/2.5c/5c/0.5c -Cbouguer.cpt ${verbose} -Ba100f20:"(mGal)": -O -K >> ${psfile}

Bueno, hemos llegado a la última parte del script. En esta estamos graficando un perfíl de la grilla. Primeramente usamos el comando project para generar un archivo track.xy que consiste de puntos intermedios entre dos ubicaciones de longitud/latitud con un ciero espaciado entre los puntos:

project -C-73/-28 -E-62/-28 -G0.01 > track.xy

En este caso estamos definiendo un perfíl a latitud 28 sur que cruce los Andes y pasa por los puntos de datos. El espaciado entre los puntos es 0.01 grados (revisen el archivo track.xy). Ahora queremos tomar los datos de las grillas de Bouguer y de topografía a los puntos a lo largo del perfíl; por eso usamos grdtrack:

grdtrack track.xy -Gbouguer2.grd > bouguertrack.xydz
grdtrack track.xy -G${region}.grd > topotrack.xydz

Estamos generando archivos con 4 columnas: longitud, latitud, distancia a lo largo del perfíl y el valor interpolado de la grilla a estos puntos. Noten que el perfíl tiene un largo de 9.709 grados o 1078 kilómetros.

Para el gráfico del perfíl, queremos definir ejes entre 0 a 1078 (kilómetros) y -500 a 500 (mgal). Los ejes deben estar abajo del mapa, entonces usamos -X0 -Y-4, con esta ubicación relativa a la ubicación del mapa original, es decir 4 centímetros abajo de él. También queremos que el ancho del perfíl es lo mismo que el mapa, entonces uasmos -JX14/2 para que el perfíl es 14 cm largo y 2 cm alto.

psbasemap -R0/1078/-500/500 -JX14/2 -X0 -Y-4 -Ba200f100 -P -O -K -V >>  ${psfile}

Ahora queremos graficar la anomalía de Bouguer a lo largo del perfíl. Tengo que multiplicar la tercera columna por 111 para convertir entre grados y kilómetros, y estoy graficando una línea con psxy:

more bouguertrack.xydz | awk '{print $3*111, $4}' | psxy -R0/1078/-500/500 -JX14/2 -B -P -W2/0/0/0t20_10:10 -V -O -K >> ${psfile}

Noten la línea (-W2/0/0/0t20_10:10): estoy tomando una línea con un ancho de 2 puntos (-W2), de color negro (0/0/0), con rayas de 20 puntos de largo t20 separadas con espacios de 10 puntos (_10), y la primera raya tiene un largo de 10 puntos solamente (:10). Prueben diferentes combinaciones de valores acá para generar diferentes líneas punteadas (prueben t20_10_5_10 también - ¿qué hace?).

El perfíl de topografía está de una línea sólida roja. Note el escalamiento de los valores de topografía ($4/20) para que podemos usar los mismos ejes. Con eso la topografía en kilómetros tiene una exageración de un factor de 50.

more topotrack.xydz | awk '{print $3*111, $4/20}' | psxy -R0/1078/-500/500 -JX14/2 -B -P -W3/255/0/0 -V -O -K >> ${psfile}

Finalmente ponemos una línea entre los puntos (0,0) y (1078,0) para que el perfíl se ve un poco mejor:

psxy << END -R0/1078/-500/500 -JX14/2 -B -P -W1/0/0/0 -O >> ${psfile}
0 0
1078 0
END

En este caso, los valores entre el comando psxy << END (literalmente, correr psxy usando las siguientes líneas de texto hasta que encuentras una línea que dice "END") y el END representan el equivalente de un archivo .xy.

Este script debe introducir una manera de poner datos tomados en terreno a un mapa que representa la región. Debe generar el imagen que sigue. El codigo completo del script viene después del imagen.


#!/bin/bash

region="ejemplo_bouguer"
bounds="-73/-62/-32/-22"

#origin:
    xshift="2.0"
    yshift="10.0"

projection="M14.0"                 #projection
palette="/home/matt/GMT/GMT4.5.2/share/cpt/GMT_ocean.cpt"   # colour palette
topography="/home/matt/GMT/TOPO_FILES/SRTM30/w100s10.Bathymetry.srtm.grd"    #topography file
ticks="a5f1g5SEWN"
illaz="225"     # illumination azimuth
portrait="-P"    # set this to be blank for the default landscape mode
verbose="-V"    # set this to be blank to suppress output to stderr

# coastline and border information --- see pscoast manpage for details
coastline="1.00p/0/0/0"
resolution="f"    # choose from (f)ull, (h)igh, (i)ntermediate,
        # (l)ow, and (c)rude
borders="-N1/1.00p/0/0/0"
#rivers="-I1/1.00p/0/0/255"

 wet="-S0/0/255"
   
 land="-G0/255/0"
       
lakes="-C50/100/200 -A0/2/4"

########################################################################################################

# output file
psfile="${region}.ps"

# make a basemap (blank b/g), insert grid image and plot coastlines
psbasemap -B${ticks} -J${projection} -R${bounds} -X${xshift} -Y${yshift} ${portrait} ${verbose} -K > ${psfile}

#grdcut input_file.grd −Goutput_file.grd −Rwest/east/south/north[r] [ −V ] [ −f[i|o]colinfo ]
#cuts the grid to the specific region
grdcut ${topography} -G${region}.grd -R${bounds} -V

    # illuminate the grid file from specified azimuth to produce an
    # intensity (.int) file
    grdgradient ${region}.grd -G${region}.int -A${illaz} -Nt -M
    # insert topography
    grdimage ${region}.grd -C${palette} -I${region}.int -B${ticks} -J${projection} -R${bounds} ${portrait} ${verbose} -O -K >> ${psfile}
   
#############################################################################################

surface bouguer.xyz -Gbouguer.grd -R-71.5/-64/-29.6/-24.9 -I1m -T0.25 -C0.1 -V
grd2cpt bouguer.grd -C/home/matt/GMT/GMT4.5.2/share/cpt/GMT_red2green.cpt -T= > bouguer.cpt
#makecpt -C/home/matt/GMT/GMT4.5.2/share/cpt/GMT_red2green.cpt -T-500/500/10 > bouguer.cpt
grdsample ${region}.int -Gtopo.int -R-71.5/-64/-29.6/-24.9 -I1m -V -F
grdsample bouguer.grd -Gbouguer2.grd -R-71.5/-64/-29.6/-24.9 -I1m -V -F
    grdimage bouguer2.grd -Cbouguer.cpt -B${ticks} -J${projection} -Itopo.int -R${bounds} ${portrait} ${verbose} -O -K >> ${psfile}



more bouguer.xyz | psxy -J${projection} -R${bounds} ${verbose} -Sc0.1c -Cbouguer.cpt -W3 -O -K >> ${psfile}

grdcontour bouguer2.grd -B${ticks} -J${projection} -C50 -A100+k0/0/0+g255/255/255 -R${bounds} ${portrait} ${verbose} -W0.5p/0/0/0 -Gd5c -O -K >> ${psfile}


pscoast -B${ticks} -J${projection} -R${bounds} ${portrait} ${verbose} -D${resolution} -W${coastline} ${borders} ${rivers}  -O -K >> ${psfile}
pscoast -B${ticks} -J${projection} -R${bounds} ${portrait} ${verbose} ${lakes} -D${resolution} -W${coastline} ${borders} ${rivers}  -O -K >> ${psfile}

    psscale  -D16c/2.5c/5c/0.5c -Cbouguer.cpt ${verbose} -Ba100f20:"(mGal)": -O -K >> ${psfile}

project -C-73/-28 -E-62/-28 -G0.01 > track.xy
grdtrack track.xy -Gbouguer2.grd > bouguertrack.xydz
grdtrack track.xy -G${region}.grd > topotrack.xydz
psbasemap -R0/1078/-500/500 -JX14/2 -X0 -Y-4 -Ba200f100 -P -O -K -V >>  ${psfile}
more bouguertrack.xydz | awk '{print $3*111, $4}' | psxy -R0/1078/-500/500 -JX14/2 -B -P -W2/0/0/0t20_10:10 -V -O -K >> ${psfile}
more topotrack.xydz | awk '{print $3*111, $4/20}' | psxy -R0/1078/-500/500 -JX14/2 -B -P -W3/255/0/0 -V -O -K >> ${psfile}
psxy << END -R0/1078/-500/500 -JX14/2 -B -P -W1/0/0/0 -O >> ${psfile}
0 0
1078 0
END


back