Grillas

En esta sección trabajamos con el set de datos bouguer.xyz que representa mediciones de terreno, en este caso, datos tomados de la anomalía de Bouguer de gravedad.
El archivo 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. Estos datos vienen de la Universidad Libre de Berlin.
Este set de datos no es continuo, es decir, las mediciones están tomadas donde existen caminos o aceso al Altiplano. Podemos usar el comando surface para generar una superficie, es decir, una grilla continua que pase 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 falsos máximos/mínimos en la superficie, y un límite de convergencia para definir cuando la superficie esta lista, esto es, que tiene buen parecer a los datos. Probemos el siguiente comando 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 gravedad 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)$^{o}$ 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, esto significa que cuando la iteración para generar la grilla esta cambiando los puntos un valor menor que 0.1, se detiene la iteración y tenemos nuestro resultado 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 puntos 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 una pero en este caso es más fácil usar una paleta ya creada para el gráfico. Usamos la paleta de colores GMT_red2green.cpt para mostrar eso. Si miran la paleta, se puede ver que la paleta varía de -1 (rojo) a +1 (verde), mientras que los datos tienen valores en los cientos de mgal, por lo cual usaramos el comando grd2cpt para generar una paleta apropriada para los datos.

grd2cpt bouguer.grd -C/usr/local/GMT4.5.8/share/cpt/GMT_red2green.cpt -T= > bouguer.cpt

Con el comando de arriba se toma una paleta de colores y se cambia su rango para que se aplique a los valores de una cierta grilla. Recuerden cambiar el camino .cpt a la ubicación de sus paletas. El -T= es una opción para que la paleta de colores generada bouguer.cpt sea simétrica, así la anomalía de Bouguer de valor 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. Descargen el script ejemplo_bouguer.sh que maneja el set de datos.
Este script es una modificación de los scripts anteriores. Empieza definiendo unas variables (hay que cambiar "palette" y "topography" para dar los caminos a sus archivos), después grafica 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 graficarla.
Después vienen las próximas 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 aqui? - Bueno, siempre la anomalía de gravedad está asociada con la topografía, recuerden los cursos de Geofísica de la Tierra Sólida o Geodesia. Aquellos que no han hecho estos cursos noten que la interpretación de la imagen no es necesaria en el curso de GMT, por los cual mostraré la anomalía de gravedad y la topografía en la misma imagen. Para esto, puedo usar el comando grdsample para forzar -F los puntos de la grilla de Bouguer bouguer.grd y la grilla de la iluminación topográfica 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 crean 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${proyeccion} -Itopo.int -R${limites} \
${portrait} ${verbose} -O -K >> ${psfile}

Las cinco líneas siguientes tienen la función de:

1.

Graficar los puntos de datos bouguer.xyz como círculos en el mapa para ver donde la grilla es confiable.

2.

Poner contornos a la superficie.

3.

Graficar la costa.

4.

Graficar los lagos en un cierto azul definido por ${lagos}.

5.

Poner una escala al lado del gráfico para la paleta de colores usada.

more bouguer.xyz | psxy -J${proyeccion} -R${limites} ${verbose} -Sc0.1c \
-Cbouguer.cpt -W3 -O -K >> ${psfile}
grdcontour bouguer2.grd -B${ticks} -J${proyeccion} -C50 -A100+k0/0/0+g255/255/255 \
-R${limites} ${portrait} ${verbose} -W0.5p/0/0/0 -Gd5c -O -K >> ${psfile}
pscoast -B${ticks} -J${proyeccion} -R${limites} ${portrait} ${verbose} \
-D${resolucion} -W${linea_costa} ${bordes} ${rios}  -O -K >> ${psfile}
pscoast -B${ticks} -J${proyeccion} -R${limites} ${portrait} ${verbose} ${lagos} \
-D${resolucion} -W${linea_costa} ${bordes} ${rios} -O -K >> ${psfile}
psscale  -D15.7c/2.5c/5c/0.5c -Cbouguer.cpt ${verbose} -Ba100f20:"(mGal)": \
-O -K >> ${psfile}

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 cierto espaciado entre los puntos.

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

En este caso estamos definiendo un perfil a latitud 28$^\circ$ sur que cruce los Andes y pase 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 perfil; pora 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 perfil y el valor interpolado de la grilla en esos puntos. Noten que el perfil tiene un largo de 9.709 grados o 1078 kilómetros. Para el gráfico del perfil, 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 perfil sea el mismo al mapa, para eso usamos -JX14/2 para que el perfil sea 14 cm en largo y 2 cm en alto.

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

Ahora graficamos la anomalía de Bouguer a lo largo del perfil. Multiplicamos la tercera columna por 111 para convertir entre grados y kilómetros, graficamos 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 - ¿ que hace?.
El perfil de topografía es de una línea sólida roja. Note el escalamiento de los valores de topografía ($4/20) para que podamos 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 vea 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 a 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 la imagen que sigue. El código completo del script viene después de la imagen.

El script completo :

#!/bin/bash

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

#Origen (posicion del grafico en la hoja):
	xshift="2.0"
	yshift="10.0"

proyeccion="M14.0"                #Tipo de proyeccion
paleta="/usr/local/GMT4.5.8/share/cpt/GMT_ocean.cpt"    #Paleta de colores a usar
topografia="/home/carlos/bin/w100s10.Bathymetry.srtm.grd" #Archivo de topografia
ticks="a5f1g5SEWN"    # Describiendo el marco del mapa
illaz="225"     # Direccion del azimut de iluminacion
portrait="-P"	# Si quiero landscape en vez de portrait, dejo esta opcion en blanco
verbose="-V"	# Dejo esta opcion en blanco si no quiero mostrar el "blah blah" de 
                # los comandos en la terminal

