Skip to content

Instantly share code, notes, and snippets.

@etpinard
Created Feb 28, 2014
Embed
What would you like to do?
IPython notebook used to make https://plot.ly/~etpinard/25/
{
"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": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment