The LearN@BSD Project

Created - Jan 20, 2008 | Modified - Jan 26, 2008

Homepage -> My GMT -> The GMT Grdimage Primer

The GMT Grdimage Primer

By Przemysław Pawełczyk
Cracow, Poland
Jan 26th, 2008
Version: 0.1

The GMT Gridimage Primer was written for GMT beginners not necessarily from scientific circles. The main purpose of the Primer is to show users how to use the GMT package without walking through GMT documentation which is very extensive as it deals with "the usage subject" in full details. Its excellent but its too fastidious for the starters at least in a hurry.

Figure 0.1. Colca Canyon topomap made using SRTM3 data and the GMT package.

Figure 0.2. Colca Canyon from Space (thanks to Google Map)

1. Why Colca Canyon?

Why did I choose Colca Canyon for the purpose of the Primer? From a simple reasons. I like the place, although I've never been to the Colca River, and because it is the place which was paddled down for the first time in 1982 by a group of Poles, my fellow countrymen.

The Polish Expedition - Websites in Polish
[1] Canion de los Polacos - in Polish (but the pictures!)
[2] Zapasy z gigantem - (Wrestling with a Giant) in Polish (pictures!)
[3] Canoandes '79 Polska Akademicka Eksploracyjna Wyprawa Kajakowa na Rzeki Górskie Argentyny
(Canoandes '79 Polish Academic Kayak Exploration Expedition to Argentinian Mountain Rivers) in Polish

The Polish Expedition - Websites in English
[1] Canoandes Incorporated 1983-2006 - articles, pictures
[2] In the belly of the Earth - direct hyperlink to Colca Kayak Expedition
[3] Paddler Magazine.com - First Descents: Expeditions for the Ages, 10 Top Contemporary River Expeditions

Colca Tourist Information
[1] Colca Canyon
[2] The Colca Valley & Canyon
[3] Colca Travel Pictures
[4] Colca Canyon Picture Gallery
[5] Colca Canyon
[6] Colca River, Peru
[7] Colca Canyon
[8] Where is the Deepest Canyon in the World?

Satellite Images
[1] Google Maps: Colca Canyon from Space
[2] NASA Earth Observatory: Deepest Canyons of the Andes
[3] Earth Science World: Image Bank
[4] Map sources/GeoHack: Colca Canyon
[5] Google Earth

Encyclopedic Information
[1] Wikipedia: Canyon
[2] Wikipedia: Colca Canyon
[3] Wikimedia Commons: Colca Valley

Maps
[1] Google Maps: Colca Canyon
[2] Perry-Castañeda Library Map Collection: Peru
[3] Peru Relief Map
[4] Colca Canyon
[5] Geology.com: Peru
[6] Peru Political Map

2. Scripts used in GMT Grdimage Primer

We will need two scripts. The first one - srtm2gmt.sh which is part of the srtm2gmt-02.tar.gz package and is used for downloading SRTM3 *.hgt.zip files then converting them to grid format, and the second one - grdimage.peru.colca.sh with which we will make topomaps from the gridded tiles.

1. grdimage.peru.colca.sh - universal script used in the GMT Grdimage Primer (shown below). - [7.4 kB]
2. grdimage.primer.www-0.1.tar.gz - The GMT Grdimage Primer with images and grdimage.peru.colca.sh script. - [4 MB]
3. grdimage.primer.www.srtm3-0.1.tar.gz - The GMT Grdimage Primer with images, grdimage.peru.colca.sh script, and Colca SRTM3 original and converted tiles. - [12 MB]
4. srtm3gmt-0.2.tar.gz - see separate description of srtm2gmt package. - [47.5 kB]

#!/bin/sh
#==============================================================================#
# grdimage.peru.colca.sh
#
# Package:   The GMT Grdimage Primer
# Licence:   BSD
# Created:   Jan 26, 2008
# Modified:
# Version:   0.1
#
# Author:    Przemysław Pawełczyk
# Contact:   pp_o2@o2.pl, pp@kv.net.pl, http://pp.kv.net.pl
#
#==============================================================================#

# ========== set your GMT variables ============================================
# here title1 placement
gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1.5c
#gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1c

# ========== enter a region griddinates and its title ==========================

# --- Eastern SRTM tile ------------------------------------
#latlon="S16W072"
#range="-72/-71/-16/-15"
#legend="f288.2/-15.9/0/20"
#title1="Colca Canyon East"
#title2="288.5 -14.9 20 0 0 CM"

# --- Western SRTM tile ------------------------------------
latlon="S16W073"
range="-73/-72/-16/-15"
legend="f287.2/-15.9/0/20"
title1="Colca Canyon West"
title2="287.5 -14.9 20 0 0 CM"

# --- Both tiles grid together -----------------------------
#latlon="S16W07372"
#range="-73/-71/-16/-15"
#legend="f287.2/-15.9/0/20"
#title1="Colca Canyon"
#title2="288 -14.9 20 0 0 CM"

# ========== define grid file ==================================================
grid="${latlon}.grd"

# ========== define grid file with illumination ================================
grid_i="${latlon}_i.grd"

# ========== define scale ======================================================
scale=10

# ========== define dots per inch ==============================================
dpi=300

# ========== define plot width =================================================
width=14

# ========== define projection =================================================
proj="M$width"

# ========== define color table position/shape =================================
# --- position X -------------------------------------------
d1=$width
d1=$(($d1+1))
d1=$d1.5
# --- middle point for position Y --------------------------
d2=5
# --- width or height --------------------------------------
d3=10
# --- bar width --------------------------------------------
d4=1 

# ========== enter character coding ============================================
# change the entry of CHAR_ENCODING in GMTDEFAULTS:
# http://gmt.soest.hawaii.edu/gmt/doc/html/gmtdefaults.html#CHAR_ENCODING
# e.g. when ISO-8859-2 then iso="lat2"
#      when ISOLatin1+ then iso="lat1p"
# the naming is up to a user, it serves outup.ps names only
iso="lat1p"

# ---------- enter ps corename ----------------------------
ps_name="peru.colca"

# ---------- enter CPT corename ----------------------------
cpt_name="topo"

# ========== make new cpt table or use existing one ============================
# if you set to ZERO make your selections of existing color tables manually
# (in ELSE part)
# make_cpt=0/1 where 0=no, 1=yes
#
make_cpt=1
if [ "$make_cpt" -eq 1 ]; then

  # ========= select method = 1,2 ==========================
  # cpt_method=1/2
  #
  cpt_method=2

  if [ "$cpt_method" -eq 1 ]; then
  # ---------- making CPT file, 1st method -----------------
    echo "-----> Creating $cpt2... method 1"
    cpt1="GMT_${cpt_name}.cpt"
    cpt2="pp_${cpt_name}_m1.cpt"
    grd2cpt $grid -C$cpt1 -E$scale > $cpt2
  else
  # ---------- making CPT file, 2nd method -----------------
    echo "-----> Creating $cpt2... method 2"
    cpt2="pp_${cpt_name}_m2.cpt"
    makecpt -C${cpt_name} -T0/6500/500 > $cpt2
  fi
else
#  cpt2="GMT_${cpt_name}.cpt"
#  cpt2="pp_${cpt_name}_m1_red.cpt"
  cpt2="pp_${cpt_name}_m2_red.cpt"
fi

# ========== show grdimage =====================================================
# grdim=gi0/gi1 where 0=no, 1=yes
grdim="gi0"
#grdim="gi1"

# ========== show illumination =================================================
# ilum=i0/i1 where 0=no, 1=yes
ilum="i0"
#ilum="i1"

# ========== show contours =====================================================
# contours=c0/c1 where 0=no, 1=yes
contours="c0"
#contours="c1"

# ========== show rivers =======================================================
# rivers=r0/r1 where 0=no, 1=yes
rivers="r0"
#rivers="r1"

# ========== show map border, Unix stamp, legend ===============================
# basemap=b0/b1 where 0=no, 1=yes
#basemap="b0"
basemap="b1"

# ========== show color table ==================================================
# ctable=ct0/ct1 where 0=no, 1=yes
ctable="ct0"
#ctable="ct1"

# ========== show second title =================================================
# subtitle=2t0/2t1 where 0=no, 1=yes
subtitle="2t0"
#subtitle="2t1"

# ========== ps output name ====================================================
psout1="$ps_name.${latlon}.${dpi}dpi.$cpt2.s$scale.$ilum."
psout2="$subtitle.$contours.$rivers.$iso"
psout="$psout1$psout2".ps

################################################################################
echo "=====> Creating $psout..."

# In case you chose to draw only small subset off GMT commands here.
# Writes nothing, intended to open the $psout file only!
# Pay attention to -K option which opens postscript plot.
echo "-----> Openning ps file..."
echo "$title2  " | pstext -R$range -J$proj                       -K     > $psout


if [ "$make_cpt" -eq 1 ]; then
  echo "-----> Using new $cpt2..."
else
  echo "-----> Using existing $cpt2..."
fi

if [ ! -e "$grid_i" ]; then
  echo "-----> Creating gradient..."
  grdgradient $grid -A0 -Nt -G$grid_i -V
else
  echo "-----> Gradient file exists..."
fi

if [ "$grdim" == "gi1" ]; then
  if [ "$ilum" == "i1" ]; then
    echo "-----> Creating grdimage with illumiation..."
    grdimage $grid -R$range -J$proj -I$grid_i -E$dpi -C$cpt2  -V -K     > $psout
  else
    echo "-----> Creating grdimage w/o illumination..."
    grdimage $grid -R$range -J$proj           -E$dpi -C$cpt2  -V -K     > $psout
  fi
fi

if [ "$contours" == "c1" ]; then
  echo "-----> Creating grdcontours..."
  grdcontour $grid -R$range -J$proj -C$cpt2 -A- -Wc0.05          -K -O >> $psout
fi

if [ "$rivers" == "r1" ]; then
  echo "-----> Creating rivers..."
  pscoast -R$range -J$proj -Df -Ir/2p,blue -Sblue                -K -O >> $psout
fi

if [ "$basemap" == "b1" ]; then
  echo "-----> Creating Legend and Unix stamp..."
  #psbasemap -R$range -J$proj -L$legend -Bf0.1a0.2g0.2:."$title1":WeSn \
  #       -U"Przemys³aw Pawe³czyk - GMT 4.2.1 - The GMT Grdimage Primer" \
  #                                                              -K -O >> $psout
  psbasemap -R$range -J$proj -L$legend -Bf0.1a0.2g0.2:."$title1":WeSn \
            -U"Przemyslaw Pawelczyk - GMT 4.2.1 - The GMT Grdimage Primer" \
                                                                 -K -O >> $psout
fi

if [ "$ctable" == "ct1" ]; then
  echo "-----> Creating color scale..."
  psscale -D$d1/$d2/$d3/$d4 -C$cpt2 -B/:m:                       -K -O >> $psout
fi

if [ "$subtitle" == "2t1" ]; then
  echo "-----> Writing second title..."
  echo "$title2 $cpt2" | pstext -R$range -J$proj -N                 -O >> $psout
fi

# In case you chose not to draw second title which
# has -O option closing $psout file.
# Writes nothing, intended to close the $psout file only! 
echo "-----> Closing ps file..."
echo "$title2  " | pstext -R$range -J$proj                          -O >> $psout

################################################################################
exit

As you can see the script was written to be the most universal tool as possible. You can turn on/off every essential GMT command without turning to (un)commenting script lines. We are not using 0/1 (zero/one) but "string" switches because they are also used to name output postscript file. Hence for example the second title (subtitle) uses "2t0/2t1" switches (see the script few lines above).

3. Gridding SRTM3 Tiles for Colca Area

GMT commands works mainly on grid data files. As GMT package is enormously versatile in dealing with geodata they can also use data files in text form. But as the grdimage command deals with gridded geodata (its name suggest so) we will have to obtain colca.grd file in the first place. We'll make use of srtm2gmt-0.2.tar.gz package to do so. Lets make the following directories:

   ${HOME}/data/gmt/gmt.scripts
   ${HOME}/data/gmt/gmt.data

Now put srtm2gmt.sh script and srtm3tiles.ftp file in the first one, both comes from srtm2gmt-02.tar.gz package. Change working directory for gmt.scripts and using any editor open srtm2gmt.sh. Set some variables accordingly and run the script (what to change and how to do it is described on srtm2gmt package documentation page). For example, set at least the four items:

Colca area griddinates.
Directory of your GMT package.
Directory where the Colca tiles will be placed.
Endian type for conversion.

ATTENTION:
1) Colca River and Colca Canyon run through area covered by 2 SRTM3 tiles at least. We must set the script to download both.
2) Put griddinates without leading ZEROs. Example: use "W72" against "W072" (don't worry the script will detect the wrong entries).

Here is a part of the srtm2gmt.sh script where we have to set the new values:

   #----------------------------------------------------------#
   # 1. Defining variables and constants
   #----------------------------------------------------------#
   ############################################################
   # Give your griddinates                                    #
   ############################################################
   lat2=S16
   lat1=S16

   lon2=W72
   lon1=W73

   #--------------------------------------#
   # Region and target area
   #--------------------------------------#
   region="peru"
   target="colca"

   descr="Getting Colca Canyon tiles."

   #--------------------------------------#
   # Intel (1) or Motorola (0)?
   #--------------------------------------#
   this_comp=1

   if [ "$this_comp" -eq 1 ]
   then
   # HGT files comes in Big-endian format and they need to be sometimes converted.
   # The author's GMT 4.2.1 was compiled under Debian (he didin't use DEB package)
   # and the endian had to be set explicit.
     endian=" B"
   else
     endian=""
   fi

   #--------------------------------------#
   # GRD database directory 
   #--------------------------------------#
   gmt_grids_dir="${HOME}/data/gmt/gmt.data/srtm3.${region}.${target}/"

   if [ ! -d $gmt_grids_dir ]
   then
     mkdir $gmt_grids_dir
   fi

   ....... [cut, See the script for details] .....
   #--------------------------------------#
   # GMT grdraster.info directory 
   #--------------------------------------#
   gmt_dbase_dir="/opt/GMT4.2.1/share/dbase/"

   ..... [cut] ......
   #------------------------------------------------------------------------------#
   # BLOCK_1 - end
   #------------------------------------------------------------------------------#

When everything goes well the new directory of ~/data/gmt/gmt.data/srtm3.peru.colca will be created, Colca tiles downloaded and then converted from HGT format to GRD format used by GMT package. The directory will contain new files shown on Fig. 3.1. displaying Colca directory contents.

Figure 3.1. Colca directory contents.

Next, we will have to copy/move the grdimage.peru.colca.sh script to the ~/data/gmt/gmt.data/srtm3.peru.colca directory. With that all arrangements to make topomaps of Colca Canyon will be fulfilled.

Figure 3.2. Colca directory contents after running grdimage.peru.colca.sh. (GMT_topo.cpt file was copied manually prior to running the script, but if GMT is installed correctly it is not needed.)

A few words on new files

.gmtcommands4

A. Documentation pages:
[1] 4.2.1 Overview and the .gmtdefaults4 file
[2] 4.2.2 Changing GMT Defaults

   # GMT common arguments shelf
   -B/:m:
   -JM14
   -R-73/-72/-16/-15
   -y1c
   -jM14
   EOF
   R-73/-72/-16/-15
   -y1c
   -jM14
   EOF

In most cases a single topomap is created by several GMT commands. The most convenient way to work with them is to put them together in a shell script. The GMT commands are "very Unix like" what means that every one of them use a lot of options e.g. -R (map range), -J (projection or xy(z) plot), -B (map's frame), etc. Have a look at our script below (some options and shell commands were removed for lucidity):

   (...)
   grdimage   $grid.grd -R$range -J$proj        -C$cpt2        -K  > $psout
   grdcontour $grid.grd -R$range -J$proj        -C$cpt2     -O -K >> $psout
   pscoast              -R$range -J$proj                    -O -K >> $psout
   psbasemap            -R$range -J$proj -Bf0.1a0.2g0.2WeSn -O -K >> $psout
   (...)

Some of the options must be repeated, here $range (-R) and projection (-J). Writing them several times, even using shell variables like $range and $proj instead of e.g. "-R-73/-72/-16/-15" and "M14c" respectively, is too much for some. The GMT package authors took the subject into consideration and invented a way around. Simply speaking the repetitive data were stored in .gmtcommands4 hidden file (with dot in its name at the beginning). Now, it suffices to define the options once (in the first invocation) and then write only option letter. The missing data are retrieved from the hidden file. Very clever.

NOTE:
In the above example options -R and -J are repeated with corresponding data ($...) due to versatility of the script. Both -R and -J must be declared once. But as the script can turn off every GMT command it could have happen that the data were not defined. Using $range and $proj variables for every GMT command we gain certitude that the script will run without errors. But of course the GMT command could be written fully in GMT way:

   (...)
   grdimage   $grid.grd -R$range -J$proj   -C$cpt2        -K  > $psout
   grdcontour $grid.grd -R       -J        -C$cpt2     -O -K >> $psout
   pscoast              -R       -J                    -O -K >> $psout
   psbasemap            -R       -J -Bf0.1a0.2g0.2WeSn -O -K >> $psout
   (...)

If we chose other variable contents e.g. adding to the variable contents also "-R" and/or "-J" GMT option tags, we would be obliged to repeat them for every GMT command demanding them. For example:

   range="-R-73/-72/-16/-15"
   proj="-JM14c"

BTW, "M" means Merkator projection, "14c" - 14 cm map frame width.

   (...)
   grdimage   $grid.grd $range $proj        -C$cpt2        -K  > $psout
   grdcontour $grid.grd $range $proj        -C$cpt2     -O -K >> $psout
   pscoast              $range $proj                    -O -K >> $psout
   psbasemap            $range $proj -Bf0.1a0.2g0.2WeSn -O -K >> $psout
   (...)

.gmtdefaults4

A. Documentation pages:
     [1] GMT man
     [2] GMTDEFAULTS man
     [3] 4.2.1 Overview and the .gmtdefaults4 file
     [4] 4.2.2 Changing GMT Defaults

The .gmtdefaults4 file is a copy of .gmtdefaults_SI file found in ${GMT}/share/conf/ directory. This file is generated in current directory when the main script contains definitions of new GMT variables. Because we defined in our script a new title placement - we moved it a bit higher on the plot to make place for the second title (subtitle) - so the .gmtdefaults4 was created. Look into the script what we have changed:

   # ========== set your GMT variables ================= ==========================
   # here title1 placement
   gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1.5c
   #gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1c

pp_topo_m2.cpt

This is our own color table created by GMT where pp defines a user, topo is the corename of default GMT_topo.cpt file, and the m2 describes a second method used for creation of the new color table. The first method makes use of GMT grd2cpt program which creates new color table taking into account 1) elevation data from Colca grid file and 2) existing GMT default cpt file (here GMT_topo.cpt). The second method utilizing GMT makecpt program makes color table taking data from 1) a user defined color scale (-T option) and 2) default GMT cpt file (also GMT_topo.cpt). You don't have to use full name with option -C, it suffices to write cpt corename (without prefix "GMT_" and suffix ".cpt"). Here's an excerpt from the Primer's main script:

   # ========== make new cpt table or use existing one ============================
   # if you set to ZERO make your selections of existing color tables manually
   # (in ELSE part)
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=1
   if [ "$make_cpt" -eq 1 ]; then

     # ========= select method = 1,2 ==========================
     # cpt_method=1/2
     #
     cpt_method=2

     if [ "$cpt_method" -eq 1 ]; then
     # ---------- making CPT file, 1st method -----------------
       echo "-----> Creating $cpt2... method 1"
       cpt1="GMT_${cpt_name}.cpt"
       cpt2="pp_${cpt_name}_m1.cpt"
       grd2cpt $grid -C$cpt1 -E$scale > $cpt2
     else
     # ---------- making CPT file, 2nd method -----------------
       echo "-----> Creating $cpt2... method 2"
       cpt2="pp_${cpt_name}_m2.cpt"
       makecpt -C${cpt_name} -T0/6500/500 > $cpt2
     fi
   else
   #  cpt2="GMT_${cpt_name}.cpt"
   #  cpt2="pp_${cpt_name}_m1_red.cpt"
     cpt2="pp_${cpt_name}_m2_red.cpt"
   fi

