Plotting data points

This tutorial will show you how to plot data points with GMT/Python.

Start by importing the gmt Python package:

In [1]:
import gmt

Cartesian plots

The gmt.Figure class has a plot method for displaying points and lines. Let’s make a Cartesian x, y plot using some random data generated using numpy:

In [2]:
import numpy as np
In [3]:
# See the random number generator so we always
# get the same numbers
np.random.seed(42)
ndata = 100
region = [150, 240, -10, 60]
# Create some fake distribution of points and a measured value
x = np.random.uniform(region[0], region[1], ndata)
y = np.random.uniform(region[2], region[3], ndata)
magnitude = np.random.uniform(1, 5, size=ndata)

Now we can plot the data using Figure.plot and the Cartesian projection X.

In [4]:
fig = gmt.Figure()
# Create a 6x6 inch basemap using the data region
fig.basemap(region=region, projection='X6i', frame=True)
# Plot using triangles (i) of 0.3 cm
fig.plot(x, y, style='i0.3c', color='black')
fig.show()
Out[4]:
../_images/tutorials_plot-data-points_6_0.png

We can make the size of the markers follow the fake “magnitude” values by passing in the argument sizes to Figure.plot. We’ll need to scale the magnitude so that it will reflect the size in centimeters.

In [5]:
fig = gmt.Figure()
fig.basemap(region=region, projection='X6i', frame=True)
fig.plot(x, y, style='ic', color='black', sizes=magnitude/10)
fig.show()
Out[5]:
../_images/tutorials_plot-data-points_8_0.png

Plotting directly from a file

Sometimes you’ll have data in a file that you just want to plot without having to load it into Python. You can use the data argument of Figure.plot to specify the file name instead x and y. GMT will take care of loading your data and plotting.

In [6]:
# Save our fake data to a file.
np.savetxt('first-steps-data.txt',
           np.transpose([x, y, magnitude]))

The columns argument controls which columns are used as x, y, color, and size, respectively. GMT allows some basic operations on the column values before plotting. For example, adding sVALUE to a column index will multiply it by that value before plotting.

In [7]:
fig = gmt.Figure()
fig.basemap(region=region, projection='X6i', frame=True)
fig.plot(data='first-steps-data.txt', style='cc', color='red',
         columns=[0, 1, '2s0.1'])
fig.show()
Out[7]:
../_images/tutorials_plot-data-points_12_0.png

Making maps using sample data

GMT shines when plotting data on a map. We can use some sample data that is packaged with GMT to try this out. They can be accessed using special file names that begin with an @ symbol, for example @tut_quakes.ngdc. You can supply these names as the data argument in Figure.plot and other plotting functions. If you don’t have the files already, they are automatically downloaded by GMT and saved to a cache directory (usually ~/.gmt/cache).

The gmt.datasets package allows easy access to these data files as Python data types. For example, we can access the sample dataset of tsunami generating earthquakes around Japan (@tut_quakes.ngdc) as a pandas.DataFrame using the datasets.load_japan_quakes function:

In [8]:
from gmt.datasets import load_japan_quakes

quakes = load_japan_quakes()
quakes.head()
Out[8]:
year month day latitude longitude depth_km magnitude
0 1987 1 4 49.77 149.29 489 4.1
1 1987 1 9 39.90 141.68 67 6.8
2 1987 1 9 39.82 141.64 84 4.0
3 1987 1 14 42.56 142.85 102 6.5
4 1987 1 16 42.79 145.10 54 5.1

The functions returns to us the data in a pandas.Dataframe object that contains the date, hypocenter coordinates, and magnitude of the earthquakes.

Let’s make a local Mercator map of the epicenter coordinates.

In [9]:
quakes_region = [quakes.longitude.min() - 1, quakes.longitude.max() + 1,
                 quakes.latitude.min() - 1, quakes.latitude.max() + 1]
In [10]:
fig = gmt.Figure()
fig.coast(region=quakes_region, projection='M6i', frame=True,
          land='black', water='skyblue')
fig.plot(x=quakes.longitude, y=quakes.latitude,
         style='c0.3c', color='white', pen='black')
fig.show()
Out[10]:
../_images/tutorials_plot-data-points_17_0.png

As before, we can map the size of the markers to the earthquake magnitude. Because the magnitude is on a logarithmic scale, it helps to show the differences by scaling the values using a power law.

In [11]:
fig = gmt.Figure()
fig.coast(region=quakes_region, projection='M6i', frame=True,
          land='black', water='skyblue')
fig.plot(x=quakes.longitude, y=quakes.latitude,
         sizes=0.02*(2**quakes.magnitude),
         style='cc', color='white', pen='black')
fig.show()
Out[11]:
../_images/tutorials_plot-data-points_19_0.png

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 [12]:
fig = gmt.Figure()
fig.coast(region=quakes_region, projection='M6i', frame=True,
          land='black', water='skyblue')
fig.plot(x=quakes.longitude, y=quakes.latitude,
         sizes=0.02*2**quakes.magnitude,
         color=quakes.depth_km/quakes.depth_km.max(),
         cmap='viridis', style='cc', pen='black')
fig.show()
Out[12]:
../_images/tutorials_plot-data-points_21_0.png

Interactive map previews

Let’s preview this map using the interactive globe. In this case, we don’t need the frame or color in the oceans. 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 [13]:
fig = gmt.Figure()
fig.coast(region=quakes_region, projection='X6id/6id', land='gray')
fig.plot(x=quakes.longitude, y=quakes.latitude,
         sizes=0.02*2**quakes.magnitude,
         color=quakes.depth_km/quakes.depth_km.max(),
         cmap='viridis', style='cc', pen='black')
fig.show(method='globe')
Out[13]:
Your browser does not support HTML5 Canvas.