Shell scripting: Generic Mapping Tools (GMT)

The way that user commands are able to interact with the operating system of the computer is done via shells. There are different Linux shells, and these can be viewed as different ways to interact with the computer. The commands are the same, but the interaction with the system can be slightly different. Shells also allow users to write scripts, or programs, that will have different syntax depending on the shell.

A. Basics system scripting

Three of the most common shells are the bourne shell (sh), the c-shell (csh), and the bourne again shell (bash). Aside from scripting, these become important when setting various parameters that impact our interaction with the operating system. As mentioned earlier, many configuration files in Linux start with a period (.), and the shell configurations are one such example.

The shell startup/configuration files optionally contain environmental variables, such as terminal type, paths to look for files, prompt display, etc. These environmental variables can be set interactively or automatically at login, and the variables/settings are maintained throughout the login (shell) session. The command env will show all of these that are currently set, while the echo command will display single ones:

could return something like this:

XDG_SESSION_ID=5

TERM=xterm-256color

SHELL=/bin/bash

SSH_CLIENT=128.171.159.89 63503 22

SSH_TTY=/dev/pts/8

USER=installation

MAIL=/var/mail/installation

PATH=/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games

PWD=/homelocal/installation

LANG=en_US.UTF-8

SHLVL=1

HOME=/homelocal/installation

LOGNAME=installation

SSH_CONNECTION=128.171.159.89 63503 128.171.153.84 22

LESSOPEN=| /usr/bin/lesspipe %s

XDG_RUNTIME_DIR=/run/user/1000

This shows that the terminal is currently set to an Xterminal (TERM=xterm-256color), the current shell is bourne again shell (SHELL=/bin/bash), etc. To display a single value, use the echo command, followed by a dollar sign ($) and the environmental variable:

installation@kilo:~$ echo $SHELL

will return /bin/bash. To change the value of an environmental variable in csh:

installation@kilo:~$ setenv TERM vt100

will change the terminal type to a vt100 (instead of an Xterminal). To do the same thing in bash:

installation@kilo:~$ TERM=vt100

installation@kilo:~$ export TERM

The files .cshrc (for csh) and .bashrc (for bash) are configuration files in the user’s home directory. These control all the environmental variables upon login, so they represent a convenient way to set variables automatically. Additionally, these files can contain aliases, or shorthand notation for various commands. Equally important, the dot-startup files are where the pathnames are set; these control where the system will “look for” commands.

Here is an example .cshrc file:

# jimp’s .cshrc file for the MacBook Air

# Set path

set path=(. /usr/local/bin /usr/sbin /opt/local/bin /usr/bin /bin /sbin /Users/jimp/Bin )

# Set terminal control stuff

stty kill ‘^U’ intr ‘^C’ echoe

# Set initial stuff

set history = 100

set filec

limit coredumpsize 0

# Define prompt

if ( (! $?ENVONLY) && $?prompt ) then

if ( -o /bin/su ) then

set prompt=”iMac \!# “

else

set prompt=”iMac \!% “

endif

endif

# GMT stuff

setenv GMTHOME /usr/local/gmt-4.5.6/data

setenv NETCDFHOME /usr/local/netcdf-3.6.3

setenv GMT_DATADIR /usr/local/gmt-4.5.6/data/dbase

# Define aliases

alias ace ‘xmgr -geometry 1000x850+25+50 -noask -nxy’

alias cl ‘clear’

alias ll ‘ls –laF’

The lines starting with a hash (#) are comments and not read by the system. The first line sets the path, so the operating system will look for commands when entered by searching the current directory (.), /usr/local/bin, /usr/sbin, etc. The next line defines how to cancel and kill commands. The next two sets of commands control how many commands are saved and what the prompt will look like. Following this are some program-specific environmental variables (in this case for a program called GMT, or Generic Mapping Tools). The final set are aliases. The syntax here is for a shorthand notation, e.g., ll, followed by the full command, ls –laF. Therefore, instead of entering “ls –laF” all the time this user can just enter ll.

An example .bashrc file would be similar:

# ~/.bashrc: executed by bash(1) for non-login shells.

# for setting history length see HISTSIZE and HISTFILESIZE in bash(1)

HISTSIZE=1000

# set variable identifying the chroot you work in (used in the prompt below)

if [ -z “${debian_chroot:-}” ] && [ -r /etc/debian_chroot ]; then

debian_chroot=$(cat /etc/debian_chroot)

fi

# set a fancy prompt (non-color, unless we know we “want” color)

case “$TERM” in

xterm-color) color_prompt=yes;;

