Category Archives: Generic Mapping Tools (GMT)

Adding real data into GMT map, making a hillshade, and positioning

Alright, let’s get down to business plotting some real data on a map. Grab the zip file here, and uncompress it in the folder you want to build a map in. The zip contains the data to map, as well as  on of my favorite color palettes from cptcity by M. Burak Yilkilmaz (file named mby.cpt).

Here’s the end product of this tutorial, we’re going to make two maps, the first will be a simple data map, and the second a hillshaded version of the same map. We’ll also learn how to position them side by side like below.

g12

We’ll start by plotting a simple map, with just the data represented as different colors (left side of the map). The data we’re plotting are from a digital elevation model (DEM), and of Eastern North America. Making this image is simple, we simply need to call
gmt grdimage real_data.nc -R278/286/36/42 -JM4 -B2WESN -Xc -Yc -Cmby.cpt -V > flat.ps
where you already know most of the switches above. grdimage is a gmt tool that you can use to plot any data stored in the netCDF format. To color the image based on the values in the .nc file, we need to use a color palette. The color palette is called with the -C switch, followed by the name of the color palette. More on color palettes in a future tutorial though.

Now, let’s talk about hillshades. A hillshade (also known as shaded relief) is an element of a map that adds depth to topography. It is essentially simulating a light source and the shadows cast on the landscape from the hills onto the valleys. Making one in GMT is straightforward but requires a few steps. The main step of the process is done first with grdgradient. We use the following command
gmt grdgradient real_data.nc -Ghillshade-grad.nc -A345 -Ne0.6 -V
In this command, -G has a different meaning than in this tutorial; another common use of the -G switch is to give the name of the output file for the command. -A is the azimuth from which the light source should be directed, and -Ne is a normalization parameter for the output from the calculation, just use 0.6.

Then, we will ensure that the shaded relief map is reasonably balanced by normalizing it again along a Gaussian distribution
gmt grdhisteq hillshade-grad.nc -Ghillshade-hist.nc -N -V

Finally, in order to ensure the values are between -1 and 1, a final step required for creating the intensity file needed for a hillshade to the image. We do this by dividing by some value just larger than the range of values in the histogram equalized grid file. Let’s find this range by looking at the info for this file.

$ gmt grdinfo hillshade-hist.nc
hillshade-hist.nc: Title: Produced by grdhisteq
hillshade-hist.nc: Command: grdhisteq hillshade-grad.nc -Ghillshade-hist.nc -N -V
hillshade-hist.nc: Remark: Normalized directional derivative(s)
hillshade-hist.nc: Gridline node registration used [Geographic grid]
hillshade-hist.nc: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
hillshade-hist.nc: x_min: 278.004166667 x_max: 286.004166667 x_inc: 0.00833333333333 name: longitude [degrees_east] nx: 961
hillshade-hist.nc: y_min: 35.9958333333 y_max: 41.9958333333 y_inc: 0.00833333333333 name: latitude [degrees_north] ny: 721
hillshade-hist.nc: z_min: -4.67873573303 z_max: 4.67873573303 name: Elevation relative to sea level [m]
hillshade-hist.nc: scale_factor: 1 add_offset: 0
hillshade-hist.nc: format: netCDF-4 chunk_size: 138,145 shuffle: on deflation_level: 3

So we will use 5 as our divisor, in the command
gmt grdmath hillshade-hist.nc 5 DIV = hillshade-int.nc
which takes our -hist file and divides every grid cell by 5, putting it within the range of -1 to 1.

Finally, we use
gmt grdimage real_data.nc -R278/286/36/42 -JM4 -B2WESN -Xc -Yc -Ihillshade-int.nc -Cmby.cpt -V > shaded.ps
the grdimage tool again, in almost exactly the same way as before, except adding the -I switch with our normalized intensity file as the input.

Now, to put the two images side by side, we need to set different -X and -Y targets for each image. These switches set the offset of the lower-left corner of the PostScript layer being produced, from the last reference point used on the PostScript page. The a argument can be added to the switches to reset the reference point to the lower-left corner of the page. Try using -Xa1 -Yc for the flat image and then overlay the hillshade image at -Xa6.5 -Yc. Remember, you’ll need to add in -K, -O, and >> in order to produce a multi-layer map.

