{
 "metadata": {
  "name": "",
  "signature": "sha256:0880e1798924fa2dc0e0453cfc09052aaedc699867bedb8fa22f8eeb2f23ed82"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Geoscience visualization with Matplotlib and CartoPy\n",
      "\n",
      "- Sophisticated 2D and 3D plotting library.\n",
      "- We'll give you a taste of the possibilities, but covering < 1%.\n",
      "- Mature and can produce publication quality graphics.\n",
      "- Static images, nothing interactive yet.\n",
      "- We'll focus on geoscience, but matplotlib can likely meet most of your scientific plotting requirements.\n",
      "- The best is to find an example of what you want in a [gallery](http://matplotlib.org/gallery.html) and emulate that.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Some CartoPy examples to whet your appetite\n",
      "\n",
      "Reference: http://scitools.org.uk/cartopy/docs/latest/gallery.html\n",
      "\n",
      "![Basic Map](http://scitools.org.uk/cartopy/docs/latest/_images/global_map_01_00.png)\n",
      "\n",
      "![Barbs Map](http://scitools.org.uk/cartopy/docs/latest/_images/barbs_01_00.png)\n",
      "\n",
      "![Contour Map](http://scitools.org.uk/cartopy/docs/latest/_images/waves_01_00.png)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Plotting netCDF data\n",
      "\n",
      "### Basic Plot\n",
      "\n",
      "- Plotting temperature as a function of depth as predicted by the RTOFS model"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Importing libraries we will need.\n",
      "import netCDF4\n",
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Open the netCDF file\n",
      "f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')\n",
      "\n",
      "# Getting the n-dimensional data\n",
      "tempv = f.variables['temperature']\n",
      "depth = f.variables['Depth']\n",
      "\n",
      "print(\"The temperature variable...\\n\")\n",
      "# Note the temperature variable has time, depth, y and x dimensions\n",
      "print(tempv)\n",
      "print(\"The dimensions...\\n\")\n",
      "print(tempv.dimensions)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Continued from previous cell..\n",
      "\n",
      "# Our goal is temperature as a function of depth so slicing along the depth axis\n",
      "# at a specific time and place\n",
      "temp = tempv[0, :, 123, 486]\n",
      "\n",
      "# Masked arrays are arrays that have bad values idenitifed by the mask array.\n",
      "print(\"The masked array containing the temperature data...\")\n",
      "print(repr(temp))\n",
      "\n",
      "# trick for filtering out good values\n",
      "x = temp[~temp.mask] \n",
      "y = depth[~temp.mask]\n",
      "\n",
      "print(\"Plotting...\")\n",
      "# plot and show data\n",
      "plt.plot(y, x)\n",
      "\n",
      "# close netCDF\n",
      "f.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Let's build upon the previous plot into something that is ready for publication.\n",
      "\n",
      "- Title\n",
      "- Axis Legends\n",
      "- Markers\n",
      "\n",
      "[Peruse matplotlib gallery](http://matplotlib.org/gallery.html) and see and \n",
      "emulate what you like."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Adding text, adjusting borders, and figure size\n",
      "# Get figure hook to manipulate our plot\n",
      "fig = plt.figure()\n",
      "desc ='Figure 1. Temperature as a function of ocean depth as\\n'\\\n",
      "    'predicted by the RTOFS model'\n",
      "# Adding our description\n",
      "plt.figtext(.5,.15,desc,fontsize=12,ha='center')\n",
      "#adjust margin\n",
      "fig.subplots_adjust(bottom=0.35)\n",
      "#adjust figure size\n",
      "fig.set_size_inches(7,7)\n",
      "\n",
      "# Improve axes\n",
      "# Get axis hook to manipulate our plot\n",
      "ax = fig.add_subplot(111)\n",
      "# Add axis labels\n",
      "ax.set_xlabel('Depth (m)', fontweight='bold')\n",
      "# \\u00b0 : degree symbol\n",
      "ax.set_ylabel(u'Temperature (\\u00b0C)', fontweight='bold')\n",
      "# Don't show top and right axis\n",
      "ax.spines[\"right\"].set_visible(False)\n",
      "ax.spines[\"top\"].set_visible(False)\n",
      "# Define ticks\n",
      "ax.tick_params(axis='both', direction='out')\n",
      "ax.get_xaxis().tick_bottom() # remove unneeded ticks \n",
      "ax.get_yaxis().tick_left()\n",
      "\n",
      "# Getting the data as we did before\n",
      "f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')\n",
      "tempv = f.variables['temperature']\n",
      "depth = f.variables['Depth']\n",
      "temp = tempv[0,:,123,486]\n",
      "x = temp[~temp.mask] #trick for getting data\n",
      "y = depth[~temp.mask]\n",
      "\n",
      "# Plotting line with triangle markers, and red line color.\n",
      "plt.plot(y, x, marker=r'^', color='r', markersize=10, clip_on=False, linewidth=2.0)\n",
      "plt.show()\n",
      "f.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## *Exercise*\n",
      "\n",
      "- Create a new notebook cell here.\n",
      "- Based on cell above, plot salinity as a function of depth.\n",
      "- Try using a different marker.\n",
      "- Try using a different line color.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "CartoPy"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- High level API for dealing with maps\n",
      "- CartoPy allows you to plot data on a 2D map.\n",
      "- Support many different map projections\n",
      "- Support for shapefiles from the GIS world"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Importing CartoPy\n",
      "import cartopy"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Works with matplotlib's built-in transform support.\n",
      "fig = plt.figure(figsize=(10, 8))\n",
      "ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Mercator())\n",
      "\n",
      "# Sets the extent to cover the whole globe\n",
      "ax.set_global()\n",
      "\n",
      "# Adds standard background map\n",
      "ax.stock_img()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Cartopy also has a lot of built-in support for a variety of map features:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(10, 8))\n",
      "ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())\n",
      "\n",
      "# Add variety of features\n",
      "ax.add_feature(cartopy.feature.LAND)\n",
      "ax.add_feature(cartopy.feature.OCEAN)\n",
      "ax.add_feature(cartopy.feature.COASTLINE)\n",
      "\n",
      "# Can also supply matplotlib kwargs\n",
      "ax.add_feature(cartopy.feature.BORDERS, linestyle=':')\n",
      "ax.add_feature(cartopy.feature.LAKES, alpha=0.5)\n",
      "ax.add_feature(cartopy.feature.RIVERS)\n",
      "\n",
      "# Set limits in lat/lon space\n",
      "ax.set_extent([-140, -60, 25, 50])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can also grab other features from the Natural Earth project: http://www.naturalearthdata.com/"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(10, 8))\n",
      "ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())\n",
      "\n",
      "ax.add_feature(cartopy.feature.LAND)\n",
      "ax.add_feature(cartopy.feature.OCEAN)\n",
      "ax.add_feature(cartopy.feature.COASTLINE)\n",
      "ax.add_feature(cartopy.feature.BORDERS, linestyle=':')\n",
      "\n",
      "# Grab state borders\n",
      "state_borders = cartopy.feature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m',\n",
      "                                                    facecolor='none') \n",
      "\n",
      "ax.add_feature(state_borders, linestyle=\"--\", edgecolor='blue')\n",
      "\n",
      "# Set limits in lat/lon space\n",
      "ax.set_extent([-140, -60, 25, 50])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plotting data works with the set projection by specifying the transform for the data when plotting. CartoPy will take care of any needed interpolation and transformation. The Plate Carree projection (Equi-rectangular projection) specifies the data in lat/lon, but not on the earth sphere. The Geodetic projection uses lat/lon, but are considered on the sphere. This example shows the difference."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(12,8))\n",
      "ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Robinson())\n",
      "ax.set_global()\n",
      "ax.stock_img()\n",
      "ax.coastlines()\n",
      "\n",
      "# Plot circle with and without specifying the transform\n",
      "ax.plot(-105.295175, 40.013176, 'o')\n",
      "ax.plot(-105.295175, 40.013176, 'ro', transform=cartopy.crs.Geodetic())\n",
      "\n",
      "# Plot line between two lat/lon points. One using equi-rectangular and one geodetic\n",
      "plt.plot([-0.08, 132.], [51.53, 43.17], 'r', transform=cartopy.crs.PlateCarree())\n",
      "plt.plot([-0.08, 132.], [51.53, 43.17], 'g', transform=cartopy.crs.Geodetic())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Exercise"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Try it yourself: Let's look at those curved flight paths on a map\n",
      "- Pick any projection\n",
      "- Plot a map of the United States with borders, coastlines, and state borders.\n",
      "- Plot lines between Seattle, WA and Orlando, FL\n",
      "    - Plot one straight line on the map (assuming lat/lon are rectangular coordinates)\n",
      "    - Plot same points assuming they come from a sphere"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Let's plot some data from RTOFS model on a map\n",
      "\n",
      "- Data are coming from the local file system"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Open and read netCDF variables\n",
      "nc = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')\n",
      "tempv = nc.variables['temperature']\n",
      "lat = nc.variables['Latitude'][:]\n",
      "lon = nc.variables['Longitude'][:]\n",
      "data = tempv[0, 0, :, :]\n",
      "\n",
      "# Set up a stereographic projection\n",
      "proj = cartopy.crs.Stereographic(central_latitude=60., central_longitude=-120.)\n",
      "\n",
      "# Construct figure\n",
      "fig = plt.figure(figsize=(10,8))\n",
      "ax = fig.add_subplot(1, 1, 1, projection=proj)\n",
      "\n",
      "# Setting the plot size and text\n",
      "plt.figtext(.5, .15, 'Figure 1. Sea surface temperature as predicted by the RTOFS model', fontsize=12, ha='center')\n",
      "\n",
      "# define color map\n",
      "cmap = plt.cm.hsv\n",
      "\n",
      "# Nice high-level, human-readable abstractions for dealing with maps.\n",
      "ax.add_feature(cartopy.feature.LAND)\n",
      "ax.add_feature(cartopy.feature.OCEAN)\n",
      "ax.coastlines(zorder=3)\n",
      "ax.gridlines()\n",
      "\n",
      "#Color-filled contour plot\n",
      "cs = ax.contourf(lon, lat, data, 50, cmap=cmap, transform=cartopy.crs.PlateCarree(), zorder=2)\n",
      "\n",
      "# Color bar\n",
      "cbar = plt.colorbar(cs)\n",
      "cbar.set_label('$^{o}C$')\n",
      "\n",
      "nc.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## *Exercise*\n",
      "\n",
      "- Create a new notebook cell here.\n",
      "- Based on cell above, create a plot showing salinity on map.\n",
      "- Try a different [projection](http://matplotlib.org/basemap/users/mapsetup.html).\n",
      "- Try a different [color map](http://matplotlib.org/api/pyplot_summary.html#matplotlib.pyplot.colormaps) [(examples)](http://matplotlib.org/examples/pylab_examples/show_colormaps.html)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Accessing and plotting data from the TDS"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "SST and ICE"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from datetime import datetime\n",
      "date = datetime(2014, 9, 15, 0) # date to plot.\n",
      "\n",
      "# open dataset.\n",
      "dataset = netCDF4.Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')\n",
      "timevar = dataset.variables['time']\n",
      "timeindex = netCDF4.date2index(date, timevar) # find time index for desired date.\n",
      "\n",
      "# read sst.  Will automatically create a masked array using\n",
      "# missing_value variable attribute. 'squeeze out' singleton dimensions.\n",
      "sst = dataset.variables['sst'][timeindex, :].squeeze()\n",
      "\n",
      "# read lats and lons (representing centers of grid boxes).\n",
      "lats = dataset.variables['lat'][:]\n",
      "lons = dataset.variables['lon'][:]\n",
      "\n",
      "# create figure, axes instances.\n",
      "fig = plt.figure(figsize=(20, 10))\n",
      "proj = cartopy.crs.Robinson(central_longitude=-105.)\n",
      "ax = fig.add_axes([0.05, 0.05, 0.9, 0.9], projection=proj)\n",
      "\n",
      "# plot sst, then ice with imshow--this relies on evenly spaced points\n",
      "bounds = (lons[0], lons[-1], lats[0], lats[-1])\n",
      "im1 = ax.imshow(sst, extent=bounds, interpolation='nearest',\n",
      "                cmap=plt.cm.jet, transform=cartopy.crs.PlateCarree())\n",
      "\n",
      "# Add gridlines and coastlines\n",
      "ax.gridlines()\n",
      "ax.coastlines('50m')\n",
      "\n",
      "# add colorbar\n",
      "plt.colorbar(im1)\n",
      "\n",
      "# add a title.\n",
      "ax.set_title('Analysis for %s' % date)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Exercise"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Modify the last example to add the sea ice analysis to the plot\n",
      "- Try a different projection"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Winds and Pressure"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# specify date to plot.\n",
      "yyyy=1993; mm=3; dd=14; hh=0\n",
      "\n",
      "# set OpenDAP server URL.\n",
      "URLbase = \"http://nomads.ncdc.noaa.gov/thredds/dodsC/modeldata/cmd_pgbh/\"\n",
      "URL = URLbase + \"{year}/{year}{month:02}/{year}{month:02}{day:02}/pgbh00.gdas.{year}{month:02}{day:02}{hour:02}.grb2\".format(\n",
      "    year=yyyy, month=mm, day=dd, hour=hh)\n",
      "data = netCDF4.Dataset(URL)\n",
      "\n",
      "# read lats,lons\n",
      "latitudes = data.variables['lat'][:]\n",
      "longitudes = data.variables['lon'][:]\n",
      "\n",
      "# get sea level pressure and 10-m wind data.\n",
      "# mult slp by 0.01 to put in units of hPa.\n",
      "slpin = 0.01 * data.variables['Pressure_msl'][:].squeeze()\n",
      "uin = data.variables['U-component_of_wind_height_above_ground'][:].squeeze()\n",
      "vin = data.variables['V-component_of_wind_height_above_ground'][:].squeeze()\n",
      "\n",
      "# make 2-d grid of lons, lats\n",
      "lons, lats = np.meshgrid(longitudes, latitudes)\n",
      "\n",
      "# make orthographic basemap.\n",
      "proj = cartopy.crs.Orthographic(central_longitude=-60., central_latitude=60.)\n",
      "\n",
      "# create figure, add axes\n",
      "fig = plt.figure(figsize=(20, 10))\n",
      "ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=proj)\n",
      "ax.set_global()\n",
      "\n",
      "# set desired contour levels.\n",
      "clevs = np.arange(960., 1061., 8.)\n",
      "\n",
      "# plot SLP contours.\n",
      "ax.contour(lons, lats, slpin, clevs, linewidths=0.5, colors='k', transform=cartopy.crs.PlateCarree())\n",
      "ax.contourf(lons, lats, slpin, clevs, cmap=plt.cm.RdBu_r, transform=cartopy.crs.PlateCarree())\n",
      "\n",
      "stride = 15\n",
      "ax.quiver(lons[::stride, ::stride], lats[::stride, ::stride],\n",
      "          uin[::stride, ::stride], vin[::stride, ::stride],\n",
      "          transform=cartopy.crs.PlateCarree())\n",
      "\n",
      "# draw coastlines and gridlines\n",
      "ax.coastlines(linewidth=1.5)\n",
      "ax.gridlines()\n",
      "\n",
      "date = datetime(yyyy, mm, dd, hh)\n",
      "ax.set_title('SLP and Wind Vectors ' + str(date))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}