esac

# uncomment for a colored prompt, if the terminal has the capability; turned

# off by default to not distract the user: the focus in a terminal window

# should be on the output of commands, not on the prompt

force_color_prompt=yes

if [ -n “$force_color_prompt” ]; then

if [ -x /usr/bin/tput ] && tput setaf 1 >&/dev/null; then

# We have color support; assume it’s compliant with Ecma-48

# (ISO/IEC-6429). (Lack of such support is extremely rare, and such

# a case would tend to support setf rather than setaf.)

color_prompt=yes

else

color_prompt=

fi

fi

# enable color support of ls and also add handy aliases

if [ -x /usr/bin/dircolors ]; then

test -r ~/.dircolors && eval “$(dircolors -b ~/.dircolors)” || eval “$(dircolors -b)”

alias ls=’ls –color=auto’

#alias dir=’dir –color=auto’

#alias vdir=’vdir –color=auto’

alias grep=’grep –color=auto’

alias fgrep=’fgrep –color=auto’

alias egrep=’egrep –color=auto’

fi

# some more ls aliases

alias ll=’ls -la’

alias l=’ls -aF’

#alias l=’ls -CF’

More information can be found on the man pages.

We will now attempt to make our first shell script. One way to think about this is whatever Linux commands we enter into the terminal can be placed into a file, then we simply run (execute) that file. For example, the command “ls -l” will list all the files in a directory. We could instead create a file, e.g., “list_files”, include the ls command, then run the file. This is a very simple example, but hopefully it illustrates the point.

First, we use leafpad to create our script file. Since it’s a shell script, let’s use the convention of appending .s to the end, i.e., create a file called list_files.s. The contents of the file would be simply:

installation@kilo:~$ ls -laF

We could now run this file, causing it in turn to execute that command, and output a listing. First, recall that in order for a file to be executable, it needs to have that set (see chmod discussion above). So, first we make it executable:

installation@kilo:~$ chmod 755 list_files.s

then we can run it:

installation@kilo:~$ ./list_files.s

Two cautionary notes here. First, the file needs to be executable to run. If it’s not, you’ll see an error like this:

installation@kilo:~$ ./list_files.s

-bash: ./list_files.s: Permission denied

Second, we need to have the script in our path. Note above we use the explicit, relative path (./ meaning “in this current directory, find the file list_files.s”). If your current working directory is not in your path, and you try run the script without specify a path, you may get:

installation@kilo:~$ list_files.s

list_files.s: command not found

These two issues, file permission and path, are usually the cause of many script failures.

Now that we have our file in the correct form, and we know how to run it, we can spruce it up a bit:

#!/bin/bash

# This script will output the long file listing to a file called “myfiles.txt”

rm myfiles.txt

ls -laF > myfiles.txt

Note here we’ve added a comment (preceded with the hash, #), which is always good practice, and made use of the file output redirect (>). We have also added a command to first remove the file.

We will not be doing much shell scripting in class, but rather some Matlab and python scripting (more on these later). One nice example of the utility of shell scripting is the Generic Mapping Tools (GMT). These are a set of programs that construct geographic maps, and it’s much easier and more efficient to run these in a script instead of interactively at the terminal (although you could do that).

B. The Generic Mapping Tools (GMT)

We will try and bring together what we’ve learned about Linux commands and scripting by using the Generic Mapping Tools (GMT). GMT was developed by Paul Wessel (here in SOEST), and is available on the SOEST site (http://www.soest.hawaii.edu Research data access). GMT is actually a series of programs that are used to make maps and other figures from input data files. Note that appears complex at first, but it is very flexible and powerful. You can download it to your PC or Mac (note version differences). The output of GMT is postscript, a format used for printing.

Like other Linux commands, GMT commands are entered in the form of some command followed by a number of options preceded with a dash (-), like we saw with ls –l. If you recall our discussion about file I/O, GMT outputs to standard out, i.e., the screen. This is not very helpful, so we will redirect standard out to a file.

The first command we will use is pscoast to draw a map. This command has many options, but we will use the bare minimum, i.e., the projection (JM) the range (R) and the pen width/color:

/usr/local/gmt/bin/pscoast -R110/300/-30/60 –JM25 -W1/0/0/0

Note that the output comes to the screen. Again, not helpful, so we redirect to a file:

/usr/local/gmt/bin/pscoast -R110/300/-30/60 -W1/0/0/0 -G255/0/255 -JM25 –Df -P > sample_map.ps

We can now use the command “evince” to view the PostScript file (note that evince comes with our Linux ubuntu installation, it’s not a generic Linux command). You should see something like this:

image0

Notes:

  1. the “path” to the GMT programs needs to be included, otherwise the operating system can’t “find” the program, so instead of just pscoast, we use /usr/local/gmt/bin/pscoast. Again, this is called an absolute path.

  2. the program pscoast takes options like all Linux commands with a dash (-), in this case we specify

    -R110/300/-30/60 longitude range (110 to 300) and latitude range (-30 to 60)

    -W1/0/0/0 pen width (1) and pen color (0/0/0)

    -G255/0/255 continent shading color

    -JM25 map projection and size (mercator, 25cm wide)

    -P portrait plot

GMT specifies colors with an “RGB” syntax, or red-green-blue values from 0 (none) to 255 (full). Also, recall that GMT programs generate PostScript, and the default is to stream the result to the terminal/screen. This is typically not very helpful, so instead we take advance of the Linux ability to redirect output to a file. Here, we specify the file “sample_map.ps”. It is good convention to note PostScript files with the extension “.ps”. The redirect is done using the greater than sign (>). This tells the operating system to redirect the output of the command “pscoast” to the file “sample_map.ps”.

Notes:

  1. in some cases, the “>” will overwrite an existing file while in others it will append the output to the file. In most cases, “>!” is used to overwrite files (e.g., if the file exists, replace it), and “>>” is used to append (e.g., if the file exists, tack this stuff at the end). Like with path setting, this is done in your .bashrc or .tcsh file.

We now want to expand upon our map-making abilities and include things like a border with lat/lon labels, maybe some symbols and lines, etc. This is done using different GMT commands, specifically “psbasemap“, “pscoast“, “psxy“, and “grdraster/grdcontour“.

Since we are using multiple GMT commands, in succession, to write to the same file, we need to add extra options. This is one common cause for error in GMT use.

The option “-O” means “overlay this feature onto my existing plot”. If you don’t have this “-O”, the result from your new commands will cover the previous one. Next, the option “-K” means “more drawing is coming, so don’t close the output file”. In this case, if you don’t include the “-K”, the PostScript file will get closing information added, and the file will appear “done”.

Here’s an example. We will take our example above and add a frame around the plot using “psbasemap”:

/usr/local/gmt/bin/psbasemap -R110/300/-30/60 -JM25 -B20/10 -K > sample_map2.ps
/usr/local/gmt/bin/pscoast -R -J -W1/0/0/0 -G255/0/255 -Df -O >> sample_map2.ps

Note the use of “>” to redirect the output to a file called sample_map2.ps, and then the use of “>>” to append new stuff to this file (first and second lines). Further, note the following GMT-related things:

  1. we include the “-K” on the first line to indicate “more follows” and the “-O” in the second line to indicate “overlay this”, but we don’t include a “-K” since we have nothing else to add

  2. we introduced the new option “-B” to add tick marks on the border, spaced 20-degrees along the x-axis (longitude) and 10-degress along the y-axis.

  3. GMT allows you to keep earlier definitions by simply not restating them. So, for example, pscoast requires “-R”. Since we already defined the ranges in the first line, we simply include “-R” to imply “keep the same definitions”. Note you don’t need to do this. You could, for example, include the pscoast on a totally different lat/lon range as the border (try do this).

You should try with different lat/lon ranges, colors, etc. to get comfortable with the GMT syntax. You will notice that typing these commands out repeatedly becomes tiresome. Instead, we want to include these in a script and then just run that. One way to thing about this is we will write a program that includes these commands, then just run that program. In this case, the program is not written in some language like FORTRAN or python or Java, but is an ascii file containing other commands, i.e., a shell scipt.

To do this, we first need to learn how to create a text file. You can use various command line editors (e.g., “vi”) but these are typically complex. You could also use something like MS Word, and save the result as a text file, but this could take a while. Instead we will use a light-weight editor called “leafpad“. You can find this via the pull-down menus or put entering the command in a terminal window. As our first example, just copy/paste the two commands that make a map into a file and then save it as a file. For example, save the file as “make_map.s“. We now want to run this command from the terminal. Does this work? If not, why not?

We now our first (hopefully working) script. Notes:

  1. we add comments to a script with a hash (#) sign. It’s always good practice to include comments in scripts and programs

  2. it’s also common, though not essential, to name shell scripts like this with a .s extension (later we’ll do matlab scripts with a .m).

  3. there are actually two ways to run it; another way is to use the bash shell

Ok, now that we have our script, let’s experiment with more GMT commands. We now will include “psxy” to dsplay symbols at points and “grdraster” to add bathymetry. Note that psxy requires a file containing the lon/lat pairs for the locations where we want to print things. We can use leafpad to create this file, save it, then access it in our script:

/usr/local/gmt/bin/psbasemap -R110/300/-30/60 -JM25 -B20/10 -K > sample_map3.ps
/usr/local/gmt/bin/pscoast -R -J -W1/0/0/0 -G255/0/255 -Df -O -K >> sample_map3.ps
/usr/local/gmt/bin/psxy places.dat -R -J -Sa0.5 -O >> sample_map3.ps

and get this:

image1

Small side note about images… PostScript works great for printing, not so good for web/wiki. Instead it’s better to use gif, jpeg or png images. Fortunately there is a Linux tool to convert between the two. It’s called imagemagick (http://www.imagemagick.org/script/index.php) and the specific program is “convert”. More examples are here: http://gmt.soest.hawaii.edu/doc/5.1.0/Gallery.html

As a quick review, we’ve gone over some basic Linux commands and then started making plots with the Generic Mapping Tools (GMT). Since GMT is essentially a set of programs, it lends itself to scripting, so we use GMT as an example for using scripts. First we used a few GMT commands including “psbasemap” to draw a map border (and other things) and “pscoast” to draw the coastlines. Since GMT typically writes output to the screen (standard output), and we want a file, we use the Linux convention of redirecting output to a file with “>”. Further, we added these commands to a file, then “ran” that file, i.e., we executed a script. Recall we had to change the file permissions on our script to make it “executable.” Here’s a quick summary of important concepts:

  • GMT commands will make “layers” or maps on top of each other. The option “-O” is needed to indicate the current layer should be an overlay on the previous one and not cover it up.

  • PostScript files require information at the end to indicate the termination of the file. In GMT this is done by default, but usually we want to include more on a plot. To do this, we use the “-K” option. Think of this as “more information to be displayed”. If you don’t include this, nothing else will be plotted.

  • To summarize the previous two points, usually the first line in a GMT file will have a -K (more to come), and the last line will have a -O (overlay) but not a -K, and all lines in between will have both -O and -K

  • On the Linux side, we use “chmod” to change the permissions on a file to make it executable (see man pages for details).

  • Also on the Linux side, we need to be careful on paths, or the directory structure to our script. The operating system needs to know where to find the program (or script) to run it. This is done by specifying the path, either “absolute” or “relative”. Here are two examples. First, this will run the script “make_map.s” in the directory /home/iniki/ges/jim/scripts/class1 (absolute path):

    installation@kilo:~$ /home/iniki/ges/jim/scripts/class1/make_map.s

    Second, this will run the script make_map.s in the current working directory (relative path):

    installation@kilo:~$ ./make_map.s

To now carry on with GMT and scripting. Note that every line in a script will be run like it were entered as a command at the prompt. Sometimes we don’t want this, for example to put notes or comments into a script. To do this, put a “#” symbol. This will tell the operating system/shell to not try and execute this line. In addition, let’s add a line to automatically remove the output file if there, and another line to automatically display the file using evince. Here’s an example from last class with comments added:

# My first GMT script

#

# Written on January 27th 2019

# 1. remove the output file if it exists

rm sample_map3.ps

# 2. draw a map frame; remember the -R is for domain

# (lon/lat), -J is for projection and -B is for tick

# marks

/usr/local/gmt/bin/psbasemap -R110/300/-30/60 -JM25 -B20/10 -K > sample_map3.ps

# 3. draw a coastline; shade the land purple (255/0/255)

# and use full resolution (-Df) /usr/local/gmt/bin/pscoast -R -J -W1/0/0/0 -G255/0/255 -Df -O >> sample_map3.ps

# 4. Display the output using evince

evince sample_map3.ps

We will now use two new GMT commands, “psxy” to plot text/symbols and “grdcontour” to plot contour levels. Both these require ascii files with information such as the location for the text/symbol (psxy) or the levels to contour (grdcontour). So, the first step will be to create an ascii text file with the lon, lat of the points we want to display. As an example, create a file called “places.dat” using the text editor “leafpad”. This file will have two columns, longitude and latitude, for every point. If we have one line with “202 24” our location will be near Hawaii. To plot this:

/usr/local/gmt/bin/psxy places.dat -R -J -Sa0.5 -O >> sample_map3.ps

Note that we specify our file places.dat and use -Sa0.5 to specify a star (Sa) that is 0.5 centimeters across. Remember the -O for overlay mode.

Next, we will add bathymetry contours as lines (you could also do these as colors; check the GMT site for examples). This is a two-step process. First, we need to create a GMT data file that contains the bathymetry. You can either use your own file or use one of the GMT-supplied files. In this case we will use an internal file, data set number 3 (see GMT notes for definitions of these).

/usr/local/gmt/bin/grdraster 3 -R -Gtopo.grd

Should create a file called “topo.grd” that can be read by GMT. This file will contain the bathymetric data from GMT dataset number 3 (which happens to be from GeoSat). Next, we use “grdcontour” to plot the depth. The file that contains the levels we wish to show is called “levels.dat” and we create this using a text editor (again, leafpad). This file will have three lines, each with one value, -500, -1000 and -2000 to plot these three isobaths. Our revised GMT script is like this:

# My first GMT script

#

# Written on January 27th 2019

# 1. remove the output file if it exists

rm sample_map4.ps

# 2. draw a map frame; remember the -R is for domain (lon/lat),

# -J is for projection and -B is for tick marks

/usr/local/gmt/bin/psbasemap -R110/300/-30/60 -JM25 -B20/10 -K > sample_map4.ps

# 3. draw a coastline; shade the land purple (255/0/255) and use

# full resolution (-Df)

/usr/local/gmt/bin/pscoast -R -J -W1/0/0/0 -G255/0/255 -Df -K -O >> sample_map4.ps

# 4. Now plot a green star (-Sa) near Hawaii (note the file

# “places.dat” has the coordinates

/usr/local/gmt/bin/psxy places.dat -R -J -W4/0/255/0 -Sa0.5 -O -K >> sample_map4.ps

# 5. Finally, we create a topography file (topo.dat) and then plot

# the contours

/usr/local/gmt/bin/grdraster 3 -R -Gtopo.grd

/usr/local/gmt/bin/grdcontour topo.grd -R -J -O -Clevels.dat >> sample_map4.ps

# 6. Convert the output from PostScript to jpg to post on the wiki

convert sample_map4.ps sample_map4.jpg

# 7. display the output using evince

evince sample_map4.ps

image2

To demonstrate the power of using a script, we can easily modify the plot to center on Hawaii but just changing the -R range to get:

image3

In this final example, we can make use of all the more common GMT commands. We would like to make an annotated map with ocean depth color-shaded (note we could shade any field, but GMT has “internal” bathymetry” files). We will also include some lines connecting a few places. We will need to make a few text files that will specify things like text, color levels and values, and we will define these as we go.

First, we use our standard basemap call, noting that we need a -K (more postscript coming) and a single redirect (>) so as not to append to existing file. We also use the -B option to specify the tick marks but also the plot title (between colons and preceded with a dot):

/usr/local/gmt/bin/psbasemap -R110/300/-30/60 -JM25 -B20/10:.”My Excellent Map”: -K > gmt_sample_map.ps

Next, we draw our map. In this example we have a red coastline (-W), grey land (-G) and cyan ocean (-S). We have several notes here as well:

  • Add -O (overlay) and -K (more coming later)

  • Add two redirects (>>) to append to previous drawing

  • We use the more simple -R -J to use the original values (can change if wanted however)

  • Specify the intermediate coastline resolution (-Di) to make the file size smaller

    /usr/local/gmt/bin/pscoast -R -J -W1/255/0/0 -Di -G125/125/125 -O -K >> gmt_sample_map.ps

We now want to plot a symbol (star in this example) at certain locations. To do this we specify the locations by the longitude/latitude in a file (in this case symbol_locs.dat):

# Honolulu

-158.0 21.4

# Tokyo

139.7 35.7

# Los Angeles

-118.2 34.0

# Samoa

-172.1 -13.8

# Seattle

-122.3 47.6

and then access this file with psxy:

/usr/local/gmt/bin/psxy symbol_locs.dat -R -J -Sa0.5 -W0/0/255 -G0/255/0 -O -K >> gmt_sample_map.ps

The next step is to label these locations with pxtext. The input file here is somewhat tricky (i.e., specific). For example, the font size should be given as point-size and as an integer.

#x, y, size, angle, fontno, justify, text

-158.0 21.4 14 0 1 1 Honolulu

139.7 35.7 14 0 1 1 Tokyo

-118.2 34.0 14 0 1 1 Los Angeles

-172.1 -13.8 14 0 1 1 Samoa

-122.3 47.6 14 0 1 1 Seattle

The angle is in degrees, while the justification is a numeric code (see man page). The text command (pstext) takes a font color with -G (alternately a background with -W):

/usr/local/gmt/bin/pstext label_locs.dat -R -J -G0/0/255 -O -K >> gmt_sample_map.ps

We will now draw a line connecting all points to Honolulu. This is done by specifying end points of the line and segments separated with a “>” in yet another file (line_locs.dat):

-158.0 21.4

139.7 35.7

>

-158.0 21.4

-118.2 34.0

>

-158.0 21.4

-172.1 -13.8

>

-158.0 21.4

-122.3 47.6

and

/usr/local/gmt/bin/psxy line_locs.dat -M -R -J -W4/255/0/255 -O -K >> gmt_sample_map.ps

Finally, we will add some color-shaded bathymetry. To do this, we will take advantage of the built-in GMT bathymetry files. In general, GMT will color-shade “raster” or gridded files of a specific format. The command grdraster will make such a file either from supplied data or the internal ones. The definition of these depends on the machine installation, and for our class the GMT install was done in /usr/local/gmt. Within this directory is a sub-directory called share/dbase. Here are the bathymetry files, along with a text file explaining each one (grdraster.info):

1 “Sea floor age from Muller et al., 1998” “Ma” -R0/360/-72/90 -I6m G i 0.01 0 999 age_1.5.i2 L

3 “ETOPO5 global topography” “m” -R0/359:55/-90/90 -I5m G i 1 0 none etopo5.i2 L

2 “US Elevations from USGS” “m” -R234/294/24/50 -I0.5m P i 1 0 -9999 us_topo30s.i2 L

4 “1 min Geoware blend Observed/Predicted Topo” “m” -R0/360/-90/90 -I1m P i 1 0 32767

world_relief_1m.i2 L

5 “2 min Geoware blend Observed/Predicted Topo” “m” -R0/360/-90/90 -I2m P i 1 0 32767

world_relief_2m.i2 L

6 “30 sec Geoware blend Observed/Predicted Topo” “m” -R0/360/-90/90 -I30c P i 1 0 32767

world_relief_30s.i2 L

The second file, curiously indexed as “3”, contains the ETOPO-5 data. To access this, we first make a raster (grid) file:

/usr/local/gmt/bin/grdraster 3 -R -Gbathy.grd

This will create a file called “bathy.grd”. We then display this using custom colors:

/usr/local/gmt/bin/grdimage bathy.grd -Cbathy_colors.dat -R -J -O -K >> gmt_sample_map.ps

The file bathy_colors.dat is another text file that contains the end-points of a range and the color for these two ends. Below is an example:

-10000 84 0 84 -5000 125 0 255

-5000 125 0 255 -1000 0 0 255

-1000 0 0 255 -100 0 145 255

-100 0 145 255 -10 0 255 255

-10 0 255 255 0 0 255 255

The first line specifies that depth of 10,000 meters (negative meaning depth below sea level) will be colored red=84, green=0 and blue=84 (purple-ish), depth of 5,000 meters will be red=125, green=0, blue=255 (more blue-ish), and values in between will be interpolated. This continues in a similar way on the subsequent lines.

Our final step is to provide a colorbar. This is done using psscale and another text file. The text file contains the color level to be annotated and the size of the block. This can be useful in the present example, since the color-scale is non-linear. In other words, we want the color-scale to show the difference between -10000 and -5000 similar to -10 and 0 (one being 5000 wide the other only ten). The colorbar_labels.dat file would like this:

0.8 -10000

0.8 -5000

0.8 -1000

0.8 -100

0.8 -10

This forces each color “block” to be 0.8 wide. Note here that this file needs to have the same levels as the colors in the plot (bathy_colors.dat). The final command is then:

/usr/local/gmt/bin/psscale -Cbathy_colors.dat -Zcolorbar_labels.dat

-D5.00/18.0/10.0/0.50h -O >> gmt_sample.ps

Again, note here we have a -O but no -K (important). Since the bathymetry was drawn last, it will obscure the rest of the plot, so all we need do is switch the order. A final script would look like this:

# A sample GMT script to plot colored bathymetry, text and such

# Jim P 01-23-2019

#

# 1. Draw the map frame:

# -R for lon/lat range

# -JM25 for Mercator projection 25cm wide

# -B20/10 for tick marks in lon/lat direction (and title in :.””.)

# -K “more postscript coming

/usr/local/gmt/bin/psbasemap -R110/300/-30/60 -JM25 -B20/10:.”My Excellent Map”: -K > gmt_sample.ps

# 2. Draw the coast

# -W1/255/0/0: draw shoreline colored red

# -Di: use intermediate coastline

# -G125/125/125: shade land grey

# -O -K: overlay and more postscript coming

/usr/local/gmt/bin/pscoast -R -J -W1/255/0/0 -Di -G125/125/125 -O -K >> gmt_sample.ps

# 3. Make a gridded bathymetry file from the ETOPO-5 data set (3)

# and shade the depths with values given in “bathy_colors.dat”

# -O -K: overlay and more postscript coming

/usr/local/gmt/bin/grdraster 3 -R -Gbathy.grd

/usr/local/gmt/bin/grdimage bathy.grd -Cbathy_colors.dat -R -J -O -K >> gmt_sample.ps

# 4. Plot symbols at places given in symbol_locs.dat (lon/lat)

# -Sa0.5: Symbol of star (a) of size 0.25 cm

# -W0/0/255: border is blue

# -G0/255/0: shading is green

# -O -K: overlay and more postscript coming

/usr/local/gmt/bin/psxy symbol_locs.dat -R -J -Sa0.5 -W0/0/255 -G0/255/0 -O -K >> gmt_sample.ps

# 5. Draw text at symbol locations (label_locs.dat)

# -O -K: overlay and more postscript coming

/usr/local/gmt/bin/pstext label_locs.dat -R -J -G0/0/255 -O -K >> gmt_sample.ps

# 6. Draw lines between points (-M multisegment)

# -W4/255/0/255: draw thick line colored purple

# -O -K: overlay and more postscript coming

/usr/local/gmt/bin/psxy line_locs.dat -M -R -J -W4/255/0/255 -O -K >> gmt_sample.ps

# 7. Draw color bar

# -O: overlay only

# -D5.00/18.0/10.0/0.50h: colorbar at x=5, y=18, width=10cm,

# height=0.5cm (and horizontal)

/usr/local/gmt/bin/psscale -Cbathy_colors.dat -Zcolorbar_labels.dat -D5.00/18.0/10.0/0.50h -O >> gmt_sample.ps

/usr/local/bin/xv gmt_sample.ps

GMT is extremely powerful for making geographic-based plots, but it can also be used to make line and contour plots. One thing GMT cannot do is “create” or modify data. For that we will need more sophisticated tools. We will next look into python, and later Matlab and QGIS.