See if you can make it work by yourself, but if you need a hint, you can get the script I used from here.

One year of graduate school travel

The 24th of August 2015 marks my one year anniversary of beginning graduate school. It has been an awesome year, filled with new science, great friends, and lots of travel. I kept a record of where I slept every night for the past year to see just how far I made it. I visited three continents (plus Central America), six states, and 32 different cities. I spent just over 30% of the year away from Houston (112 nights). I got to see some really amazing places on the planet. It puts into perspective for me how just how great of a gig being a grad student is — get paid to do cool science.

Below, I made some maps to show where I went. The color of a dot represents when I traveled there during the year (Aug 24 2014 to Aug 24 2015) according to traditional rainbow color order. The maps were made in GMT using this script.

combine

Blow ups:
global map
US map
Gulf Coast map

Cheers to what next year brings!

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

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 > APP_con.ps

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 >> APP_con.ps
-75.379    40.607167
-95.402778 29.716944
EOF

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!

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 > output.ps
gmt commandtwo -K -O >> output.ps
gmt commandthr -O >> output.ps

In the first command, we initialize the file output.ps 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 “firstscript.sh” and then open it to edit it.

cd ~/firstscript
mkdir firstscript
cd firstscript
touch firstscript.sh

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:

#!/bin/bash
# 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 firstscript.sh. Execute the file: ./firstscript.sh.

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.

A brief introduction to the Generic Mapping Tools (GMT 5.1.x)

GMT operates on a UNIX style command line interface. If you don’t know what that means, that’s fine, I’m going to break it down for you here. Basically, it means that GMT operates by you typing a string of characters that tell the system what program to run and how exactly to run that program. So, the first thing in any command is the call to the GMT suite within your system, followed by the name of the GMT program you wish to run. This is followed by a series of switches that give the GMT program specific instructions on what to do. Each switch has a series of arguments associated with it that further detail to the program how exactly you want the program executed. Let’s take a look at the following command:

gmt pscoast -R-87/-69/31/45 -JM10c -Xc -Yc -B2 -Dl -A0/0/2 -Ggrey -W > map.ps

Wow! Looks complicated, huh? Well let’s break it down and work through it.

