513372 - Linux, Scripts y GMT

CLASE 8 - scripts y GMT

back

8.0

En esta clase, trabajaremos con la grilla de topografía de SRTM. Esta tiene una resolución de 30 segundos. ¿Cuántos metros cubre un píxel entonces? .
Se puede bajar ese archivo aca, o está disponible en un pen drive en el laboratorio:

w100s10.Bathymetry.srtm.grd.gz

8.1 Un script para un mapa en GMT

Espero que en la última clase hallan notado tres cosas:

    • GMT es básicamente una serie de comandos en un terminal (en bash).
    • Cuando se hace una modificación pequeña a un mapa, es una molestia correr todos los comandos de nuevo.
    • Cuando esta modificación cambia un parámetro que está en todas los comandos usados, esta modificación tendrá efecto en todos las
       líneas de comando que lleven este parametro.
Por eso existen scripts para generar los mapas. Bajen el script siguiente: ejemplo_mapa1.sh

Este script define los parámetros usados para generar el mapa, y después corre la serie de comandos en GMT para ejecutar el dibujo.
El script esta  abajo, en verde se ve bastante feo, pero si lo abres con un editor de texto como gedit se separa en parámetros, valores, comentarios y comandos.

Note los siguientes parámetros (disculpa, están en inglés. Pueden cambiarlos si quieren):

    • region - define el nombre del región. Note que el gráfico generado (output file) se llama "${region}.ps".
    • bounds - define la latitud y longitud de la región.
    • xshift, yshift - mueven el mapa horizontalmente y verticalmente, respectivamente.
    • projection - en este caso uso Mercator de 14cm.
    • palette, topography - hay que cambiar estos variables para dar el camino hacia sus archivos.
    • ticks - El mapa tiene lados con números cada 2 grados, barras cada 0.5 grados, y líneas cada 1 grado.
    • illaz - azimut de la iluminación

También  hay parámetros para la costa, su resolución, las fronteras graficadas, y los ríos.

Note que mi parámetro para los ríos (#rivers="-I1/1.00p/0/0/255") tiene un # antes, significando que es un comentario nada más. Con eso, no voy a graficar los ríos, pero si quiero graficarlos, simplemente tengo que sacar el "#" para que la línea sea un comando y no un comentario.

#!/bin/bash

region="ejemplo_mapa1"
bounds="-74.5/-71/-38.5/-36"

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

projection="M14.0"                 #projection
palette="/home/matt/GMT/GMT4.5.2/share/cpt/GMT_relief.cpt"   # colour palette
topography="/home/matt/GMT/TOPO_FILES/SRTM30/w100s10.Bathymetry.srtm.grd"    #topography file
ticks="a2f0.5g1SEWN"
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"
    # comment this out if you intend to show bathymetry
# land="-G0/255/0"
        # comment this out if you intend to show topography

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

# 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}
psscale  -D16c/2.5c/5c/0.5c -C${palette} -Ba5000f1000/a5000f1000:"(m)": -O -K >> ${psfile}

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

Si ustedes corren este script, (con los párametros de la paleta y la topografiá cambiados a sus respectivas rutas  en su pc), deben generar esta imagen:



Antes de salir de esta sección, revisa que entienden lo que hace el script y que pueden modificar el script para hacer lo siguiente:

    1) Cambiar el nombre de la imagen generada a bio_bio.ps
    2) Cambiar la región graficada a  una que muestre el altiplano de sudamerica (llamado altiplano.ps, o algo así)
    3) Cambiar la paleta. Prueben globe, ocean, sealand, topo. Noten que algunas de esas son mejores que otros!
    4) Cambiar los ejes: números y barras cada 1 grado, y líneas cada 0.2 grado
    5) Cambiar el azimut de la iluminación para que sea en la dirección opuesta del ejemplo mostrado aquí
    6) Cambiar la resolucion de la costa a full (máximo)
    7) Poner ríos
    8) Sacar del gráfico la escala de la topografía con un gato (#)

8.2 Contornos en GMT

El archivo .grd para la topografía es una grilla, por lo tanto podemos graficar sus contornos usando el comando grdcontour. Por ejemplo, puedo poner

grdcontour ${region}.grd -B${ticks} -J${projection} -C200 -A1000+k0/0/0 -R${bounds} ${portrait} ${verbose} -W0.1p/0/0/0 -Gd10c -O -K >> ${psfile}

dentro del script, entre el grdimage y el psscale.

    • -C200                   significa contornos cada 200m
    • -A1000+k0/0/0     significa anotación cada 1000m, color negro
    • -W0.1p/0/0/0       significa el tamaño, y el color de la lapiz con que dibujamos los contornos
    • -Gd10c               significa 10 cm entre la anotación en cada contorno

El script que uso para generar el siguiente gráfico es ejemplo_mapa2.sh, es bastante similar al primero, solamente cambio la paleta y pongo contornos.



8.3 Personalización del mapa con condicionales

El mapa de arriba esta bastante feo. ¿Sera mejor si elimino la iluminación de la topografía?

Puedo hacer eso borrando la opción -I en el comando grdimage, pero con eso hay que recordar como ponerlo de nuevo. Es mejor  definir unas variables al inicio del script para decir exactamente que queremos, por ejemplo ,
si quiero algo   (1)
no quiero algo (0)         

Donde "algo" puede ser topografía , iluminación , escala , etc.

    topo_option=1
    illumination_option=0
    contour_option=1
    topo_scale_option=0

Despúes en la sección  confección del mapa ( cuerpo del script), solamente grafico de acuerdo a las opciones declaradas al inicio de script , donde el condicional más apropiado para esta selección es  if.
Note además que las opciones que declara al inicio del script ( 1 o 0)  tambien las puede declarar en forma general, usando el condicional case .Este le permitirá declarar las opciones de gráfica como argumentos de entrada , facilitando la tarea de modificar el script poniendo 1 o 0 según lo que desee en
cada ocación.

if [ "${topo_option}" = 1 ]; then
#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

    if [ "${illumination_option}" = 1 ]; then
    # 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}
    else
    grdimage ${region}.grd -C${palette} -B${ticks} -J${projection} -R${bounds} ${portrait} ${verbose} -O -K >> ${psfile}
    fi

    if [ "${contour_option}" = 1 ]; then
    grdcontour ${region}.grd -B${ticks} -J${projection} -C200 -A1000+k0/0/0 -R${bounds} ${portrait} ${verbose} -W0.1p/0/0/0 -Gd10c -O -K >> ${psfile}
    fi

    if [ "${topo_scale_option}" = 1 ]; then
    psscale  -D16c/2.5c/5c/0.5c -C${palette} -Ba5000f1000/a5000f1000:"(m)": -O -K >> ${psfile}
    fi

fi

El script entero se puede encontrar aquí: ejemplo_mapa3.sh. Pruebe este script con varias combinaciones de 1s y 0s en las opciones iniciales para la topografía, y asegurense de entender el gráfico que resulta.

En mi caso (topo_option=1, illumination_option=0, contour_option=1, topo_scale_option=0) Si !  el mapa se ve mejor (soy seco!)



8.4 Concusión

En esta clase vemos que con un script podemos facilmente generar mapas de topografía con opciones de gráfico ( región,  colores usados, contornos, ejes, etc.) que son fácil de cambiar.

En la próxima clase  trabajamos con el mismo script, y veremos los comandos para poner símbolos en los mapas.

back