Below is the pp_topo_m2.cpt contents (the first and the fifth columns sets height range, the remaining triple numbers define colors in RGB format). We will return to the subject later.

   # cpt file created by: makecpt -Ctopo -T0/6500/500
   #COLOR_MODEL = RGB
   #
      0	173	138	230	500	173	138	230
    500	138	176	230	1000	138	176	230
   1000	140	243	217	1500	140	243	217
   1500	154	243	137	2000	154	243	137
   2000	234	230	133	2500	234	230	133
   2500	232	168	136	3000	232	168	136
   3000	116	163	179	3500	116	163	179
   3500	189	206	118	4000	189	206	118
   4000	238	224	175	4500	238	224	175
   4500	252	238	219	5000	252	238	219
   5000	255	250	246	5500	255	250	246
   5500	255	252	249	6000	255	252	249
   6000	255	254	253	6500	255	254	253
   B	236	140	255
   F	255	255	255
   N	128	128	128

Thus set grdimage.peru.colca.sh yields the following plot (run it after setting according to the text above):

Figure 3.3. The first topomap of Colca Canyon made with the help of GMT Grdimage Primer's grdimage.peru.colca.sh script.

IMPORTANT NOTE:
The resulting postscript filename:
peru.colca.S16W073.300dpi.pp_topo_m2.cpt.s10.i1.2t1.c0.r0.lat1p.ps
is very long but it was designed intentionally. Every name delimited with dots shows one of the script's option. Sometimes it is a number e.g. s10 (scale 10), sometimes string option e.g. c0 (no contours). Compare the dot delimited strings to the Primer script's variables.

