A modern Python interface for the Generic Mapping Tools - Demo

Thank you for giving GMT/Python and GMT6 a try!

This is a Jupyter notebook. It's an interactive computing environment where you can mix text (like this), code, and figures. The notebook is organized into cells. This is a Markdown cell (click on it to see the source) and it can contain text, hyperlinks, images, and even Latex equations.

To execute the code cells, select it and type shift+enter or click on the Run button above.

The following are examples of what you can currently do with GMT/Python. We'll also make the image for the poster background. There are some empty cells for you to experiment on the bottom of the notebook along with an example using the command line GMT6 in the new "modern mode".

Loading the library

The GMT modules are available as functions and classes in the gmt Python package. So we'll start by importing it:

In [ ]:
import gmt

A first map

All figure generation in GMT/Python is handled by the gmt.Figure class. It has methods to add layers to your figure, like a basemap, coastlines, and data.

We start a new figure by creating an instance of gmt.Figure:

In [ ]:
fig = gmt.Figure()

We add elements to the figure using its methods. For example, lets add the coastlines of Central America to a 6 inch wide map using the Mercator projection (M). Our figure will also have a nice frame with automatic ticks.

In [ ]:
fig.coast(region=[-90, -70, 0, 20], projection='M6i', land='chocolate', 
          frame=True)

You can see a preview of the figure directly in the Jupyter notebook using fig.show().

In [ ]:
fig.show()

A note for experienced GMT users

You'll probably have noticed several things that are different from classic command-line GMT. Many of these changes reflect the new GMT modern execution mode that will be part of the future 6.0 release. A few are GMT/Python exclusive (like the long argument names).

  1. The name of method is coast instead of pscoast. As a general rule, all ps* modules had their ps removed. The exceptions are: psxy == plot, psxyz == plot3d, and psscale == colorbar.
  2. The arguments don't use the GMT 1-letter syntax (R, J, B, etc). Those are still available as aliases and the methods will accept them (see below).
  3. Arguments like region can take lists instead of strings like 1/2/3/4. You can still use the string form but the list form is easier in Python.
  4. If a GMT argument has no options (like -B instead of -Baf), use a True value instead. An empty string would also be acceptable.
  5. There is no output redirecting to a PostScript file. The figure is generated in the background and will only be shown or saved when you ask for it.

We could have generated the figure above using the classic GMT argument names (but not the module names):

In [ ]:
fig_alias = gmt.Figure()
fig_alias.coast(R='-90/-70/0/20', J='M6i', G='chocolate', B=True)
fig_alias.show()

Saving the figure

Unlike the GMT command-line interface, no figure file was generated until you ask for one.

Use method fig.savefig (based on the matplotlib function) to save your figure to a file:

In [ ]:
fig.savefig('first-steps-central-america.png')

If you're running a Python script, you can tell fig.savefig to open the figure in an external viewer:

fig.savefig('first-steps-central-america.png', show=True)

Plot point data

We can use gmt.Figure.plot to plot data on our map.

First, lets create some fake data to plot using numpy:

In [ ]:
import numpy as np
In [ ]:
# See the random number generator so we always 
# get the same numbers
np.random.seed(42)

ndata = 30
region = [150, 240, -30, 60]
lon = np.random.uniform(region[0], region[1], ndata)
lat = np.random.uniform(region[2], region[3], ndata)
magnitude = np.random.uniform(1, 9, ndata)
depth = np.random.uniform(0, 1, ndata)

There are 3 ways to pass data into Figure.plot:

  1. x and y coordinates as 1d numpy arrays using the x and y arguments.
  2. A whole data table as a 2d numpy array using the data argument.
  3. A file name using the data argument.

Let's explore all of these options.

Using x and y

Now we can plot the data using Figure.plot and passing the x and y coordinate arrays:

In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='G200/30/6i', frame='ag', 
          resolution='i', area_thresh=5000, land='white', 
          water='DarkTurquoise')
# Plot using circles (c) of 0.3 cm
fig.plot(x=lon, y=lat, style='c0.3c', color='black')
fig.show()

We can make the size of the markers follow a data value by passing in the argument sizes:

In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='G200/30/6i', frame='ag', 
          resolution='i', area_thresh=5000, land='white', 
          water='DarkTurquoise')
fig.plot(x=lon, y=lat, sizes=0.1*magnitude, style='cc', color='black')
fig.show()

We can also map the colors of the markers to data by passing color as an array and providing a colormap name (cmap). We can even use the new matplotlib colormap "viridis":

In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='G200/30/6i', frame='ag', 
          resolution='i', area_thresh=5000, land='white', 
          water='DarkTurquoise')
fig.plot(x=lon, y=lat, style='cc', sizes=0.1*magnitude, 
         color=depth, cmap='viridis')
fig.show()

Using a matrix

Sometimes, you'll have data that you loaded from a file in the form of a single numpy 2d array (like a table). You can plot these data by passing in the data argument.

In [ ]:
table = np.transpose([lon, lat, magnitude, depth])

We can use the columns argument to specify the columns that correspond to x, y, color, and size, respectively. We'll use it to specify that we want to set the colors using the depths:

In [ ]:
fig = gmt.Figure()
fig.coast(region=region, projection='M6i', frame='ag', 
          resolution='i', area_thresh=5000, land='grey')
fig.plot(data=table, style='c0.8c', cmap='viridis', columns='0,1,3')
fig.show()

Using a file name

If you have data in a file and you don't need to do any calculations on it, you can pass the file name to plot directly. The syntax for plotting will be the same as the example for a data matrix.

First, we must save our sample data to a file:

In [ ]:
np.savetxt('first-steps-data.txt', table)
In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='N-165/6i', frame='ag', 
          resolution='i', area_thresh=5000, land='chocolate')
fig.plot(data='first-steps-data.txt', style='c0.4c', 
         cmap='viridis', columns='0,1,3')
fig.show()

Using sample data from GMT

GMT ships some sample data that can be accessed by passing @data_file as the data argument. For example, we can plot the earthquake epicenters from the tut_quakes.ngdc sample dataset:

In [ ]:
fig_quakes = gmt.Figure()
fig_quakes.coast(
    region=[130, 150, 35, 50], projection='M6i', 
    frame='afg', shorelines=True, land='gray', 
    water='lightblue')
fig_quakes.plot(data='@tut_quakes.ngdc', style='c0.5c', 
                color='blue', pen='faint', columns=[4, 3])
fig_quakes.show()

Making the background plot for the poster

I'll use the USGS quake data that comes with GMT. I downloaded it and saved to a numpy-friendly format in file usgs_quakes.txt.

In [ ]:
lon, lat, magnitude = np.loadtxt('usgs_quakes.txt', unpack=True)

Make a global Mercator map without any borders. Plot the quakes mapping the marker size and color to the magnitude.

In [ ]:
fig = gmt.Figure()
fig.coast(region=[-270, 90, -70, 70], projection='M10i', land='#aaaaaa', 
          water='white', resolution='l')
fig.plot(lon, lat, style='cc', sizes=0.005*2**magnitude, 
         color=magnitude/magnitude.max(), cmap='ocean')
fig.show(width=900)

If it looks OK, save it to a high resolution PNG.

In [ ]:
fig.savefig('poster_background.png', dpi=1200)

Experiment for yourself

Type anything you want in the cells below.

GMT6 demo

You can even try the GMT6 command line programs. The following are some of the plots from Paul's talk "The Generic Mapping Tools 6: Classic versus Modern Mode".

Cells that start with %%bash run their code through the bash shell (use them to execute GMT commands). We'll use the IPython.display.Image class to display in the notebook the PNG figures that you generate on the %%bash cells.

In [ ]:
from IPython.display import Image

Example 1

In [ ]:
%%bash
gmt begin Chile png
    gmt pscoast -RCL+r2 -JM15c+ -BWSne -B -Gbeige -Sblue -N1,1p
gmt end
In [ ]:
Image('Chile.png', width=120)

Example 2

In [ ]:
%%bash
gmt begin map png
    gmt grdimage @earth_relief_05m -RMG+r2 -Cgeo -I+
    gmt coast -Wthin -BWSne -B
    gmt colorbar -DJTC -B -C+Uk
gmt end
In [ ]:
Image('map.png', width=300)

Example 3

In [ ]:
%%bash
gmt begin islands png
  gmt set MAP_ANNOT_OBLIQUE 34 FORMAT_GEO_MAP ddd:mmF
  gmt subplot begin 2x2 -M0.05i -Fs3i/3i -BWSne -A+gwhite+p0.5p
    gmt grdimage @earth_relief_03s -R-30/30/-30/30+uk -JA159:32W/22:03N -Ba20mf5m -Csrtm -c1,1
    gmt grdimage @earth_relief_03s -R-15/15/-15/15+uk -JA109:20W/27:07S -Ba20mf5m -c1,2
    gmt grdimage @earth_relief_03s -R-30/30/-30/30+uk -JA149:22W/17:43S -Ba20mf5m -c2,1
    gmt grdimage @earth_relief_03s -R-10/10/-10/10+uk -JA138:39W/10:29S -Ba20mf5m -c2,2
  gmt subplot end
  gmt colorbar -B -Dx0/0.4i+jBC+w5i+h -Xw/2 -Yh
gmt end
In [ ]:
Image('islands.png', width=500)

Example 4

In [ ]:
%%bash
gmt begin
    gmt figure hotspots png
    gmt grdimage @earth_relief_10m -JG200/30/6i -Cgeo -I+
    gmt coast -W -Dc -Bafg
    gmt plot @hotspots.txt -Sc0.2c -Gred
gmt end
In [ ]:
Image('hotspots.png', width=500)

Experiment for yourself

Try running different things in the cells below or edit and rerun the cells above.

In [ ]:
%%bash
In [ ]:
Image('XXX.png', width=500)
In [ ]:
%%bash
In [ ]:
Image('XXX.png', width=500)
In [ ]:
%%bash
In [ ]:
Image('XXX.png', width=500)