#Lineas de costa e informacion de los limites fronterizos. 
#(Para mas detalles, escribe en la terminal: man pscoast)
linea_costa="1.00p/0/0/0"
resolucion="f"	  # Resolucion de las lineas de costa. Puedes elegir:
                # (f)ull, (h)igh, (i)ntermediate, (l)ow, y (c)rude
bordes="-N1/1.00p/0/0/0"  #Limites politicos de los paises

#Tambien se puede añadir:
#rios="-I1/1.00p/0/0/255"

agua="-S0/0/255" #Esto va a graficar de un solo color todo lo que son masas de agua
tierra="-G0/255/0" #Esto va a graficar de un solo color todo lo que es tierra 
lagos="-C50/100/200 -A0/2/4"

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

# Archivo de salida:
psfile="${region}.ps"

# Construyendo el marco del mapa:
psbasemap -B${ticks} -J${proyeccion} -R${limites} -X${xshift} -Y${yshift} ${portrait} \
${verbose} -K > ${psfile}

# Cortando la grilla original a los limites establecidos. 
# grdcut archivo_entrada.grd -Garchivo_salida.grd -Roeste/este/sur/norte [-V] 
grdcut ${topografia} -G${region}.grd -R${limites} -V 

# Iluminando la grilla topografica generada en el grdcut desde un azimut 
# especifico para producir un archivo de intensidad (.int) de iluminacion:
grdgradient ${region}.grd -G${region}.int -A${illaz} -Nt -M

# Graficando la topografia de la region seleccionada con su iluminacion:
grdimage ${region}.grd -C${paleta} -I${region}.int -B${ticks} -J${proyeccion} -R${limites} \
${portrait} ${verbose} -O -K >> ${psfile}
	
###########################################################################################

# Generando la grilla (.grd) de anomalias de Bouguer a partir de los datos de 
# gravedad proporcionados por bouguer.xyz
surface bouguer.xyz -Gbouguer.grd -R-71.5/-64/-29.6/-24.9 -I1m -T0.25 -C0.1 -V

# Creando la paleta de colores asociada a la nueva grilla de gravedad:
grd2cpt bouguer.grd -C/usr/local/GMT4.5.8/share/cpt/GMT_red2green.cpt -T= > bouguer.cpt

# Generando nuevas grillas de iluminacion topografica (topo.int) y de bouguer 
# (bouguer2.grd), de tal manera que se inserten en los mismos puntos de una 
# misma region (-R): 
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

# Graficando la nueva grilla de bouguer, pero con iluminacion asociada a la 
# topografia:
grdimage bouguer2.grd -Cbouguer.cpt -B${ticks} -J${proyeccion} -Itopo.int -R${limites} \
${portrait} ${verbose} -O -K >> ${psfile}

# Graficando los datos de bouguer.xyz como circulos en el mapa para ver donde la 
# grilla es confiable: 
more bouguer.xyz | psxy -J${proyeccion} -R${limites} ${verbose} -Sc0.1c -Cbouguer.cpt \
-W3 -O -K >> ${psfile}

# Insertando contornos a la nueva superficie bouguer2.grd que fue creada:
grdcontour bouguer2.grd -B${ticks} -J${proyeccion} -C50 -A100+k0/0/0+g255/255/255 \
-R${limites} ${portrait} ${verbose} -W0.5p/0/0/0 -Gd5c -O -K >> ${psfile}

# Graficando las costas y limites del mapa:
pscoast -B${ticks} -J${proyeccion} -R${limites} ${portrait} ${verbose} -D${resolucion} \
-W${linea_costa} ${bordes} ${rios}  -O -K >> ${psfile}

# Agregando al grafico los lagos con un color determinado en la variable $lagos:
pscoast -B${ticks} -J${proyeccion} -R${limites} ${portrait} ${verbose} ${lagos} \
-D${resolucion} -W${linea_costa} ${bordes} ${rios} -O -K >> ${psfile}

# Agregando una escala de magnitudes la grilla de bouguer en base a su paleta de 
# colores:
psscale  -D15.7c/2.5c/5c/0.5c -Cbouguer.cpt ${verbose} -Ba100f20:"(mGal)": -O \
-K >> ${psfile}

# Generando un perfil a los a lo largo de la latitud 28S entre las longitudes 
# 73W y 62W, el cual se guarda en un archivo track.xy
project -C-73/-28 -E-62/-28 -G0.01 > track.xy

# Tomando los datos de las grillas de bouguer y de topografia que se encuentran 
# a lo largo del perfil creado y guardandolos en archivo (.xydz):
grdtrack track.xy -Gbouguer2.grd > bouguertrack.xydz
grdtrack track.xy -G${region}.grd > topotrack.xydz

# Creando los bordes de un nuevo mapa, con el objetivo de graficar en el los 
# perfiles creados:
psbasemap -R0/1078/-500/500 -JX14/2 -X0 -Y-4 -Ba200f100 -P -O -K -V >>  ${psfile}

# Graficando la anomalia de bouguer a lo largo del perfil como una curva discontinua:
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}

# Graficando la topografia a lo largo del perfil como una curva roja:
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}

# Agregando una linea recta para ver mejor las variaciones de la anomalia y 
# topografia a lo largo del perfil seleccionado:
psxy << END -R0/1078/-500/500 -JX14/2 -B -P -W1/0/0/0 -O >> ${psfile}
0 0
1078 0
END