4. The GMT Commands

4.0. The GMT Basics

A. Documentation pages:
     [1] GMT man
     [2] GMTDEFAULTS man
     [3] GMTSET man
     [4] GRDINFO man
     [5] 4.4.6 Plot Overlays: The -K -O options

All GMT output postscript files are created with the help of operating system's shell scripts and the GMT commands (programs). Every shell script comprises of several blocks containing one or a few GMT commands and auxiliary shell/GMT commands.

The output postscript files are build layer after layer. For example, the first layer may contain picture made by grdimage command, the second layer can contain lat/lon grid lines accompanied by Unix timestamp for example, the third layer can create color scale bar, and the fourth layer can add a finishing text.

All the layers are imposed one onto another with -O and -K options (see GMT documentation [5]). The process goes like this:

   layer 1 -K    >  output.ps
   layer 2 -K -O >> output.ps
   layer 3 -K -O >> output.ps
   layer 4    -O >> output.ps

This chapter presents several GMT programs. We will shortly described them and will show their results. In other words we will describe step by step what to do to get partial output plots and full blown topographic visuallizations. We will outline the usage of the following programs (in order of appearance): psbasemap, pstext, makecpt, grd2cpt, psscale, grdimage, grdgradient, and grdcontour.

And please, do not forget to look into GMT documentation pertaining to reading topics of the Primer.

