GMT and TCSH Shell Scripting Tutorial

GMT and TCSH shell scripting tutorial

Background

The Generic Mapping Tools (GMT) software package is designed to make various plots using the method of “Paint by numbers”. It is widely used in earth science due to its ease of application, quality of plots, and origin of being written by oceanographers at the University of Hawaii (Paul Wessel and Walter Smith).

The tools take command line options and formatted files to produce post-script output. Post-scripts are then easily converted to other standard images such as .pdf, .png, .tif, .bmp, etc… with the imagemagick library or can usually be directly plotted from a printer or viewed with ghost-viewer.

While the tools can be independently called on a terminal command line, the most typical way of creating plots is by creating a shell script. There are many shells available to use, but the most common in seismology (or at least at Berkeley) is tcsh. The advantage of shell scripting is you use the same syntax interactively that you can also call in the script. This allows you to string together a series of commands which together result in the desired final result.

Before you get started, make sure you have the zip archive, plot_stations_for_orientation.zip which contains this document. You can get it from the IRIS website in the resources for interns section or from a given directory (on board during session).

Objectives

By the end of this session, students should be able to:

• Create simple tcsh scripts
• Create presentation quality location maps
• Overlay various plots to create informative figures

Getting started

Open a terminal to begin. This can be done in different ways varying by system, but the most universal method is to right click somewhere on the desktop and then select “Open/New terminal”. You should get something similar to this:

This is a screenshot of a default terminal from MAC OS X. I’ve customized it slightly to give the time, name of computer, and current directory in the terminal prompt. The grey bar following the $ is a placeholder to start entering commands. You can also see in this image the title bar is “Terminal – bash – 80x8”. This tells me the terminal is currently running the Bourne Again SHell (BASH). If your terminal does not say this, you may be able to find out what shell you’re running by typing:

You can see the terminal is running /bin/bash (which is the full location of the bash shell). You can change to tcsh by entering tcsh at the prompt. All shells have their advantages and disadvantages, but for this tutorial it will be most straightforward if you’re using tcsh.

Having set the shell, we have a little file system work. First create a new directory to work in:

mkdir gmt_tutorial

and copy the files from the zip there:

