How cold has this winter been?

Let's check using IPython and plot.ly!

I will plot a world map of the average daily temperature anomalies (i.e. deviation from climatology) for December 2013 and Februray 2014.

This notebook will use:

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).
However, plot.ly's 'heatmap' plot type currently only accepts 1-dimensional 'x' and 'y' coordinate vectors.
So, my plot will use the simplest of all cylindrical projections (called the equirectangular projection) where all meridians have constant spacing.

I divided the tasks as:

  • Setup plot.ly
  • Import the temperature data
  • Get the coastline and country boundary data
  • Setup plot layout and color map
  • Plot!

1) Setup plot.ly

Import plot.ly with python api [see https://plot.ly/api/python/getting-started ], sign in, create shortcut to plot.ly object and check version:

In [1]:
import plotly

username='etpinard'       # my username
api_key = 'a35m7g6el5'    # my api key

py = plotly.plotly(username,api_key)  # shortcut to plot.ly object
plotly.__version__                    # check ploty.ly version
Out[1]:
'0.5.7'

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 ]:

In [2]:
from IPython.display import HTML

def show_plot(url, width=1000, height=500):
    s = '<iframe height="%s" id="igraph" scrolling="no" seamless="seamless" src="%s" width="%s"></iframe>' %\
    (height+50, "/".join(map(str,[url, width, height])), width+50)
    return HTML(s)

Import all other packages needed for executing this notebook:

In [3]:
import numpy as np                         # numpy ... of course
from scipy.io import netcdf                # scipy NetCDF i/o module
from mpl_toolkits.basemap import Basemap   # Basemap model of the matplotlib toolkit

2) Import temperature data

I got data on this http://www.esrl.noaa.gov/psd/data/composites/day/ .

This website does not allow to code your output demand and/or use wget to download the data.
However, the data I used can be found in a only a few clicks:

  • Select Air Temperature in Varaibles
  • Select Surface in Analysis level?
  • Select Dec | 1 and Jan | 31
  • Enter 2014 in the Enter Year of last day of range prompt
  • Select Anomaly in Plot type?
  • Select All in Region of globe
  • Click on Create Plot and a temporary NetCDf file can be downloaded on that new.

The data represents the Dec2013-Jan2014 average daily surface air temperature anomaly (in deg. C) with respect to 1981-2010 climatology.

That is, our temperature variable of interest $T$ is

$$ T = \frac{1}{62 \; \text{days}} \sum_{i=1}^{62} \bigg( T_i - \bar{\hskip2pt T}_i \bigg) $$

where $\bar{\hskip2pt T}_i$ is daily surface air temepature 1981-2010 climatological mean and $i$ is the day index starting from 1 December 2013.

In the NetCDF file, $T$ is refered to as air.
The NetCDF file also includes a longitude (lon), latitude (lat) and time singleton (time).
Note also that this dataset is of quite low resolution 2.5 deg x 2.5 deg (about 300 km) resolution.

Now, let's extract these arrays [inspired by http://earthpy.org/interpolation_between_grids_with_basemap.html ]:

In [4]:
f_tas_path = 'compday.qQbv_4_mY1.nc'         # file must be in working directory
f_tas = netcdf.netcdf_file(f_tas_path, 'r')  # will most likely have a different name if you downloaded it yourself

lon = f_tas.variables['lon'][::]             # copy as an array
lat = f_tas.variables['lat'][::-1]           # invert the latitude vector, now South to North
air = f_tas.variables['air'][0,::-1,:]       # squeeze out the time dimension, invert latitude index

# shift 'lon' from [0,360] to [-180,180] for better-looking plots
lon_tmp = lon.copy()
for n, l in enumerate(lon_tmp):
    if l >= 180:
        lon_tmp[n] = lon_tmp[n] - 360. 
lon = lon_tmp                         # gives [0,180]U[-180,-2.5]
mid_lon = lon.shape[0]/2              # index of second of the -180 lon
lon_east = lon[0:mid_lon]             # [0,180]
lon_west = lon[mid_lon:]              # [-180,0]
lon = np.hstack((lon_west, lon_east)) # stack the 2 halves

# correspondingly shift the 'air' array
air_east = air[:,0:mid_lon]
air_west = air[:,mid_lon:]
air = np.hstack((air_west, air_east))

f_tas.close()   # close NetCDF file

#print lon
#print lat

3) Get coastline and country boundary data

The Basemap module includes data for drawing coastlines and country boundaries onto world maps.
Adding coastlines and/or country boundaries on a regular Python (i.e, non-plot.ly) figure can easily be done with the .drawcoaslines() or .drawcountries() Basemap methods.

Here, I retrive the Basemap plotting data (or polygons) and convert them to longitude/latitude arrays compatible with plot.ly
[inspired by http://stackoverflow.com/questions/14280312/world-map-without-rivers-with-matplotlib-basemap ]

In other words, the goal is to plot each continuous coastline and country boundraies as 1 plot.ly trace.

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/ ].

In [5]:
m = Basemap()     # shortcut to Basemap object (no need to specify projection type)

# Function assigning plot.ly plotting options for each coastline/country trace
def MakePlotData(lon_cc,lat_cc):
    ''' lon_cc: lon. of coastline/country, lat_cc: lat of coastline/country '''
    return {
            'x': lon_cc,
            'y': lat_cc,  
            'mode': 'lines',
            'line': {'color': 'rgb(0,0,0)'},
            'name': ''
           }