4.1. psbasemap

A. Documentation pages:
     [1] 1.4.1 Linear projection
     [2] 4.2.1 Overview and the .gmtdefaults4 file
     [3] 4.4.3 Map frame and axes annotations: The -B option
     [4] 4.4.3.1 Geographic basemaps
     [5] 4.4.3.2 Cartesian linear axes
     [6] PSBASEMAP man
     [7] Timestamp

Let's change the grdimage.peru.colca.sh script in a way to achieve only a map's frame, namely leaving valid only psbasemap GMT command. For example:

B. Script modifications:

Dominish HEADER offset to 1 cm:

   # ========== set your GMT variables ============================================
   # here title1 placement
   #gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1.5c
   gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1c

Exclude all options setting corresponding variables to "x0" (where "x" means a character, see the script contents) but do not exclude psbasemap:

   # ========== show map border, Unix stamp, legend ===============================
   # basemap=b0/b1 where 0=no, 1=yes
   #basemap="b0"
   basemap="b1"

Get rid of grdimage too.

   # ========== show grdimage =====================================================
   # grdim=gi0/gi1 where 0=no, 1=yes
   grdim="gi0"
   #grdim="gi1"

Run the script

We are interested with the options standing near to psbasemap command. They need to be explained.

   if [ "$basemap" == "b1" ]; then
     echo "-----> Creating Legend and Unix stamp..."
     #psbasemap -R$range -J$proj -L$legend -Bf0.1a0.2g0.2:."$title1":WeSn \
     #       -U"Przemys³aw Pawe³czyk - GMT 4.2.1 - The GMT Grdimage Primer" \
     #                                                              -K -O >> $psout
     psbasemap -R$range -J$proj -L$legend -Bf0.1a0.2g0.2:."$title1":WeSn \
               -U"Przemyslaw Pawelczyk - GMT 4.2.1 - The GMT Grdimage Primer" \
                                                                    -K -O >> $psout
   fi

