Hello! Thank you for giving GMT/Python 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 (double click on it to see the source) and it can contain text, hyperlinks, images, and even Latex equations.

To execute any cell, click on it and type Shift + Enter or click on the "Run" button in the menu bar. Executing a Markdown cell will render its content. Code execution can happen non-linearly, so you can change and rerun a cell without running all of the ones that came before it. But you'll still need to run cells that define a variable/import a module before you can use the variable/module in another cell. You can restart and clear the notebook at any time using the options in the "Kernel" menu.

This is an example of what you can currently do with GMT/Python. Feel free to experiment, change the code, and create new cells.

If you run into any problems or bugs, please create a Github issue explaining what went wrong.

For more information and if you'd like to get involved, visit the official website: https://www.gmtpython.xyz

Loading the library

You can load GMT/Python by importing the gmt Python package. Most GMT processing modules will be avialable as functions in this package. The plotting modules are methods of the gmt.Figure class.

In [ ]:
import gmt

Making a 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, etc.

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()

Saving the figure

Unlike the GMT command-line interface, no figure file was generated until you ask for one.
That means that fig.show won't produce a figure file.

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

In [ ]:
fig.savefig('central-america.png')
In [ ]:
!ls

Sample data

GMT includes sample data that are downloaded automatically to a custom folder (usually ~/.gmt/cache). You can load these datasets into Python using the functions in the gmt.datasets package.

In [ ]:
quakes = gmt.datasets.load_usgs_quakes()

For tabular data, the data are loaded as a pandas.DataFrame.

In [ ]:
quakes.head()

Plotting point data

Use the gmt.Figure.plot method to add points and lines to your map. By default, it will use the same projection that you used previously to setup your map. Let's setup a map using a Mollweide projection (W) to plot our earthquake locations. The point style will be 0.1 inch circles (c0.1i) colored yellow with black outlines.

In [ ]:
fig = gmt.Figure()
fig.coast(region="g", projection="W0/10i", land="grey", frame=True)
fig.plot(x=quakes.longitude, y=quakes.latitude, 
         style="c0.1i", color="yellow", pen="black")
fig.show()

We can make the size of the markers follow the earthquake magnitude by passing in the argument sizes to Figure.plot. We'll need to scale the magnitude so that it will reflect the size in inches.

In [ ]:
fig = gmt.Figure()
fig.coast(region="g", projection="W0/10i", land="grey", frame=True)
fig.plot(x=quakes.longitude, y=quakes.latitude, 
         sizes=0.005*2**quakes.mag, 
         style="ci", color="yellow", pen="black")
fig.show()

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

In [ ]:
fig = gmt.Figure()
fig.coast(region="g", projection="W0/10i", land="grey", frame=True)
fig.plot(x=quakes.longitude, y=quakes.latitude, 
         sizes=0.005*2**quakes.mag, color=quakes.depth/quakes.depth.max(),
         style="ci", pen="black", cmap="viridis")
fig.show()

We still don't have support for adding a colorbar or customizing the colormaps. We're working on it.

Plotting grids

GMT uses netCDF as its default grid format. In GMT/Python, we adopted the xarray DataArray to represent grids. Let's load a grid of Earth relief at 30 arc-minute resolution.

In [ ]:
grid = gmt.datasets.load_earth_relief(resolution="30m")
grid

Regular grids can be plotted using the Figure.grdimage method.

In [ ]:
fig = gmt.Figure()
fig.basemap(region="g", projection="W0/10i", frame=True)
fig.grdimage(grid, cmap="geo")
fig.show()

We can also layer on the earthquakes on this map.

In [ ]:
fig = gmt.Figure()
fig.basemap(region="g", projection="W0/10i", frame=True)
fig.grdimage(grid, cmap="geo")
fig.plot(x=quakes.longitude, y=quakes.latitude, 
         sizes=0.005*2**quakes.mag, color=quakes.depth/quakes.depth.max(),
         style="ci", pen="black", cmap="viridis")
fig.show()

Using the special file names

The Earth relief data can be accessed using GMT's special file names: @earth_relief_XX. We can give this file name to grdimage instad of a grid.

In [ ]:
fig = gmt.Figure()
fig.basemap(region="g", projection="W0/10i", frame=True)
fig.grdimage("@earth_relief_30m", cmap="geo")
fig.show()

Automatic hillshading

GMT can do automatic hillshading based on the input grid. Use shading=True to get the default shading.

In [ ]:
fig = gmt.Figure()
fig.basemap(region="g", projection="W0/10i", frame=True)
fig.grdimage("@earth_relief_30m", cmap="geo", shading=True)
fig.show()

Interactive visualization

You can visualize the GMT/Python generated figure in an interactive virtual globe using NASA's Web WorldWind. In this case, we don’t need the frame or color in the continents. We must also use a Cartesian projection (X) and degrees (d) for plot units so that the figure can be aligned with the globe.

In [ ]:
fig = gmt.Figure()
fig.grdimage("@earth_relief_10m", cmap="magma", shading=True,
             region="BR", projection="X10id/10id")
fig.show(method="globe")