try-gmt-python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hello! Thank you for giving [GMT/Python](http://gmtpython.xyz/) a try.\n",
    "\n",
    "\n",
    "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 (double click on it to see the source) and it can contain text, hyperlinks, images, and even Latex equations.\n",
    "\n",
    "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.\n",
    "\n",
    "This is an example of what you can currently do with GMT/Python. Feel free to **experiment, change the code, and create new cells**. \n",
    "\n",
    "If you run into **any problems or bugs**, please [create a Github issue](https://github.com/GenericMappingTools/gmt-python/issues/new) explaining what went wrong.\n",
    "\n",
    "For more information and **if you'd like to get involved**, visit the official website: https://www.gmtpython.xyz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the library\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gmt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making a map\n",
    "\n",
    "All figure generation in GMT/Python is handled by the `gmt.Figure` class. \n",
    "It has methods to add layers to your figure, like a basemap, coastlines, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start a new figure by creating an instance of `gmt.Figure`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.coast(region=[-90, -70, 0, 20], projection='M6i', land='chocolate', \n",
    "          frame=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see a preview of the figure directly in the [Jupyter notebook](http://jupyter.org) using `fig.show()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving the figure\n",
    "\n",
    "Unlike the GMT command-line interface, **no figure file was generated until you ask for one**.  \n",
    "That means that `fig.show` won't produce a figure file.\n",
    "\n",
    "Use method `fig.savefig` (based on the [matplotlib](https://matplotlib.org/) function) to save your figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig('central-america.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample data\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quakes = gmt.datasets.load_usgs_quakes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For tabular data, the data are loaded as a `pandas.DataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quakes.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting point data\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.coast(region=\"g\", projection=\"W0/10i\", land=\"grey\", frame=True)\n",
    "fig.plot(x=quakes.longitude, y=quakes.latitude, \n",
    "         style=\"c0.1i\", color=\"yellow\", pen=\"black\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.coast(region=\"g\", projection=\"W0/10i\", land=\"grey\", frame=True)\n",
    "fig.plot(x=quakes.longitude, y=quakes.latitude, \n",
    "         sizes=0.005*2**quakes.mag, \n",
    "         style=\"ci\", color=\"yellow\", pen=\"black\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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\"`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.coast(region=\"g\", projection=\"W0/10i\", land=\"grey\", frame=True)\n",
    "fig.plot(x=quakes.longitude, y=quakes.latitude, \n",
    "         sizes=0.005*2**quakes.mag, color=quakes.depth/quakes.depth.max(),\n",
    "         style=\"ci\", pen=\"black\", cmap=\"viridis\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We still don't have support for adding a colorbar or customizing the colormaps. We're working on it.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting grids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GMT uses netCDF as its default grid format. In GMT/Python, we adopted the [xarray](http://xarray.pydata.org/) `DataArray` to represent grids. Let's load a grid of Earth relief at 30 arc-minute resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = gmt.datasets.load_earth_relief(resolution=\"30m\")\n",
    "grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regular grids can be plotted using the `Figure.grdimage` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.basemap(region=\"g\", projection=\"W0/10i\", frame=True)\n",
    "fig.grdimage(grid, cmap=\"geo\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also layer on the earthquakes on this map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.basemap(region=\"g\", projection=\"W0/10i\", frame=True)\n",
    "fig.grdimage(grid, cmap=\"geo\")\n",
    "fig.plot(x=quakes.longitude, y=quakes.latitude, \n",
    "         sizes=0.005*2**quakes.mag, color=quakes.depth/quakes.depth.max(),\n",
    "         style=\"ci\", pen=\"black\", cmap=\"viridis\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using the special file names\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.basemap(region=\"g\", projection=\"W0/10i\", frame=True)\n",
    "fig.grdimage(\"@earth_relief_30m\", cmap=\"geo\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic hillshading\n",
    "\n",
    "GMT can do automatic hillshading based on the input grid. Use `shading=True` to get the default shading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.basemap(region=\"g\", projection=\"W0/10i\", frame=True)\n",
    "fig.grdimage(\"@earth_relief_30m\", cmap=\"geo\", shading=True)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive visualization\n",
    "\n",
    "You can visualize the GMT/Python generated figure in an **interactive virtual globe** using NASA's [Web WorldWind](https://worldwind.arc.nasa.gov/web/). 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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gmt.Figure()\n",
    "fig.grdimage(\"@earth_relief_10m\", cmap=\"magma\", shading=True,\n",
    "             region=\"BR\", projection=\"X10id/10id\")\n",
    "fig.show(method=\"globe\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:try-gmt-python]",
   "language": "python",
   "name": "conda-env-try-gmt-python-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}