Colca Canyon embraces two tiles. Both covers an area defined in geographical griddinates. They are declared in $range variable. The options of our interest are "-B" and "-U" mainly. But we will describe the "-L" option as well. And do not forget to look into documentation listed above as well! The "-B" string (in our example) comprises of three parts:

Option -B

f0.1a0.2g0.2 (see documentation [3],[4],[5]): "f" means frame tick spacing, marking black and white frame parts every 0.1 degree span, "a" means annotation tick spacing (one every frame, if there were more annotation texts they would overlapped on themselves), "g" stands for grid lines spacing (one every frame, "covering" annotation span. (Set f0.1a0.1g0.1 as an exercise and run our script once more.)
:."$title1": - this string defines plot main title. The dot denotes the title, see documentation [3].
WeSn - from [4]: "By default, all 4 boundaries are plotted (referred to as W, E, S, N). To change the default, append the code for only those axes you want (e.g., WS for standard lower-left x- and y-axis system). Upper case (e.g., W) means draw axis/tickmarks AND annotate it, whereas lower case (e.g., w) will only draw axis/tickmarks."

Option -U

-U"Przemyslaw Pawelczyk - GMT 4.2.1 - The GMT Grdimage Primer" - it makes GMT Unix system timestamp with the attached string at the left bottom part of output plot.

Option -L

f287.2/-15.9/0/20 - (legend="f287.2/-15.9/0/20") creates fancy map scale showing distance in km. "f" - means fancy scale, "287.2" - lon of the center point of the scale, "-15.9" - lat of the center point of the scale, "0" - lat where the scale is calculated, "20" - show length in km (default)

Figure 4.1. Basic map frame with timestamp at the bottom.

4.2. pstext

A. Documentation pages:
     [1] 2.2 Plotting text strings
     [2] PSTEXT MAN
     [3] 4.16 Character escape sequences
     [4] pstext Example: 7.17 Images clipped by coastlines
     [5] Example: 7.7 A simple location map

Let's add subtitle to our ps plot now. We have to make two changes to our script.

B. Script modifications:

Bring back HEADER offset to 1.5 cm:

   # ========== set your GMT variables ============================================
   # here title1 placement
   gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1.5c
   #gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1c

And include subtitle (second title):

   # ========== show second title =================================================
   # subtitle=2t0/2t1 where 0=no, 1=yes
   #subtitle="2t0"
   subtitle="2t1"

Run the script

Figure 4.2. Basic map frame with timestamp and a subtitle.

4.3. makecpt and grd2cpt

A. Documentation pages:
     [1] 4.1 Cpt files
     [2] M. Built-in color palette tables
     [3] Palettes from GMT
     [4] 4.15 Color palette tables
     [5] I. Color Space -- The final frontier
     [6] MAKECPT MAN
     [7] GRD2CPT MAN
     [8] CPT City

Now we try to add color table using psscale command. But to do so we must create the table. There are two GMT commands to do that, I wrote about that in the chapter: 3. Gridding SRTM3 Tiles for Colca Area, in part entitled "A few words on new files", in section "pp_topo_m2.cpt". So let's have a quick look on them.

   # ========== make new cpt table or use existing one ============================
   # if you set to ZERO make your selections of existing color tables manually
   # (in ELSE part)
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=1
   if [ "$make_cpt" -eq 1 ]; then

     # ========= select method = 1,2 ==========================
     # cpt_method=1/2
     #
     cpt_method=2

     if [ "$cpt_method" -eq 1 ]; then
     # ---------- making CPT file, 1st method -----------------
       echo "-----> Creating $cpt2... method 1"
       cpt1="GMT_${cpt_name}.cpt"
       cpt2="pp_${cpt_name}_m1.cpt"
       grd2cpt $grid -C$cpt1 -E$scale > $cpt2
     else
     # ---------- making CPT file, 2nd method -----------------
       echo "-----> Creating $cpt2... method 2"
       cpt2="pp_${cpt_name}_m2.cpt"
       makecpt -C${cpt_name} -T0/6500/500 > $cpt2
     fi
   else
   #  cpt2="GMT_${cpt_name}.cpt"
   #  cpt2="pp_${cpt_name}_m1_red.cpt"
     cpt2="pp_${cpt_name}_m2_red.cpt"
   fi

grd2cpt

We do not want reinvent a wheel so let's read a short excerpt from [7]: "grd2cpt reads a grid file and writes a color palette (cpt) file to standard output. The cpt file is based on an existing master cpt file of your choice,...". What we need to explain are its options used in our script.

$grid - grid file of plotted area
-C$cpt1 - default GMT cpt file used as a color template for new color table
-E$scale - number of colors related to number of elevation spans (max height divided by $scale).

makecpt

Another short excerpt from official GMT documentation [6]: "makecpt is a utility that will help you make color palette tables (cpt files). You define an equidistant set of contour intervals or pass your own z-table, and create a new cpt file based on an existing master cpt file."

-C${cpt_name} - default GMT cpt file used as a color template for new color table
-T0/6500/500 - three data defining new color table: height range in meters from minimum (here 0), to maximum (6500), and height span (500). When relative height (6500m-0m) is divided by the height span it will define a number of steps in new color table.

How they look like on the plots we will be able to see after using GMT psscale command.

4.4. psscale

A. Documentation pages:
     [1] Example: 7.2 Image presentations
     [2] Example: 7.16 Gridding of data, continued
     [3] Example: 7.18 Volumes and Spatial Selections
     [4] PSSCALE MAN

Time to present psscale. It requires few options:

   if [ "$ctable" == "ct1" ]; then
     echo "-----> Creating color scale..."
     psscale -D$d1/$d2/$d3/$d4 -C$cpt2 -B/:m: -O -K >> $psout
   fi

-D$d1/$d2/$d3/$d4 - we plot vertical color bar: d1 - defines color bar starting point in X direction (left-right), d2 - defines middle point of color bar in Y direction (up-down), d3 - defines color bar height, d4 - points to color bar width.
-C$cpt2 - nothing to explain, color table to plot.
-B/:m: - color table label. You can use -B:title:/:label: option as well.

Time to introduce changes to our script to obtain plots containing color scale bar. Pay attention to the scales in the two color bars.

B. Script modifications:

Let's start to use new cpt file:

   # ========== make new cpt table or use existing one ============================
   # if you set to ZERO make your selections of existing color tables manually
   # (in ELSE part)
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=1

Let's select first method with grd2ecpt command:

     # ========= select method = 1,2 ==========================
     # cpt_method=1/2
     #
     cpt_method=1

Let's enable script option to plot the color table:

   # ========== show color table ==================================================
   # ctable=ct0/ct1 where 0=no, 1=yes
   #ctable="ct0"
   ctable="ct1"

Run the script

Figure 4.4.1. grd2cpt(GMT_topo.cpt) -> pp_topo_m1.cpt

C. Script modifications:

Let's select now the second method with makecpt command:

     # ========= select method = 1,2 ==========================
     # cpt_method=1/2
     #
     cpt_method=2

Run the script

Figure 4.4.2. makecpt(GMT_topo.cpt) -> pp_topo_m2.cpt

Changing manually created cpt files.

If you dont like the new color tables, for example, they have too many close to white colors for higher elevations, you can make your own modifications. Let's do it to pp_topo_m2.cpt file. Listing below presents its contents:

   #	cpt file created by: makecpt -Ctopo -T0/6500/500
   #COLOR_MODEL = RGB
   #
      0	173	138	230	500	173	138	230
    500	138	176	230	1000	138	176	230
   1000	140	243	217	1500	140	243	217
   (...)
   5000	255	250	246	5500	255	250	246
   5500	255	252	249	6000	255	252	249
   6000	255	254	253	6500	255	254	253
   B	236	140	255
   F	255	255	255
   N	128	128	128

Copy the file to pp_topo_m2_red.cpt and make changes to the heights of 5500-6000 and 6000-6500m.

   (...)
   5500  255 120 120  6000  255 120 120
   6000  255   0   0  6500  255   0   0
   (...)

D. Script modifications:

We must also change a little our script - "make_cpt" set to zero, and uncomment (unhash) cpt2="pp_${cpt_name}_m2_red.cpt":

   # ========== make new cpt table or use existing one ============================
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=0
   if [ "$make_cpt" -eq 1 ]; then
   (......)
   else
   #  cpt2="GMT_${cpt_name}.cpt"
   #  cpt2="pp_${cpt_name}_m1_red.cpt"
     cpt2="pp_${cpt_name}_m2_red.cpt"
   fi

Run the script

Figure 4.4.3. Color bar based on pp_topo_red.cpt file.

4.5. grdimage

A. Documentation pages:
     [1] 4.3 Color images
     [2] Cook-book, Examples 02: 7.2 Image presentations
     [3] Cook-book, Examples 17: 7.17 Images clipped by coastlines
     [4] GRDIMAGE man
     [5] B.2 Grid files
     [6] 4.17 Grid file format specifications

We come to the main subject of the Primer, to grdimage GMT command/program. It requires few options:

if [ "$grdim" == "gi1" ]; then
  if [ "$ilum" == "i1" ]; then
    echo "-----> Creating grdimage with illumiation..."
    grdimage $grid -R$range -J$proj -I$grid_i -E$dpi -C$cpt2  -V -K     > $psout
  else
    echo "-----> Creating grdimage w/o illumination..."
    grdimage $grid -R$range -J$proj           -E$dpi -C$cpt2  -V -K     > $psout
  fi
fi

$grid - grid file of plotted area
-R$range - geographical range
-J$proj - map projection and width of resulting map
-E$dpi - dots per inch definition of postscript plot
-C$cpt2 - CPT color table used to color elevation spans

B. Script modifications:

We use grdimage to see the difference between the two color scales desrcibed in former section (4.4). As we have all variables needed by grdimage settled, the only thing left to do is enabling grdimage in our script:
```
   # ========== show grdimage =====================================================
   # grdim=gi0/gi1 where 0=no, 1=yes
   #grdim="gi0"
   grdim="gi1"
```

Check also if "make_cpt" is set to zero, and uncomment (unhash) cpt2="pp_${cpt_name}_m2_red.cpt":

   # ========== make new cpt table or use existing one ============================
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=0
   if [ "$make_cpt" -eq 1 ]; then
   (......)
   else
   #  cpt2="GMT_${cpt_name}.cpt"
   #  cpt2="pp_${cpt_name}_m1_red.cpt"
     cpt2="pp_${cpt_name}_m2_red.cpt"
   fi

Run the script

Figure 4.5.1. Colca Canyon topomap based on pp_topo_m2_red.cpt file.

C. Script modifications:

After running our script with color table set to pp_topo_m2_red.cpt file, change the make_cpt variable to 1 (see the listing below) and run the script once more.

   # ========== make new cpt table or use existing one ============================
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=1
   if [ "$make_cpt" -eq 1 ]; then
   (......)
   else
   #  cpt2="GMT_${cpt_name}.cpt"
   #  cpt2="pp_${cpt_name}_m1_red.cpt"
     cpt2="pp_${cpt_name}_m2_red.cpt"
   fi

Run the script

Figure 4.5.2. Colca Canyon topomap based on pp_topo_m2.cpt file.

4.6. grdgradient

A. Documentation pages:
     [1] 4.2 Illumination and intensities
     [2] 4.3 Color images - (-Ne and -Nt options explained)
     [3] Cook-book, Example 02: 7.2 Image presentations - (grdgradient and grdimage)
     [4] H.4 Hints
     [5] GRDGRADIENT man

All previous plots seemed "flat" devoid of topographic details. GMT is able to "up" the flatness into third dimension. What is it all about and how to do it is well explained in [1] and [5]. We will simplify the issue and confine ourselves to describing the subject in a few sentences. The 3D data, called also shade/illumination data, are created by grdgradient program. Grdgradient allows to define light source's (e.g. Sun's) shading direction and its elevation plus other options like gradient computational methods. We will use the most default ones. Pay attention to the shell commands, once the illumination file was created there's no point in its recreation every time our script is running.

   if [ ! -e "$grid_i" ]; then
     echo "-----> Creating gradient..."
     grdgradient $grid -A0 -Nt -G$grid_i -V
   else
     echo "-----> Gradient file exists..."
   fi