The first part is gmt pscoast; this is the name of the program we are calling (it is a program for accessing and plotting shoreline, territory border, and other datasets) The next part is -R-87/-69/31/45. There are several GMT switches that exist as options in many GMT programs (e.g. -R,-J,-G) as the same switch, however this is not always the case. In this instance, the switch (-R is followed by the argument (-87/-69/31/45), in this case the boundaries of the area you want to plot in the order left/right/bottom/top, in decimal degrees.

Next we see -JM10c. Okay, this is a set of one switch, and two arguments for that switch. The -J tells pscoast that the following text input is about what map projection to use when plotting the data. The M indicates to use the Mercator Projection type, and lastly the 11c indicates that the entire plot should have the longest axis measure 11 centimeters. The -Xc and -Yc switches and arguments tell the program to (c)enter the plot in the output page.

The next part, -B2 tells the pscoast program how to draw the border of the map plot. The switch is -B, and the arguments indicate to break the border every 2 decimal degrees with a mark.

-D is the switch for selecting the level resolution to use in plotting the datasets. Five levels exist: (f)ull, (h)igh, (i)ntermediate, (l)ow, and (c)rude. We will use the low resolution, because we are plotting a large area.

The -A switch indicates what level of coastlines to plot. Without going into too much detail, explained here, the selection we made tells the program to plot the coastline and all lakes.

The next switch, -G is a tricky one, because depending on the program you are using has different meanings. In this program it means the color fill of land (or “dry”) areas; we have selected grey. The other common meaning for this switch is to specify an output file, however this can be ignored for the moment. Also, we will not do so here, but another popularly included switch is to use -S to specify the fill color for “wet” areas.

-W specifies to trace the shorelines. Typically there are “pen-specifying” arguments to follow the switch, but the default will be fine for our purposes now.

Now the last part, > map.ps, simply tells pscoast where to write the output file. The name is chosen by you, but needs the file extension “.ps”.

Now look back at the command we gave, that long string of seemingly arbitrary characters should seem pretty simple now, each character with a special and important meaning.

The command we dissected produces the output below. It sure is not beautiful yet, but we can see that the command worked.

GMT_basic_map

The figure output from the command dissected in this article.

Maybe you have noticed this is a map of a section of the East Coast of the United States of America. As an America, I’ll probably draw a lot of the content I post from this area, but I promise I’ll try to keep it interesting.

Note that we have just scratched the surface of the options that the pscoast program has, but this should give you a good idea of the basic syntactical nature of GMT program commands. In the next chapter of this series, we’ll make a slightly more complicated map and continue building our Generic Mapping Tools skill set!

 

Click here for the next tutorial.

Installing the Generic Mapping Tools 5 (GMT 5.1.x) on Ubuntu Linux

The Generic Mapping Tools (GMT) are a suite of tools for making maps. The graphics generated with default settings in GMT exceed in many ways those which would come out of most users’ ArcGIS work. Admittedly, GMT does come with a steep learning curve, not excluding the installation process. For earlier versions of GMT (4.x.x) installation is more simple and direct, as you can just install from the Ubuntu Software Center (or other package managers), however, these managers typically don’t have the most recent versions — in this case GMT version 5.1.x.

Without further adieu, in the following steps, we will collect the necessary components of the GMT install, compile, and install. Most of the guide comes from the GMT website, but is supplemented with some of what I think are helpful details for new users to GMT.

READ: This guide uses the terminal and a file explorer, but the entire install can be done from the terminal; in some steps incomplete commands to do operations via terminal are listed in square brackets […] where items in carrot brackets must be replaced with system/your-install specific items <…>. All terminal commands listed without brackets are required for successful install.

1) Begin by installing several dependency packages. Running the following command will be sufficient, as any already installed packages will be skipped.
sudo apt-get install subversion ghostscript build-essential cmake libnetcdf-dev libgdal1-dev libfftw3-dev libpcre3-dev

2) Download the latest stable version of GMT using
svn checkout svn://gmtserver.soest.hawaii.edu/gmt5/trunk gmt5-dev

3) Visit the GMT download page and download the latest versions of the packages titled “gshhg-gmt-x.x.x.tar.gz” and “dcw-gmt-x.x.x.tar.gz”.
[wget http://gmt.soest.hawaii.edu/files/download?name=dcw-gmt-x.x.x.tar.gz]

4) Copy each compressed folder into the directory downloaded via subversion — this should be located at ~/gmt5-dev by default. Uncompress the folders here.
[cd ~/gmt5-dev]
[cp ./ && cp ./]
[tar -zxvf gshhg-gmt-x.x.x.tar.gz]

5) In the ~/gmt5-dev folder, enter the cmake folder, make a copy of the ConfigUserTemplate.cmake file and rename the copy to ConfigUser.cmake. Open this file in an editor.
[cp ConfigUserTemplate.cmake ./ConfigUser.cmake]
[gedit ConfigUser.cmake]

6) You need to edit the following lines of this file:
a. enable (uncomment) line 107 (set (GSHHG_ROOT…) and replace the path name with the absolute path to the gshhg-gmt-x.x.x folder.
b. enable copy in line 110
c. enable line 113 and replace the path nae with the absolute path to the dcw-gmt-x.x.x folder.
d. enable copy in line 116
Save the file and return to the terminal. NOTE: the line numbers may change with updates by the GMT devs to the .cmake file. Thanks @jerrodwalker for the most recent file lines.

7) cd into the ~/gmt5-dev folder and execute the following commands, waiting to finish each time.
mkdir build
cd build
cmake ..
make
sudo make install

That should complete your install of GMT 5.1.x! To test the install, try the following command into a new terminal window
gmt pscoast -R-130/-30/-50/50 -Jm0.025i -B30g30:.Mercator: -Di -W > mercator.ps

Check out the first article of my series on making maps with GMT here!

UPDATE: These instructions were tested on 08/25/16 on Ubuntu 16.04 LTS and remain successful when followed correctly (thanks Luis M.).