Monthly Archives: February 2015

Making simple context maps with GMT

In this Generic Mapping Tools (GMT) tutorial, we will focus on two new things, 1) a different type of map projection and 2) a method to plot data on a map. The map we will make is produced below. I’m going to give some hints about how to make the specific plot, and the script used to produce it is linked below, but I recommend you try and write your own script to make it — this really is the only way to learn.


Final product context map from this tutorial.

1) Another projection type

Let’s start with the following, basic command. It may look similar to a command we used in a previous tutorial, but something should strike you as different.

gmt pscoast -R-126/16/-54/48r -JA-100/30/5i -Xc -Yc -B20a0g20 -Dl -Ggrey -A1000/0/0 >

The -R command has a different style, and there is an A argument following the -J switch we haven’t seen before. Simply put, in the previous post, we were plotting on a Mercator projection, and in this map, we are plotting on an Azimuthal projection. For an Azimuthal projection, it is necessary to specify (in addition to the range) the center of the image, around which all other objects plotted will be scaled. If you need the map to be perfect with respect to minimizing distortion, you will want to be very careful about this, but we don’t really care that much for a context map, so I just picked something close. Also, not that the -R command is specified differently; here we are specifying with a rectangle (lower-left-lon/lower-left-lat/upper-right-lon/upper-right-lat, followed by an r) for the plot area. This is required for an Azimuthal projection, because some parts of the “lines” you specify in the other format for the range may be cut off.

All the projection types with examples can be found here, and the man page description of each type is linked here.

2) Plot data on the map

There are many ways to plot data on a map.  In this case, we want to show where two locations (Lehigh University, and Rice University) are in relation to one another. We’ll also draw a line between them, and create a box that shows the extent of a future tutorial.

Below is the part of the script required to make the two red stars on the map above.

gmt psxy -R -JA-100/30/5i -Sa0.15 -Gred -W0.1p,black << EOF >>
-75.379    40.607167
-95.402778 29.716944

psxy is a GMT tool you will use extensively for plotting anything with point data, e.g. latitude-longitude pairs. You should recognize most of the switches used in this command. The new two switches are the -S and -W; the EOF will be new too.

Now, I want to make two red stars at the two locations for which I have lat-lon data. The -Sa0.15 prints a star of size 0.15 inches. -Gred simply makes it red (similarly to the gray fill we used to make the basemap). The -W switch is a common switch used to denote the pen to use to make a line — in this case it is the outline of the stars, and it is a thin (0.1pt), black pen.

The EOF is sort of a trick you can use when scripting that I find very useful. Essentially, we are using the EOF in the first line to tell the prompt to “hold on”, as we’re going to provide more information. The text from the next line, until a subsequent EOF is read into the command, and then the remainder of the first line is processed. More information here.

Explore a bit!

In all honesty, through this point, you’ve been exposed to most of the basics of building maps with GMT. You now have every tool I used to create the map above. I encourage you to try and write your own script to make the map, but if you get stuck, I’ll provide the script below.

When you build your map, think about layers and what you want to put on top. Don’t forget about using -K and -O properly. Make sure that you include necessary switches, and format arguments correctly. You’ll want to look at the some of the options for specifying pens to make the red and black, dashed/dotted lines.

All the coordinates you will need for the map are below; you’ll need to reformat them a bit to use them in the GMT commands, hint: does latitude or longitude represent an x-coordinate?

Location latitude longitude
Lehigh University 40.607167 -75.379
Rice University 29.716944 -95.402778
Lower left of box 33 -87
Lower right of box 33 -74
Upper left of box 41.5 -87
Upper right of box 41.5 -74

All that said, if you should get stuck and need a hint, or want to see how I did it, you can see my script and final product here.

Good luck!


Click here for the next tutorial about adding real data into your GMT script and making hillshades!

Be inspired while coding! — Matlab script

Here is a joke script I wrote a while ago that returns an inspirational quote when called. I made it as a joke to send to my lab group during a particularly grueling week of coding. You can simply call the function at the beginning of your script to have a quote printed to the stdout or you can wrap it inside a wbar if you want.

Grab the code from my GitHub, here.


example output from the inspire.m script.

Discussion on flexural deformation of crust at passive margins

Recently in class, I led a discussion on the flexural deformation of the crust due to unloading of sediment from mountains, and the concomitant deposition of sediment in passive margin basins. We discussed 3 papers, with a little bit of additional information.

The main point I wanted to make to the group was that flexural deformation is a big player in the long-term decay of an orogen.