cp –r plot_stations_for_orientation/* gmt_tutorial

This will make identical copies of all the files and directories from the zip directory to the newly created directory. You can now change directory to the gmt_tutorial directory and list the contents with the following two commands:

cd gmt_tutorial

ls

Now we can start to make a script by opening a new text file with a simple editor.

gedit plot_stations.tcsh &

This will open the program “gedit” and open/create a file “plot_stations.tcsh”. The “&” allows the terminal to continue to accept commands by placing the gedit program in the background. Add the following three lines to the new file so it looks similar to the image below:

#!/bin/tcsh

echo “getting_started”

exit

Select save from the menu bar and then go to the terminal. You can run this script by first changing the permissions:

chmod +x plot_stations.tcsh

(This only ever needs to be done once for each file you make)

And then run the script with

./plot_stations.tcsh

This script will just display some text and exit (note that in the image, I’ve changed to the tcsh shell which has a different prompt on my computer).

Creating a basemap

Creating a basemap is a simple way to generate a new plot file and lay out boundary information. Remove the line “echo “getting_started”” from the script and replace it with:

psbasemap –R0/1/0/1 –JM3.5i –Ba1f1/a1f1nsew >! station_map.ps

save and run the script to produce a new ps file station_map.ps (this could also be run as is on the command line outside of the script)

This command will create a nearly blank image. The gmt tool psbasemap operates with the options “-R0/1/0/1” “-JM3.5i” and “–Ba1f1/a1f1nsew” and outputs post-script file information to what is called standard out. This standard out is redirected with “>!” to create a new file station_map.ps. If you run the command without “>! station_map.ps” on the command line you will see the post-script text put into the terminal.

There is extensive documentation on each of the gmt commands and options. The easiest way to read is by using the man command:

man psbasemap

This will return necessary information to use the command. You can also check out the GMT cookbook at this website:

The power of shell scripts can be more fully utilized by creating some shell variables. If we set some variables before the psbasemap command we can use them repeatedly in a relatively simple manner.

set region = "-R-106.93451/-106.90916/34.06863/34.08084"

set boundaries = "-BSWnea0.005f.1m/.0025f.05m"

set projection = "-JM9.5i"

set misc = "-K –V –X1i –Y1.5i"

set post_script_name = station_map.ps

so the script is now like:

The region variable (-R option) defines the longitude and latitude box

The boundaries variable (-B option) defines how often to plot tick marks, annotations, and gridlines. You can also add titles and axis labels in this option.

The projection (-J option) defines how to map the region of the spherical earth or general data to a flat 2D image. This instance uses a Mercator (M) projection with a horizontal length of 9.5 inches. In this case the vertical axis is defined automatically.

In the misc variable I put various GMT options. –K allows more post script text to be added later. –V turns on verbose mode. –X1i and –Y1.5i move the plot 1 inch in the horizontal direction and 1.5 inches in the vertical respectively.

The post_script_name variable is more commonly chosen to be simply ps or PS. This is a variable commonly used in many scripts. This way you can copy and paste from one script to another and not have to worry about changing the name of the output.

By using shell variables you can add many more commands and re-use the variables. Using long named meaningful names can add significantly to the readability of a script, but its also a good idea to add comments. Any text following the pound sign or hash (#) is considered a comment. Ie:

# this is a comment

Plotting topography

In the directory, there’s a file called PIC.xyz. This contains data in the form of

longitude latitude elevation (m)

in three columns. You can see the limits by using the gmt command minmax:

minmax PIC.xyz

You can change this to a grid file ready for gmt plot with the commands surface or xyz2grd. I prefer surface as it creates a smooth output:

surface PIC.xyz –GPIC.grd $region –I.0001/.0001

In this command we re-use the region variable (which is very handy and will be used significantly). We also add the –GPIC.grd and –I.0001/.0001 options. –G{filename} creates a new grid file named {filename} or in this case PIC.grd. –I{x_increment}/{y_increment} uses an increment of .0001 degrees between each x point and between each y point. This value is a trade off between speed and accuracy. However, there’s a point of diminishing returns where surface is just heavily interpolating.

In this instance surface only needs to be used once. Therefore you could put it on the command line once to create the grid file. Or you can put it the script and comment it out if things get slow.

It is also helpful to have a color palette designed specifically for the data we are about to plot. This can be done with grd2cpt:

set misc = “-D –Z –Chaxby”

grd2cpt PIC.grd $msic >! PIC.cpt

This creates a file (PIC.cpt) with information GMT can use to create a color palette. Like in surface above, the first argument after the command is the filename to use. –D means to set the back and fore-ground values to be the same as the min and maxes. –Z creates a continuous palette, which helps make smoother plots. –Chaxby chooses the “haxby” default palette to be used. GMT has several options you can browse by typing grd2cpt without any arguments.

Finally we are ready to plot with grdimage (after resetting the misc variable)

set misc = “-O –CPIC.cpt –V”

grdimage PIC.grd $projection $region $misc > $post_script_name

Your script should now look like:

This should yield a plot like:

If you wish you can now comment out the surface and grd2cpt lines and associated misc variables by adding the hash sign at the front of the line.

Adding illumination with a gradient file

The GMT tool grdgradient creates a file with the spatial gradients of the 2D grid file created earlier. So add the following commands between the surface command and the grdimage command:

set misc = “-A0 –M –Nt”

grdgradient PIC.grd –GPIC.grad $misc

In this case the –A defines the angle (azimuth) in which the gradient is computed. –M tells grdgradient to compute in meters instead of lat/long. –N is for normalization – and the t option with –N uses a cumulative Cauchy distribution (see man grdgradient)

Now add “-IPIC.grad” to the misc command prior to the grdimage command. Again, after a single run the grdgradient and associated misc can be commented out.

Adding contours

Contours can be added easily once you have a grid with the grdcontour command. But before adding this, update the misc for grdimage to include –K. Also, this is usually better to have after the grdimage command.

set misc = “-C10 –A20f14 –O –V”

grdcontour PIC.grd $region $projection $misc > $post_script_name

The options –C decides the interval of how many meters between each contour. –A defines the meters between each labeled contour and the f option at the end sets the font size. –O means overlay this image over previous images.

Important side note about –O and –K. Commands which write post-script files with either “>! $post_script_name” or “> $post_script_name” usually will need one or both. The first command (in this case psbasemap) needs only –K. The very last command needs to NOT have –K. All others need both –O and –K. If you get frustrated having problems opening the post-script, go back and check the usage of –O and –K.

Overlaying points and text

Plotting x-y data is done with the psxy command. It can receive input from standard in (similar to standard out mentioned before) or from simple data files with x in the first column and y in the second (with the option of some z in the third).

set misc = “-O –K –Sa0.2i -W1p/0 –G0”

psxy $region $projection $misc < END > $post_script_name

-106.91911 34.07344

END

This invokes the psxy command to plot:

a star with radius .2 inches (-Sa0.2i)

a black border around that star with a 1 point black pen (-W1p/0)

fill in that star with black (-G0)

centered at the point -106.91911E longitude 34.07344N latitude

The syntax:

“command < END

….

END”

is the standard in method of saying run this command until you see the text string “END”

We can repeat this for the stations:

set misc = “-O –K –Sa0.2i –W1p/0 –G0/0/255”

psxy $region $projection $misc < END > $post_script_name

-106.92027 34.073959

-106.92382 34.76003

-106.91941 34.075728

END

This will plot all stations as blue (-G0/0/255) stars. If you wished to make each station a different color you could repeat the command with only the latitude and longitude of the station with its own –Gred/green/blue value where each of red/green/blue is a number from 0 to 255. You can alter the symbols by changing –Sa0.2i to something such as –Sc0.1i, which will make circles of radius 0.1 inches.

Text is plotted in a similar fashion with pstext.

set misc = “-O –K –G0”

pstext $region $projection $misc < END > $post_script_name

-106.91911 34.07294 14 0 1 TL PIC

END

In this case the input sent to pstext is a bit more complicated. The single line of standard input is:

Longitude latitude font_size text_rotation font_type text_start_location “text to plot”

In many cases this can be such a pain its easier to just annotate with your final mark up tool such as Powerpoint or Adobe Illustrator.

Adding a scale

Adding a scale bar is good practice. It may not always be necessary on topographic maps as they are usually intuitively colored.

set misc = “-D4.5i/-0.5i/9i/0.25ih –O –K –B10/:meters:”

psscale –CPIC.cpt $misc > $post_script_name

The –D option is required for psscale. It tells the x-center/y-center/x-length/y-length and I’ve append “h” to make it horizontal while default is vertical. The –B flag tells psscale to label every 10 on the scale and put the text “meters” at the end.

Making a location map

You will be presenting your work to a large audience and many people may not know where Socorro, the PIC, or your study area is. Therefore adding a larger scale location map can help orient people. This can be done in four steps: new variables, new basemap on top of the current plot, psxy to on the new plot, and pscoast to show state/country boundaries and lakes/oceans, etc…

set region_inset = “-R-130/-70/26/50”

set projection_inset = “-JM3i”

set misc = “-O –K –G255 –Bnews”

psbasemap $region_inset $projection_inset $misc > $post_script_name

Now reset misc and do a psxy

set misc = “-O –K –Sa0.2 –W1p/0 –G0”

psxy $region_inset $projection_inset $misc < END >$post_script_name

-106.903 34.0663

END

Finally we use pscoast to draw state boundaries

set misc = “-Dl –Na –W1p/0 –O –K”

pscoast $region_inset $projection_inset $misc > $post_script_name

In this case the –Dl selects a “low” resolution dataset. –Na means to draw all political boundaries.

At this point your map should look similar to below:

Plotting a SAC waveform (ADVANCED)

You’ve now created a presentation quality map, which will be useful to orient others to your study location. You can also use a similar process to plot any type of xyz data (such as a tomographic model) or you can use the tools to plot any number of figures. This section goes a little bit further into preparing programs and then uses a simple c program to prepare data and send the output to psxy to plot a seismic waveform at the top of the map.

The first step is to compile the program “sac2stdout”. There are two versions of this code currently available because SAC data files have slightly different formats between 32bit and 64bit systems. The SAC data file IIPC.PC.LHZ.SAC included with this was created on a 64bit system. Therefore you should install the 64bit version of sac2stdout:

cd sac2stdout

cp Makefile.64 Makefile

make

cp sac2stdout ..

cd ..

This will create a copy of the exectuable function “sac2stdout” in the same working directory you’ve been plotting in. Generally a best practice would be to move the exectuable to a standard location (such as ~/bin) within your linux/unix search path so you only have one version of the executable.

Now we’re going to set up some variables. I often make use of the syntax:

command `other_unix_command`

This will first execute whatever is in the ` ` and then return the value to the outside command similar to parentheses in math.

set sac_data = IIPC.PC.LHZ.SAC

set misc = “-O –K –G255 –Bnews –Y6i”

set projection_sac_data = “-JX9.5i/.7i”

set region_sac_data = `./sac2stdout $sac_data | minmax –I1”

Several options and syntax were introduced here:

-Y6i: Move the plot axis 6 inches vertically

-JX9.5i/.7i: The projection is flat x-y with an x axis of 9.5 inches and y axis of 0.7 inches

|: the “pipe” redirects output from the first command to become the input to the second command.

-I1 on minmax returns the minimum and maximum of the input to the form:

-Rx_min/x_max/y_min/y_max

Now we can create a basemap

psbasemap $region_sac_data $projection_sac_data $misc > $post_script_name

We can also use minmax with the –C option to return columns with each of the limits:

set limits = `sac2stdout $sac_data | minmax –C`

This command creates a variable array called limits with 4 entries. We can then extract them with the echo command:

set xmin = `echo $limits[1]`

set xmax = `echo $limits[2]`

set ymin = `echo $limits[3]`

set ymax = `echo $limits[4]`

With those variables in hand we can call on awk (which is another programming language in itself) to do some simple math and calculate reasonable spacing for the tick marks on our x-y seismic data plot.

set xtick = `echo $xmax $xmin | awk ‘{print ($1-$2) / 20}’`

set ytick =`echo $ymax $ymin | awk ‘{print ($1-$2) / 5}’`

Finally we can setup a new misc variable and plot

set misc = “-O –W2p/20/120/200 –Ba${xtick}f${xtick}/a${ytick}f${ytick}nSeW”

sac2stdout $sac_data | psxy $region_sac_data $projection_sac_data $misc > $post_script_name

This will be the last thing we plot in this exercise so the misc variable must contain –O but it must not contain –K.

Finishing and exporting

If the ImageMagick library is properly installed, two common conversion tools are available. ps2pdf will convert your freshly created post-script file to a more portable pdf file (syntax: ps2pdf file.ps). convert will freely convert any input image format to any desired output. It works on the last letters after the last “.” So be careful about naming files. To make a .png you can add:

set png_name = station_map.png

convert $post_script_name $png_name

Finally, I like to finish my scripts with a message to the user.

echo “$post_script_name created and converted to $png_name”

Your final image should be similar to: