#!/usr/bin/env python # coding: utf-8 # # 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](http://jupyter.org). 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](http://jupyter.org) 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](http://gmt.soest.hawaii.edu/projects/gmt/wiki/Modernization) 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](https://matplotlib.org/) 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: # # ```python # 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. # # In[ ]: # In[ ]: # In[ ]: # # 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](https://agu.confex.com/agu/fm17/meetingapp.cgi/Paper/233558)". # # 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[ ]: get_ipython().run_cell_magic('bash', '', 'gmt begin Chile png\n gmt pscoast -RCL+r2 -JM15c+ -BWSne -B -Gbeige -Sblue -N1,1p\ngmt end\n') # In[ ]: Image('Chile.png', width=120) # ## Example 2 # In[ ]: get_ipython().run_cell_magic('bash', '', 'gmt begin map png\n gmt grdimage @earth_relief_05m -RMG+r2 -Cgeo -I+\n gmt coast -Wthin -BWSne -B\n gmt colorbar -DJTC -B -C+Uk\ngmt end\n') # In[ ]: Image('map.png', width=300) # ## Example 3 # In[ ]: get_ipython().run_cell_magic('bash', '', 'gmt begin islands png\n gmt set MAP_ANNOT_OBLIQUE 34 FORMAT_GEO_MAP ddd:mmF\n gmt subplot begin 2x2 -M0.05i -Fs3i/3i -BWSne -A+gwhite+p0.5p\n gmt grdimage @earth_relief_03s -R-30/30/-30/30+uk -JA159:32W/22:03N -Ba20mf5m -Csrtm -c1,1\n gmt grdimage @earth_relief_03s -R-15/15/-15/15+uk -JA109:20W/27:07S -Ba20mf5m -c1,2\n gmt grdimage @earth_relief_03s -R-30/30/-30/30+uk -JA149:22W/17:43S -Ba20mf5m -c2,1\n gmt grdimage @earth_relief_03s -R-10/10/-10/10+uk -JA138:39W/10:29S -Ba20mf5m -c2,2\n gmt subplot end\n gmt colorbar -B -Dx0/0.4i+jBC+w5i+h -Xw/2 -Yh\ngmt end\n') # In[ ]: Image('islands.png', width=500) # ## Example 4 # In[ ]: get_ipython().run_cell_magic('bash', '', 'gmt begin\n gmt figure hotspots png\n gmt grdimage @earth_relief_10m -JG200/30/6i -Cgeo -I+\n gmt coast -W -Dc -Bafg\n gmt plot @hotspots.txt -Sc0.2c -Gred\ngmt end\n') # 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[ ]: get_ipython().run_cell_magic('bash', '', '') # In[ ]: Image('XXX.png', width=500) # In[ ]: get_ipython().run_cell_magic('bash', '', '') # In[ ]: Image('XXX.png', width=500) # In[ ]: get_ipython().run_cell_magic('bash', '', '') # In[ ]: Image('XXX.png', width=500)