# Functions converting coastline/country plot polygons to lon/lat arrays
def polygons_to_lonlat(N_poly,poly_paths):
    ''' N_poly: number of polygon to convert, poly_paths: paths to polygons '''
    data_tmp = []  # init. plotting list 
    for i_poly in range(N_poly):
        poly_path = poly_paths[i_poly]
        
        # get the Basemap coordinates of each segment
        coords_cc = np.array([(vertex[0],vertex[1]) 
                              for (vertex,code) in poly_path.iter_segments(simplify=False)])
        
        # convert coordinates to lon/lat by 'inverting' the Basemap projection
        lon_cc, lat_cc = m(coords_cc[:,0],coords_cc[:,1], inverse=True)
        
        # add plot.ly plotting options
        data_tmp.append(MakePlotData(lon_cc,lat_cc))
        
    return data_tmp
        
# 1) 'draw' coastlines using the Basemap built-in methods
m_draw = m.drawcoastlines()
poly_paths = m_draw.get_paths()  # get coastline polygon paths
N_poly = 91                      # use of the 91 biggest coastlines (i.e. don't show rivers)
data_cc1 = polygons_to_lonlat(N_poly,poly_paths)

# 2) Similarly 'draw' countries boundaries
m_draw = m.drawcountries()
poly_paths = m_draw.get_paths() # get country boundaries polygon paths
N_poly = len(poly_paths)     # use all countries boundaries
data_cc2 = polygons_to_lonlat(N_poly,poly_paths)

# Concatenate coastlines and country boundaries list
data_cc = data_cc1+data_cc2

#print data_cc
#print len(data_cc), len(data_cc1), len(data_cc2)

4) Setup plot layout and color map

Layout options:

In [6]:
layout_style = {
                'autosize': False,
                'width': 1000,
                'height': 500,
                'margin': {'t':100, 'b':50, 'l':50, 'r':0, 'pad':4},
                "plot_bgcolor": "rgb(255,255,255)",                     # white plot bg
                "paper_bgcolor": "rgb(222,222,222)",                    # gray frame bg
                'showlegend': False,
                'hovermode': 'closest',
                'title': "Average daily surface air temperature anomalies [in deg. C] <br>\
                          from 2013-12-01 to 2014-01-31 with respect to 1981-2010 climatology",   
                'titlefont': {'color':'rgb(0,0,0)', 'size':18},
                'xaxis': {
                          #'title': 'longitude',
                          'range': [lon[0],lon[-1]],       # show the whole globe
                          'showgrid': False,
                          'zeroline': None,
                          'ticks': None,
                          'showticklabels': False
                         },
                'yaxis': {
                          #'title': 'latitude',
                          'range': [lat[0],lat[-1]],        # show the whole globe
                          'showgrid': True,
                          'zeroline': False,
                          'ticks': None,
                          'showticklabels': False
                         },
                'annotations': [{
                                'text': 'Data courtesy of <a href="http://www.esrl.noaa.gov/psd/data/composites/day/">NOAA Earth System Research Laboratory </a>',
                                'xref': 'paper',
                                'yref': 'paper',
                                'x': 0,
                                'y': -0.1,
                                'align': 'left',
                                'showarrow': False 
                                }]
                                
                }

By default, plot.ly linearly maps [min(z), max(z)] to [0,1] to build its color map.
Instead, let's build a symmetric color map (and color bar) around zero, for a more elegant plot.

So, first find the smaller possible symmetric range with integer bounds that includes all z values.
Then, using the default 'scl' map (coded below as default_scl_map()), map this range for our plot.ly call.

To put emphasis on the temperature anomalies, let's use a variant of ColorBrewer2.0's divergent 'RdBu' color scheme
with 10 levels [ref. http://colorbrewer2.org ], with the 2 centermost levels in white.

In [7]:
# The default 'scl' map that plot.ly utilizes
def default_scl_map(x):
    """ x: some vector to be map """
    a = 1./(air.max() - air.min())     # the slope 
    return a*(x - air.min())           # the image 

N_scl = 10   # number of colors 
lwb_scl = np.min([np.floor(air.min()),-np.ceil(air.max())])  # lower bound
upb_scl = np.max([-np.floor(air.min()),np.ceil(air.max())])  # upper bound
i_scl = default_scl_map(np.linspace(lwb_scl,upb_scl,N_scl))  # scaled levels to be used
#print lwb_scl, upb_scl

# colorbrewer2.0 RdBu 
cmap = np.array([ [103, 0, 31], [178, 24, 43], [214, 96, 77], [244, 165, 130],
                  [253, 219, 199], [209, 229, 240], [146, 197, 222], 
                  [67, 147, 195], [33, 102, 172], [5, 48, 9] ])

# invert colormap so that negative temperature anomalies are in blue  
cmap = cmap[::-1,:]

# with white at 2 centermost positions
cmap[4,:] = [255,255,255]
cmap[5,:] = [255,255,255]

# now we are ready to make a plot.ly dictionary for our temperature variable
data_air = {
             'x': lon,
             'y': lat,
             'z': air, 
             'type': 'heatmap',
             'scl': # scaled level, rgb('cmap') , add ',' in-between 'cmap' elements
                    [[i_scl[i], "rgb("+','.join(map(str,cmap[i,:]))+")"] for i in range(N_scl)], 
            }
# 'text' option does not work when using the 'heatmap' type.

5) Plot the results!

In [8]:
# Concatenate the coastline/country list (above) and temperature list (below)
data = data_cc+[data_air]
    
py.ioff()         # turn off new tab pop-up in browser
    
# call plot.ly, save fig as 'polar_vortex'
plot_out = py.plot(
                   data,
                   layout=layout_style,
                   filename='polar_vortex',
                   fileopt='overwrite'
                  )
    
# show plot in notebook using 'plot_out'
print('\n This plot is available at '+plot_out['url'])
show_plot(plot_out['url'])
 This plot is available at https://plot.ly/~etpinard/25
Out[8]:

Conclusion: Cold in Eastern North America, quite warm in Alaska, Greenland and most Europe.