I’m putting my powerpoint slides up here as a pdf, and the citations for the discussion are below.

  1. Kusznir et al., 1991. A flexural cantilever model of continental lithosphere extension.
  2. Pazzaglia and Gardner, 1994. Late Cenozoic flexural deformation of the middle US Atlantic passive margin.
  3. Figure, modified from Pazzaglia and Gardner, 2000.
  4. Driscoll and Karner, 1994. Flexural deformation due to Amazon Fan loading.
  5. Rouby et al., 2013. Long term stratigraphic evolution of Atlantic-type passive margins: a numerical approach of interactions between surface processes, flexural isostasy, and 3D thermal subsidence.

Writing scripts to build maps with the Generic Mapping Tools (GMT 5.1.x)

In the last tutorial, we made a simple map, which is what we’ll call “one layer”. What I mean by that is that the entire map was created with one command, and it is just an outline of the eastern United States. In reality, we will want to develop complex maps, with many layers that pointing to and highlight interesting features, places, or anything else you can think of.

Maps with many layers need to be built “from the bottom up”, i.e. later layers will cover the underlying layers (not unlike Steno’s law of superposition!). This is done in the way you probably already figured out, command after command. It may be fine to build a basic map by entering a few commands into the terminal in a specific order, but it quickly becomes cumbersome when maps have many layers and you have to re-enter each command again and again if you want to change anything from an early layer.

For that reason, scripting with GMT is a must. Essentially, we will write a “program” that can be executed line by line by your computer to build the map — this way, each time you want to change something, all you have to do is re-run the program and the map will update! Before we get into the scripting I need to explain the manner in which layering is done with GMT.

GMT uses PostScript to build maps, which means that you need to explicitly state whether you are going to have more layers following the present command, and that you want to overlay the output, not overwrite it. For example take the following three theoretical commands that would, together, produce the theroetically desired map:

gmt commandone -K >
gmt commandtwo -K -O >>
gmt commandthr -O >>

In the first command, we initialize the file and tell the file that there is more PostScript coming (the -K does this). In the second line, we create some more output, overlay (-O >>) the output, and state that more output is to follow; note that the -O switch does not need to be directly adjacent to the >> but both are necessary to properly overlay PostScript. The third line, we again tell the output to be overlayed, but this time, omit the -K switch to finalize the PostScript file after overlaying. I remember the necessary switches for what I want to do by saying in my head that I want to Keep the PostScript file open for more output and/or Overlay the output.

So back to script building, this part will depend on your system, but if you are running on Ubuntu like me, you will need to write the script for a Bash shell. To make a script:

cd into a new directory, wherever you desire, I’m going to make a folder in my home directory called “firstscript”. In there, I’ll make a script called “” and then open it to edit it.

cd ~/firstscript
mkdir firstscript
cd firstscript

Inside the script file, we have to tell the file what shell to execute the script with: make the first line of you script #!/bin/bash. Let’s use the simple command echo to print whatever text we want from the script. We can create comments, or lines of text not to be executed, by beginning the line with a pound sign (#). Insert a comment describing what this script will do above the echo line. This makes the script read:

# this script prints some text
echo "Hello World!"

Now, save the file. We need to change permissions to make it an executable file. Do this with the terminal command chmod 755 Execute the file: ./

You should see your text, “Hello World!” show up in the terminal window. Let’s apply this
to the Generic Mapping Tools with a very basic example.

Below, you will see a script that makes a map with GMT. Copy and paste it into a new shell script you create, or grab it from here (note: make sure the file has executable permissions!).

Make note of the different -K;, -O, and >> positions. I’m not going to detail what each switch and argument does here, but if you want more information check out the manual page for pscoast here.


Click here for the next tutorial

Quick script for connecting to University VPN

I can be pretty lazy at times, so much that I will go out of my way a bit to write a bit of code to simplify my life. I don’t often connect to my university’s VPN, but when I eventually do want to, I can never remember the command needed to do it. So I have to take the time and look it up. Well I cut that out today with a basic little bash script with a name I can remember when I need it. I titled my script “” but anything memorable to you would be fine. The script takes one argument (“on” or “off”) to either connect or disconnect with the proxy. Stick the script in your bin folder and it will execute simply with on.

The script is reproduced below, or you can grab it from here.

# Andrew J. Moodie
# Feb 2015
# easily remembered interface for vpn connection
if [ "$1" = "on" ]
    sudo /usr/sbin/vpnc RiceVPN.conf
elif [ "$1" = "off" ]
    echo " input error -- one argument must be given"
    echo " valid inputs: 'on' or 'off'"

Note that you could further automate this if you wanted to (e.g. auto-enter passwords, connect on startup), but this is a bad idea for both security and for internet speed.