$grid - grid file of plotted area
-G$grid_i - new grid file with illumination data
-A0 - 0 degree azimuthal direction what means a light source is set to the North. It happens that 180 degree value creates upside-down illusion of plotted surface - what was down (e.g. river valleys) are then plateaus. Will show the "brain deficiency" effect.
-Nt - see [2], some texts say to set the "-N" option to t0.75; [2] suggests -Ne0.8, [5] points to -Ne0.6 value.
-V - say something about the process

Let's check that our script have set the following options:

B. Script modifications:

   # ========== set your GMT variables ============================================
   #
   # here title1 placement
   #gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1c
   gmtset HEADER_FONT_SIZE 28 HEADER_OFFSET 1.5c

   # ========== make new cpt table or use existing one ============================
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=1
   if [ "$make_cpt" -eq 1 ]; then
     # ========= select method = 1,2 ==========================
     # cpt_method=1/2
     #
     cpt_method=2
     (...)
   else
     (...)
   fi

   # ========== show grdimage =====================================================
   # grdim=gi0/gi1 where 0=no, 1=yes
   #grdim="gi0"
   grdim="gi1"
   
   # ========== show illumination =================================================
   # ilum=i0/i1 where 0=no, 1=yes
   #ilum="i0"
   ilum="i1"

   # ========== show contours =====================================================
   # contours=c0/c1 where 0=no, 1=yes
   contours="c0"
   #contours="c1"

   # ========== show rivers =======================================================
   # rivers=r0/r1 where 0=no, 1=yes
   rivers="r0"
   #rivers="r1"

   # ========== show map border, Unix stamp, legend ===============================
   # basemap=b0/b1 where 0=no, 1=yes
   #basemap="b0"
   basemap="b1"

   # ========== show color table ==================================================
   # ctable=ct0/ct1 where 0=no, 1=yes
   #ctable="ct0"
   ctable="ct1"

   # ========== show second title =================================================
   # subtitle=2t0/2t1 where 0=no, 1=yes
   #subtitle="2t0"
   subtitle="2t1"

Let's try new "-N" option value. Please set it to -Ne0.6.

   if [ ! -e "$grid_i" ]; then
     echo "-----> Creating gradient..."
     grdgradient $grid -A0 -Ne0.6 -G$grid_i -V
   else
     echo "-----> Gradient file exists..."
   fi

Run the script

Figure 4.6.1. Colca Canyon topomap based on pp_topo_m2.cpt file and with shading/illumination supplement.

C. Script modifications:

Please, change one option now. Set make_cpt to "0", and uncomment our modified CPT file: cpt2="pp_${cpt_name}_m2_red.cpt".

   # ========== make new cpt table or use existing one ============================
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=0
   if [ "$make_cpt" -eq 1 ]; then
     (...)
   else
   #  cpt2="GMT_${cpt_name}.cpt"
   #  cpt2="pp_${cpt_name}_m1_red.cpt"
     cpt2="pp_${cpt_name}_m2_red.cpt"
   fi

It concerns pp_topo_m2_red.cpt file. Writing the Primer I decided to change our modified color table (pp_topo_m2_red.cpt). The last two red ranges between 5500-6000-6500 meters were poorly descerned. I changed the last but one range's color from 150 to 200, making its red hue much paler.
```
   (...)
   5000	255	250	246	5500	255	250	246
   5500	255	200	200	6000	255	200	200
   6000	255	0	0	6500	255	0	0
   (...)
```
Run the script

