mapping_seismic_stations.ipynb

{
 "metadata": {
  "name": "mapping_seismic_stations"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Mapping seismic stations in the Himalayas with Numpy and Matplotlib\n",
      "## Or reading datasets with custom dtypes and plotting Earth-based data with basemap"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this exercise, we consider loading measurement files with the format:\n",
      "\n",
      "<pre>\n",
      "# Station  Lat    Long   Elev \n",
      "BIRA\t26.4840\t87.2670\t0.0120\n",
      "BUNG\t27.8771\t85.8909\t1.1910\n",
      "etc...\n",
      "</pre>\n",
      "\n",
      "These are seismic measurement stations in the Himalaya, with the elevation indicated in km.  Data with a structure such as this is common in many disciplines, and because we have a combination of text and numerical fields, we can't directly load it into a regular numpy array.\n",
      "\n",
      "But we can use numpy's ability to [define custom data types (dtypes)](http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html) to compactly describe our data in a single array, which we can then manipulate.\n",
      "\n",
      "If you have the basemap matplotlib toolkit installed, at the end of this example we will show a real Earth map and overlay the station locations on top of that.\n",
      "\n",
      "We start by configuring pylab support and loading the required modules."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we need to describe this dataset.  There are several ways of declaring a dtype, in this simple case we show two equivalent ones.  See the [numpy reference docs](http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html) for more details:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# Data descriptor to make a proper array.\n",
      "dt = [('station', (str, 4)), ('lat', np.float64), ('lon', np.float64), ('elev', np.float64) ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we load the data using this dtype we've constructed, and view it as a recarray for convenient named-field access:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "tab = np.loadtxt('data/stations.txt', dt)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ptitle = 'Seismic stations in the Himalaya'\n",
      "print ptitle\n",
      "print 'Stations:', tab['station']\n",
      "print 'Elevations (km):', tab['elev']\n",
      "print 'First station:', tab[0]\n",
      "print 'Mean latitude:', tab['lat'].mean()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f1, ax = plt.subplots(figsize = (8,5))\n",
      "\n",
      "# Make the size of the circles proportional to the elevation\n",
      "sizes = 40*(tab['elev'] + 1)\n",
      "s = ax.scatter(tab['lon'], tab['lat'], s=sizes, c=tab['elev'])\n",
      "\n",
      "# The colorbar must be associated with the return value of scatter()\n",
      "f1.colorbar(s)\n",
      "ax.set_title(ptitle)\n",
      "# Now add text labels for all the stations.  \n",
      "\n",
      "for record in tab:\n",
      "    ax.text(record['lon'] + 0.1, record['lat'] + 0.1, record['station'], weight='bold')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we find the matplotlib basemap toolkit, we can show an even better plot by\n",
      "overlaying the stations on top of a map of Earth at that location.  But we\n",
      "check this import so the code runs even without basemap."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "try:\n",
      "    from mpl_toolkits.basemap import Basemap\n",
      "except ImportError:\n",
      "    pass\n",
      "else:\n",
      "    # Draw the stations on a real map of the Earth.\n",
      "    # Find boundaries \n",
      "    lon0 = 0.995*tab.lon.min()\n",
      "    lon1 = 1.01*tab.lon.max()\n",
      "    lat0 = 0.995*tab.lat.min()\n",
      "    lat1 = 1.01*tab.lat.max()\n",
      "    # Geographic grid to draw\n",
      "    parallels = np.linspace(lat0, lat1, 5)\n",
      "    meridians = np.linspace(lon0, lon1, 5)\n",
      "\n",
      "    # Resolution of the basemap to load ('f' is *very* expensive)\n",
      "    resolution = 'i' # intermediate resolution for map info\n",
      "\n",
      "    f2, ax2 = plt.subplots(figsize=(10,6))\n",
      "    m = Basemap(lon0, lat0, lon1, lat1, resolution=resolution, ax=ax2)\n",
      "    m.drawcountries(color=(1,1,0))  # country boundaries yellow\n",
      "    m.drawrivers(color=(0,1,1))  # rivers in cyan\n",
      "    m.bluemarble()  # NASA bluemarble image\n",
      "    m.drawparallels(parallels, labels=[1,0,0,0], fmt='%.2f')\n",
      "    m.drawmeridians(meridians, labels=[0,0,0,1], fmt='%.2f')\n",
      "    s = m.scatter(tab.lon, tab.lat, s=sizes, c=tab.elev, zorder=10, alpha=0.6)\n",
      "    f2.colorbar(s)\n",
      "    ax2.set_title(ptitle)\n",
      "    for record in tab:\n",
      "        ax2.text( record['lon']+0.05, record['lat']+0.05, record['station'], \n",
      "                weight='bold', color='yellow', zorder=10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}