These are a Few of My Favorite Things - Mapping with GMT Part 1
So you want to make a map ... a pretty map, a useful map, a map that helps people understand your data. Where do you go? Matlab? Python? ArcGIS? All of these platforms have pros and cons (ok, Matlab's plotting pros are pretty slim but it'll get the job done) but for me there is no comparison with the hardcore, no frills utility aptly named "Generic Mapping Tools".
GMT (current version 5.2.1) is an open source compilation of 80+ command line tools used to produce Cartesian and geographic post-script figures. The beautiful power of GMT is in its simplicity. With just a few quick commands or even just one you can make a decent looking map. Add on a couple more commands and you can create an impeccable map that displays your data in whatever way you want! In the following post I will share with you the beginnings of what I know, love, and have learned better of about GMT.
First thing first you need to think about what you want to show with your map. What relationships are you trying to highlight? What point do you want driven home? Without a doubt you likely will want to include some kind of elevation/bathymetry in the background. If you're looking for elevation grids especially on a continental to global scale there's no better place to look than ETOPO1. Long story short, this is a 1-arc minute global relief model that includes both elevation and oceanic bathymetry. I created the above map with ETOPO1 and you can see the amazing detail included in this dataset. If you are interested in a smaller region I also suggest checking out the Shuttle Radar Topography Mission (SRTM) datasets which offer incredible 30m resolution across many parts of the globe.
Once you've got your datasets in hand (topography grids, text files of volcano locations, ...) it's time to actually start coding with GMT - the whole point of this post! Because GMT is a set of command line operations it can be run solely from the terminal but that is by no means an effective tac to take. Instead I strongly encourage you to embrace either bash or csh scripting. Using either of these languages you can construct a script composed of all the GMT commands you need to create your perfect map.
My regular GMT scripts follow a fairly regular format:
Let's walk through each of the above components. The first block of text sets a bunch of qualitative defaults for your figure with gmtset like page size, frame width (the boundary around the map), and the format of the latitude and longitude labels. The number of defaults you can set is EXTENSIVE! I would suggest start off with just a few and then add as you need. I commonly use the same gmt defaults in every script I use.
The next section sets up several environmental variables that make the actual scripting easier when you need to use the same values in multiple commands. I commonly set a variable for the name of the postscript file (PS1), the names of elevation grids and illumination files (GRID & INT), the map scale and area (SCALE & REGION), and the color file (CPT). Note that you will be forcing the output of all commands to your postscript file so if nothing else SET a variable for the name of your postscript file!
The chunk of text setting variables for colors is fairly self explanatory and if you haven't seen my post about colors then go there now!
Now we finally get to the meat of the script where we actually call the gmt commands that plot things. First is the command grdgradient which we use to create the illumination file for our elevation grid. Illumination files allow us to shade topography (think hillshade) and should always be used unless you have a really good reason not to. I've commented out this first command as you only need to do it once to create the file and then you wont need to call it again.
Next we plot the elevation grid with shading from the illumination file via grdimage. Because this is the first command to write to our postscript file we need to set some map defaults like the map scale and projection (-J flag) and the region (-R flag). We set the labeling and tick marks for the lat/lon axes with the -B flag and the color scale to use with the -C flag. Also because we are not done plotting we include the -K flag which tells GMT that more commands are upcoming. Finally to make sure we don't plot on top of a pre-existing file we use a single > symbol to send the output to our new postscript. Phew that was a lot about the different flags and I could talk your ear off for hours about all the options but instead just head over to the GMT wiki and check it out. There is extensive documentation that goes into great detail about all the options.
Now let's plot some coastlines and political boundaries with the command pscoast. Because we have already set up the map region and scale we can just include the -J and -R flags without including the actual values. The -W command always sets the thickness of stroke objects (i.e. lines) and the stroke color which are separated by a comma. The -D and -N flags sets the resolution of the boundaries plotted and the items plotted (lakes, rivers, ...), respectively. Again check out wiki for details. Now because we want to place this output on top of our previous output we include the -O flag (overlay) and use >> rather than > when sending the output to the postscript file. We include another -K flag as more commands are coming!
The next command is rather redundant but sets up another series of tick marks for our current map. Also this is a great command to replace our grdimage command if you'd like to make a quick map without a topography background. Sometimes the grdimage command can take quite a while to run especially if you have a high resolution grid. Thus you may want to trouble shoot your code without having to wait for the topography to plot each time.
Our final command is an example of the ubiquitous psxy command. This is what we use to plot any and all symbols from lines to stars to triangles to self defined polygons. Like most GMT commands you can get quite complex with this tool if you'd like but the example above is very straight forward. In this example I have a plain text file with the latitude and longitude of volcanoes included from the Smithsonian Global Volcanism Program. The -S flag sets the type of symbol and the size. In this case we are plotting a triangle (-St) with a diameter of 0.15 inches. We set the fill color (think paint bucket) of the triangle with the -G flag. I also include the -: flag which tells GMT that my file is written latitude - longitude rather than longitude - latitude. If you ever are trying to plot something and no symbols appear but no errors are occurring I'll bet you it's because you've flipped the lat and lon columns. Thus throw a -: flag in and see if that helps! Lastly note that here I have left out the -K flag as we expect this to be the last command sent to the postscript file.
Phew! That was a lot of information about GMT. But you know what they say - you can't run before you walk! In upcoming posts we'll get into coloring and sizing symbols by different values, making color scales, and plotting multiple maps on the same page.
Until then, get out there, download GMT, and make a map! I know you can do it :)