Figure 4.6.2. Colca Canyon topomap based on pp_topo_m2_red.cpt file and with shading/illumination supplement.

D. Script modifications:

The normalization option is very intriguing. Let's cite documentation [2] (at the page bottom): "Both -Ne and -Nt yield well behaved gradients. Personally, we prefer to use the -Ne option; the value of norm is subjective and you may experiment somewhat in the 0.5-5 range." As we used -Ne0.6 what was close to the lower limit let's now select upper limit value of -Ne5.0. But do not forget to remove $grid_i file before running the script!
```
   if [ ! -e "$grid_i" ]; then
     echo "-----> Creating gradient..."
     grdgradient $grid -A0 -Ne5.0 -G$grid_i -V
   else
     echo "-----> Gradient file exists..."
   fi
```
Run the script

Figure 4.6.3. Colca Canyon topomap based on pp_topo_m2_red.cpt file with shading/illumination supplement and -Ne5.0 option value.

E. Script modifications:

I promised to show the "brain deficiency" effect. Set then -A0 option to -A180 (and -Ne5.0 to -Ne0.6). Like in previous example remove $grid_i file prior to run our script.

   if [ ! -e "$grid_i" ]; then
     echo "-----> Creating gradient..."
     grdgradient $grid -A180 -Ne0.6 -G$grid_i -V
   else
     echo "-----> Gradient file exists..."
   fi

Run the script

Figure 4.6.4. Colca Canyon topomap based on pp_topo_m2_red.cpt file with shading/illumination supplement and with -A180 option.

IMPORTANT NOTE:
Try to look onto canyons and valleys (green), not the peaks (red). Depending to individual brain activity some parts of the picture flip-flops from concave to convex shapes. The illusion strongly hinges on topography as well. Here the central peak is seen properly but the river valleys evince tendency to be seen this-or-that way. Other plots with -A180 option create very strong impression of looking onto inner surface of topo relief. Look onto the two pictures below:


Figure 4.6.5. Poland's SRTM3 tile N49E019 with northern light source. Down on the right lies Tatra Mountains.	Figure 4.6.6. Poland's SRTM3 tile N49E019 with southern light source. Do you see the detestable pale blue tentacles on the top of the picture?

4.7. grdcontour

A. Documentation pages:
     [1] 3.1 Contouring gridded data sets
     [2] O. Annotation of Contours and ``Quoted Lines''
     [3] Cook-book, Example 01: 7.1 The making of contour maps
     [4] Tsunami travel times from the Canary Islands
     [5] GRDCONTOUR man

This one of the most complex GMT program. It has numerous options with dozens of possibilities within the options. True Minotaur's maze. We will use it to plot most simplistic contours as possible without any annotations. The data to plot the contours is taken from CPT file.

Let's introduce Mr. Grdcontour.

   if [ "$contours" == "c1" ]; then
     echo "-----> Creating grdcontours..."
     grdcontour $grid -R$range -J$proj -C$cpt2 -A- -Wc0.05          -K -O >> $psout
   fi

-A - the "-" character standing behind letter "A" defines no annotation.
-W - "0.05" declares pen thickness, "c" regular contours. The thinnest values are better suited for higher dpi (dots per inch) plots. They are also better for complex contour lines found mostly on mountainous areas.

B. Script modifications:

As usual, we must enable the program first.

   # ========== show contours =====================================================
   # contours=c0/c1 where 0=no, 1=yes
   #contours="c0"
   contours="c1"

Put grdgradient to former status (-A180 to -A0) and do not forget to remove $grid_i file.

   if [ ! -e "$grid_i" ]; then
     echo "-----> Creating gradient..."
     grdgradient $grid -A0 -Ne0.6 -G$grid_i -V
   else
     echo "-----> Gradient file exists..."
   fi

Run the script

Figure 4.7.1. Colca Canyon with contour lines.

Figure 4.7.2. Fig 4.7.1 magnified 50 times. The thinner the contour pen the better the picture which has to present details.

4.8. sd-a.cpt from CPT City WWW

A. Documentation pages:
[1] CPT City
[2] Palettes for topography: sd-a

GMT_topo.cpt wasn't designed for every topo plot. There are better cpt files like "sd-a" which can be found on cpt city pages [2].

B. Script modifications:

Change cpt corename to "sd-a"

   # ---------- enter CPT corename ----------------------------
   cpt_name="sd-a"

Set "cpt creation" on (method 2).

   # ========== make new cpt table or use existing one ============================
   # make_cpt=0/1 where 0=no, 1=yes
   #
   make_cpt=1
   if [ "$make_cpt" -eq 1 ]; then
     # ========= select method = 1,2 ==========================
     # cpt_method=1/2
     #
     cpt_method=2
     (...)
   else
   (...)
   fi

Turn off contours.

   # ========== show contours =====================================================
   # contours=c0/c1 where 0=no, 1=yes
   contours="c0"
   #contours="c1"

Rename sd-a.cpt file to GMT_sd-a.cpt
```
   # mv sd-a.cpt GMT_sd-a.cpt
```
Run the script

Figure 4.8.1. Colca Canyon with new color table.

ATTENTION:
The dark spots (areas), which could be seen on another Primer's plots too, indicate areas with no data. I downloaded Colca Canyon tiles several times and got the same results. Nothing can be done.

The End (of The GMT Grdimage Primer)

Przemysław Pawełczyk - Cracow, Poland

Static website generated with Regina-Rexx/Open Object Rexx scripts