IPython notebook used to make https://plot.ly/~etpinard/25/

polar_vortex

{
 "metadata": {
  "name": "plotly_polar_vortex"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "How cold has this winter been?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's check using IPython and plot.ly!\n",
      "\n",
      "I will plot a world map of the average daily temperature anomalies (i.e. deviation from climatology) for December 2013 and Februray 2014.\n",
      "\n",
      "This notebook will use:\n",
      "\n",
      "* plot.ly's 'heatmap' plot type for colored contour maps compatible with plot.ly [see https://plot.ly/api/python/docs/heatmaps ]\n",
      "* NetCDF data imports which are very common in Atmospheric/Oceanic sciences [ref. http://docs.scipy.org/doc/scipy/reference/generated/scipy.io.netcdf.netcdf_file.html ]\n",
      "* Basemap module for coastlines and countries data [ref. http://matplotlib.org/basemap/ ]\n",
      "\n",
      "\n",
      "Note that, I would have loved to use Basemap also to project the data to a curvilinear grid (e.g. the Robinson projection for aesthestics). <br>\n",
      "However, plot.ly's 'heatmap' plot type currently only accepts 1-dimensional 'x' and 'y' coordinate vectors. <br>\n",
      "So, my plot will use the simplest of all cylindrical projections (called the equirectangular projection) where all meridians have constant spacing.\n",
      "\n",
      "I divided the tasks as:\n",
      "\n",
      "- Setup plot.ly\n",
      "- Import the temperature data\n",
      "- Get the coastline and country boundary data\n",
      "- Setup plot layout and color map\n",
      "- Plot!"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "1) Setup plot.ly"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Import plot.ly with python api [see https://plot.ly/api/python/getting-started ], \n",
      "sign in, create shortcut to plot.ly object and check version:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import plotly\n",
      "\n",
      "username='etpinard'       # my username\n",
      "api_key = 'a35m7g6el5'    # my api key\n",
      "\n",
      "py = plotly.plotly(username,api_key)  # shortcut to plot.ly object\n",
      "plotly.__version__                    # check ploty.ly version"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "'0.5.7'"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Define a function to embed plot.ly graph in IPython notebook [ref. http://nbviewer.ipython.org/github/plotly/IPython-plotly/blob/master/Plotly%20Quickstart.ipynb ]:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import HTML\n",
      "\n",
      "def show_plot(url, width=1000, height=500):\n",
      "    s = '<iframe height=\"%s\" id=\"igraph\" scrolling=\"no\" seamless=\"seamless\" src=\"%s\" width=\"%s\"></iframe>' %\\\n",
      "    (height+50, \"/\".join(map(str,[url, width, height])), width+50)\n",
      "    return HTML(s)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Import all other packages needed for executing this notebook:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np                         # numpy ... of course\n",
      "from scipy.io import netcdf                # scipy NetCDF i/o module\n",
      "from mpl_toolkits.basemap import Basemap   # Basemap model of the matplotlib toolkit"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "2) Import temperature data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I got data on this http://www.esrl.noaa.gov/psd/data/composites/day/ .\n",
      "\n",
      "This website does not allow to *code* your output demand and/or use `wget` to download the data. <br>\n",
      "However, the data I used can be found in a only a few clicks:\n",
      "\n",
      "- Select `Air Temperature` in `Varaibles`\n",
      "- Select `Surface` in `Analysis level?`\n",
      "- Select `Dec` | `1` and `Jan` | `31` \n",
      "- Enter `2014` in the `Enter Year of last day of range` prompt\n",
      "- Select `Anomaly` in `Plot type?`\n",
      "- Select `All` in `Region of globe`\n",
      "- Click on `Create Plot` and a temporary NetCDf file can be downloaded on that new."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The data represents the Dec2013-Jan2014 average daily surface air temperature anomaly (in deg. C) with respect to 1981-2010 climatology.\n",
      "\n",
      "That is, our temperature variable of interest $T$ is\n",
      "\n",
      "$$ T = \\frac{1}{62 \\; \\text{days}} \\sum_{i=1}^{62} \\bigg( T_i - \\bar{\\hskip2pt T}_i \\bigg) $$\n",
      "\n",
      "where $\\bar{\\hskip2pt T}_i$ is daily surface air temepature 1981-2010 climatological mean\n",
      "and $i$ is the day index starting from 1 December 2013.\n",
      "\n",
      "In the NetCDF file, $T$ is refered to as `air`. <br>\n",
      "The NetCDF file also includes a longitude (`lon`), latitude (`lat`) and time singleton (`time`). <br>\n",
      "Note also that this dataset is of quite low resolution 2.5 deg x 2.5 deg (about 300 km) resolution.\n",
      "\n",
      "Now, let's extract these arrays [inspired by http://earthpy.org/interpolation_between_grids_with_basemap.html ]:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f_tas_path = 'compday.qQbv_4_mY1.nc'         # file must be in working directory\n",
      "f_tas = netcdf.netcdf_file(f_tas_path, 'r')  # will most likely have a different name if you downloaded it yourself\n",
      "\n",
      "lon = f_tas.variables['lon'][::]             # copy as an array\n",
      "lat = f_tas.variables['lat'][::-1]           # invert the latitude vector, now South to North\n",
      "air = f_tas.variables['air'][0,::-1,:]       # squeeze out the time dimension, invert latitude index\n",
      "\n",
      "# shift 'lon' from [0,360] to [-180,180] for better-looking plots\n",
      "lon_tmp = lon.copy()\n",
      "for n, l in enumerate(lon_tmp):\n",
      "    if l >= 180:\n",
      "        lon_tmp[n] = lon_tmp[n] - 360. \n",
      "lon = lon_tmp                         # gives [0,180]U[-180,-2.5]\n",
      "mid_lon = lon.shape[0]/2              # index of second of the -180 lon\n",
      "lon_east = lon[0:mid_lon]             # [0,180]\n",
      "lon_west = lon[mid_lon:]              # [-180,0]\n",
      "lon = np.hstack((lon_west, lon_east)) # stack the 2 halves\n",
      "\n",
      "# correspondingly shift the 'air' array\n",
      "air_east = air[:,0:mid_lon]\n",
      "air_west = air[:,mid_lon:]\n",
      "air = np.hstack((air_west, air_east))\n",
      "\n",
      "f_tas.close()   # close NetCDF file\n",
      "\n",
      "#print lon\n",
      "#print lat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "3) Get coastline and country boundary data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Basemap module includes data for drawing coastlines and country boundaries onto world maps. <br>\n",
      "Adding coastlines and/or country boundaries on a regular Python (i.e, non-plot.ly) figure can easily be done\n",
      "with the `.drawcoaslines()` or `.drawcountries()` Basemap methods. \n",
      "\n",
      "Here, I retrive the Basemap plotting data (or polygons) and convert them to longitude/latitude arrays compatible with plot.ly <br>\n",
      "[inspired by http://stackoverflow.com/questions/14280312/world-map-without-rivers-with-matplotlib-basemap ]\n",
      "\n",
      "In other words, the goal is to plot each *continuous* coastline and country boundraies as 1 plot.ly *trace*.\n",
      "\n",
      "Note that there are other ways to this, such as import shapefile file off the web [e.g. http://www.naturalearthdata.com/downloads/110m-cultural-vectors/ ]."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = Basemap()     # shortcut to Basemap object (no need to specify projection type)\n",
      "\n",
      "# Function assigning plot.ly plotting options for each coastline/country trace\n",
      "def MakePlotData(lon_cc,lat_cc):\n",
      "    ''' lon_cc: lon. of coastline/country, lat_cc: lat of coastline/country '''\n",
      "    return {\n",
      "            'x': lon_cc,\n",
      "            'y': lat_cc,  \n",
      "            'mode': 'lines',\n",
      "            'line': {'color': 'rgb(0,0,0)'},\n",
      "            'name': ''\n",
      "           }\n",
      "\n",
      "# Functions converting coastline/country plot polygons to lon/lat arrays\n",
      "def polygons_to_lonlat(N_poly,poly_paths):\n",
      "    ''' N_poly: number of polygon to convert, poly_paths: paths to polygons '''\n",
      "    data_tmp = []  # init. plotting list \n",
      "    for i_poly in range(N_poly):\n",
      "        poly_path = poly_paths[i_poly]\n",
      "        \n",
      "        # get the Basemap coordinates of each segment\n",
      "        coords_cc = np.array([(vertex[0],vertex[1]) \n",
      "                              for (vertex,code) in poly_path.iter_segments(simplify=False)])\n",
      "        \n",
      "        # convert coordinates to lon/lat by 'inverting' the Basemap projection\n",
      "        lon_cc, lat_cc = m(coords_cc[:,0],coords_cc[:,1], inverse=True)\n",
      "        \n",
      "        # add plot.ly plotting options\n",
      "        data_tmp.append(MakePlotData(lon_cc,lat_cc))\n",
      "        \n",
      "    return data_tmp\n",
      "        \n",
      "# 1) 'draw' coastlines using the Basemap built-in methods\n",
      "m_draw = m.drawcoastlines()\n",
      "poly_paths = m_draw.get_paths()  # get coastline polygon paths\n",
      "N_poly = 91                      # use of the 91 biggest coastlines (i.e. don't show rivers)\n",
      "data_cc1 = polygons_to_lonlat(N_poly,poly_paths)\n",
      "\n",
      "# 2) Similarly 'draw' countries boundaries\n",
      "m_draw = m.drawcountries()\n",
      "poly_paths = m_draw.get_paths() # get country boundaries polygon paths\n",
      "N_poly = len(poly_paths)     # use all countries boundaries\n",
      "data_cc2 = polygons_to_lonlat(N_poly,poly_paths)\n",
      "\n",
      "# Concatenate coastlines and country boundaries list\n",
      "data_cc = data_cc1+data_cc2\n",
      "\n",
      "#print data_cc\n",
      "#print len(data_cc), len(data_cc1), len(data_cc2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "4) Setup plot layout and color map"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Layout options:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "layout_style = {\n",
      "                'autosize': False,\n",
      "                'width': 1000,\n",
      "                'height': 500,\n",
      "                'margin': {'t':100, 'b':50, 'l':50, 'r':0, 'pad':4},\n",
      "                \"plot_bgcolor\": \"rgb(255,255,255)\",                     # white plot bg\n",
      "                \"paper_bgcolor\": \"rgb(222,222,222)\",                    # gray frame bg\n",
      "                'showlegend': False,\n",
      "                'hovermode': 'closest',\n",
      "                'title': \"Average daily surface air temperature anomalies [in deg. C] <br>\\\n",
      "                          from 2013-12-01 to 2014-01-31 with respect to 1981-2010 climatology\",   \n",
      "                'titlefont': {'color':'rgb(0,0,0)', 'size':18},\n",
      "                'xaxis': {\n",
      "                          #'title': 'longitude',\n",
      "                          'range': [lon[0],lon[-1]],       # show the whole globe\n",
      "                          'showgrid': False,\n",
      "                          'zeroline': None,\n",
      "                          'ticks': None,\n",
      "                          'showticklabels': False\n",
      "                         },\n",
      "                'yaxis': {\n",
      "                          #'title': 'latitude',\n",
      "                          'range': [lat[0],lat[-1]],        # show the whole globe\n",
      "                          'showgrid': True,\n",
      "                          'zeroline': False,\n",
      "                          'ticks': None,\n",
      "                          'showticklabels': False\n",
      "                         },\n",
      "                'annotations': [{\n",
      "                                'text': 'Data courtesy of <a href=\"http://www.esrl.noaa.gov/psd/data/composites/day/\">NOAA Earth System Research Laboratory </a>',\n",
      "                                'xref': 'paper',\n",
      "                                'yref': 'paper',\n",
      "                                'x': 0,\n",
      "                                'y': -0.1,\n",
      "                                'align': 'left',\n",
      "                                'showarrow': False \n",
      "                                }]\n",
      "                                \n",
      "                }"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By default, plot.ly linearly maps [`min(z)`, `max(z)`] to [0,1] to build its color map. <br>\n",
      "Instead, let's build a symmetric color map (and color bar) around zero, \n",
      "for a more elegant plot.\n",
      "\n",
      "So, first find the smaller possible symmetric range with integer bounds \n",
      "that includes all `z` values. <br>\n",
      "Then, using the default 'scl' map (coded below as `default_scl_map()`), \n",
      "map this range for our plot.ly call.\n",
      "\n",
      "To put emphasis on the temperature anomalies, \n",
      "let's use a variant of ColorBrewer2.0's divergent 'RdBu' color scheme <br>\n",
      "with 10 levels [ref. http://colorbrewer2.org ], with the 2 centermost levels in white."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The default 'scl' map that plot.ly utilizes\n",
      "def default_scl_map(x):\n",
      "    \"\"\" x: some vector to be map \"\"\"\n",
      "    a = 1./(air.max() - air.min())     # the slope \n",
      "    return a*(x - air.min())           # the image \n",
      "\n",
      "N_scl = 10   # number of colors \n",
      "lwb_scl = np.min([np.floor(air.min()),-np.ceil(air.max())])  # lower bound\n",
      "upb_scl = np.max([-np.floor(air.min()),np.ceil(air.max())])  # upper bound\n",
      "i_scl = default_scl_map(np.linspace(lwb_scl,upb_scl,N_scl))  # scaled levels to be used\n",
      "#print lwb_scl, upb_scl\n",
      "\n",
      "# colorbrewer2.0 RdBu \n",
      "cmap = np.array([ [103, 0, 31], [178, 24, 43], [214, 96, 77], [244, 165, 130],\n",
      "                  [253, 219, 199], [209, 229, 240], [146, 197, 222], \n",
      "                  [67, 147, 195], [33, 102, 172], [5, 48, 9] ])\n",
      "\n",
      "# invert colormap so that negative temperature anomalies are in blue  \n",
      "cmap = cmap[::-1,:]\n",
      "\n",
      "# with white at 2 centermost positions\n",
      "cmap[4,:] = [255,255,255]\n",
      "cmap[5,:] = [255,255,255]\n",
      "\n",
      "# now we are ready to make a plot.ly dictionary for our temperature variable\n",
      "data_air = {\n",
      "             'x': lon,\n",
      "             'y': lat,\n",
      "             'z': air, \n",
      "             'type': 'heatmap',\n",
      "             'scl': # scaled level, rgb('cmap') , add ',' in-between 'cmap' elements\n",
      "                    [[i_scl[i], \"rgb(\"+','.join(map(str,cmap[i,:]))+\")\"] for i in range(N_scl)], \n",
      "            }\n",
      "# 'text' option does not work when using the 'heatmap' type."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "5) Plot the results!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Concatenate the coastline/country list (above) and temperature list (below)\n",
      "data = data_cc+[data_air]\n",
      "    \n",
      "py.ioff()         # turn off new tab pop-up in browser\n",
      "    \n",
      "# call plot.ly, save fig as 'polar_vortex'\n",
      "plot_out = py.plot(\n",
      "                   data,\n",
      "                   layout=layout_style,\n",
      "                   filename='polar_vortex',\n",
      "                   fileopt='overwrite'\n",
      "                  )\n",
      "    \n",
      "# show plot in notebook using 'plot_out'\n",
      "print('\\n This plot is available at '+plot_out['url'])\n",
      "show_plot(plot_out['url'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        " This plot is available at https://plot.ly/~etpinard/25\n"
       ]
      },
      {
       "html": [
        "<iframe height=\"550\" id=\"igraph\" scrolling=\"no\" seamless=\"seamless\" src=\"https://plot.ly/~etpinard/25/1000/500\" width=\"1050\"></iframe>"
       ],
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<IPython.core.display.HTML at 0x3c0c2d0>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Conclusion:** Cold in Eastern North America, quite warm in Alaska, Greenland and most Europe.  "
     ]
    }
   ],
   "metadata": {}
  }
 ]
}