Arctic sea ice and snow cover and are two of the most prominent features of the cryosphere of the northern hemisphere and can be seen with a naked eye from the Moon. Measurements of sea ice area started around 1979 and snow cover a bit earlier. Both show a strong seasonal signal and area also a decline over the three decades. Even stronger is the decline calculated by a sea ice model run by the Polar Science Center, Washington. PIOMAS outputs daily sea ice volume for the same time allowing a good comparison of the three data sets. Instead of the usual line charts this notebook translates the daily data into a color range. Changes within a dataset are far better visible, while maintaining comparability.

A preview of the final results:¶

In [72]:

%pylab inline
import matplotlib.pyplot as plt
from matplotlib._png import read_png

fg, axes = plt.subplots(figsize=(20, 15), dpi=100, nrows=3, ncols=1)
for im, ax in zip(["piomas-colored.png", "area-colored.png", "snow-colored.png"], axes) :
    ax.axes.get_yaxis().set_visible(False)
    ax.axes.get_xaxis().set_visible(False)
    ax.set_frame_on(False)
    ax.imshow(read_png(im))

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['axes']
`%pylab --no-import-all` prevents importing * from pylab and numpy

Preparation, imports and helpful functions¶

In [73]:

%pylab inline

import gzip, urllib, time
from StringIO import StringIO
from datetime import datetime, timedelta

import numpy as np
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
from matplotlib._png import read_png
from matplotlib.colors import LinearSegmentedColormap
from IPython.core.display import Image

now = datetime.utcnow()

print
print "   now: %s" % now
print "python: %s" % str(sys.version_info[:])
print " numpy: %s" % numpy.__version__
print "   mpl: %s" % matplotlib.__version__

Populating the interactive namespace from numpy and matplotlib

   now: 2014-01-14 20:19:07.230332
python: (2, 7, 4, 'final', 0)
 numpy: 1.7.1
   mpl: 1.2.1

WARNING: pylab import has clobbered these variables: ['axes']
`%pylab --no-import-all` prevents importing * from pylab and numpy

In [74]:

## little helper to pixel-perfectly render a figure as an inline image
def showFigure(fg, filename, dpi=100) :
    plt.savefig(filename, dpi=dpi)
    plt.close(fg) 
    return Image(filename=filename)

None of the build-in color palettes really picks up the details of the data. Here is how to make an optimized color palette for this purpose:

In [75]:

## the final result:
Image(filename="colorbar.png")

Out[75]:

In [76]:

# human readable version
palette  = [
    (0.00, "000000"), # black
    (0.04, "880000"), # dark red
    (0.14, "ff0000"), # red
    (0.24, "ff8800"), # orange
    (0.34, "ffff00"), # yellow
    (0.44, "33ff00"), # green
    (0.57, "00ccff"), # light blue
    (0.74, "4444cc"), # dark blue
    (0.84, "aa00aa"), # violett
    (0.94, "660066"), # dark violett
    (1.00, "000000"), # black
]

# converter to mpl
def makeColormap(pal) :

	class Color(object) :
		def __init__(self, color) :
			self.r = float(int(color[0:2], 16))/255
			self.g = float(int(color[2:4], 16))/255
			self.b = float(int(color[4:6], 16))/255
		def __getitem__(self, rgb) :
			return self.__dict__[rgb]

	p = {}
	
	for col in ['red', 'green', 'blue']:
		
		akku = []; l = col[0]; 
		
		for i, e in enumerate(pal) :

			c = Color(e[1])
			if   i == 0:           akku.append((e[0],  0.0, c[l]))
			elif i == len(pal)-1 : akku.append((e[0], c[l],  1.0))
			else :                 akku.append((e[0], c[l], c[l]))
				
		p[col] = tuple(akku)
		
	return LinearSegmentedColormap('cryo', p, 255)

## convert and display
cmap = makeColormap(palette)

fg, ax = plt.subplots(figsize=(4.00, 0.5), dpi=100)
bar = mpl.colorbar.ColorbarBase(ax, cmap=cmap, orientation='horizontal')
ax.set_position((0., 0., 1., 0.99))

showFigure(fg, "colorbar.png")

Out[76]:

The plotting function for all three data sets. The only differences are axis decoration and scaling.

In [105]:

def plotData(data, vmax, xlocs, xlbls, ylocs=[], ylbls=[]) :
    
    fg, ax = plt.subplots(figsize=(5, 3), dpi=100, nrows=1, ncols=1)
    
    ax.set_position([0.08, 0.12, 0.81, 0.78])
    
    ## here the magic happens
    img = ax.imshow(data, interpolation="nearest", aspect='auto', 
        vmin=0, vmax=vmax, cmap=cmap)
    
    ## only decoration left over
    ax.set_xticks(xlocs)
    ax.set_xticklabels(xlbls)
    ax.tick_params(axis='x', labelsize=8)
    
    ylbls = ylbls if len(ylbls) else \
        [time.strftime("%b", datetime(2000, m, 1).timetuple()) for m in range(1, 13)]
    ylocs = ylocs if len(ylocs) else \
        [int(time.strftime("%j", datetime(2000, m, 1).timetuple())) for m in range(1, 13)]
    ax.set_yticks(ylocs)
    ax.set_yticklabels(ylbls)
    ax.tick_params(axis='y', labelsize=10)
    
    ax.grid(axis='y')
    
    cax = fg.add_axes([0, 0.12, 1.08, 0.78], frameon=False)
    cax.axes.get_xaxis().set_visible(False)
    cax.axes.get_yaxis().set_visible(False)
    bar = plt.colorbar(img, ax=cax, orientation='vertical')
    fts = [t.set_fontsize(8) for t in bar.ax.get_yticklabels()]
    
    return fg, ax

Volume First¶

PIOMAS volume comes as a three column text file, the import is simple. The data is usually updated once per month, well, except during a government shutdown.

In [78]:

url = "http://psc.apl.washington.edu/wordpress/wp-content/uploads/schweiger/ice_volume/PIOMAS.vol.daily.1979.2013.Current.v2.dat.gz"
rawVolu = np.genfromtxt(gzip.GzipFile( \
    fileobj=StringIO(urllib.urlopen(url).read())), skiprows=1)

# check
print rawVolu[-1], "Last day", \
    datetime(int(rawVolu[-1,0]), 1, 1) + timedelta(int(rawVolu[-1,1]) -1)
    
tmp = plot(rawVolu[:,2])

[ 2013.      365.       15.225] Last day 2013-12-31 00:00:00

I hope you'll like my tendency to spend similar variables the same amount of characters. Readability tops orthography ;)

In [79]:

# shape and resize the data, nan for missing
yirs = rawVolu[:, 0]
days = rawVolu[:, 1]
volu = rawVolu[:, 2].copy()

yirs = np.arange(int(yirs.min()), int(yirs.max() +1), 1)
days = np.arange(int(days.min()), int(days.max() +1), 1)
volu.resize((len(yirs), len(days))); 
volu[volu==0] = np.nan
volu = volu.transpose()

# PIOMAS ignores leap years, which makes days a little bit longer in a leap year.
print "shape: %s, min: %s, max: %s" % \
    (volu.shape, np.nanmin(volu), np.nanmax(volu))

shape: (365, 35), min: 3.261, max: 33.035

Now same data with volume encoded as color.

In [106]:

xlbls = np.arange(1980, 2011, 5)
xlocs = [lbl - 1979 for lbl in xlbls]

fg, ax = plotData(volu, 34, xlocs, xlbls)

plt.figtext(0.5, 0.98, u"PIOMAS - Arctic Sea Ice Volume (1000km³)",
    ha='center', va="top", fontsize=12)

plt.figtext(0.98, 0.00, "Credit: Zhang, Jinlun and D.A. Rothrock\n" + 
    "http://psc.apl.washington.edu/", ha='right', va='bottom', fontsize=5)

showFigure(fg, "piomas-colored.png")

Out[106]:

Sea Ice Area¶

Sea ice area ignores the third dimension and reduces the ice pack to a flat feature. Will this data expose a similar trend like volume?

In [86]:

url = "http://arctic.atmos.uiuc.edu/cryosphere/timeseries.anom.1979-2008"
rawArea = np.genfromtxt(url)
print rawArea[-1:]

[[  2.01403280e+03  -8.46257400e-01   1.19673138e+01   1.28135710e+01]]

Data is organized in four columns: time of year, anomaly, actual area and climatology.

In [91]:

## check
fg, axes = plt.subplots(figsize=(7, 10), dpi=100, nrows=3, ncols=1)
for dt, ax in zip([rawArea[:,1], rawArea[:,2], rawArea[:,3]], axes) :
    ax.plot(dt)

In [92]:

## shape and resize
nans = np.empty(365); nans.fill(np.nan)
yirs = range(1979, 2014)

## Only the first 365 entries of a year will be used, dropping Dec 31th in case of a leap year.
def splitYearsArea(a, y, f) :
    d = a[(a[:, 0] >= y) & (a[:, 0] < y +1)]
    return np.hstack((d[:, f], nans))[:365]

area = np.transpose(np.array([splitYearsArea(rawArea, y, 2) for y in yirs]))

print "shape: %s, min: %s, max: %s" % (area.shape, np.nanmin(area), np.nanmax(area))

shape: (365, 35), min: 2.2340095, max: 15.074604

Plotting of area is similar to thickness

In [107]:

xlbls = np.arange(1980, 2011, 5)
xlocs = [lbl - 1979 for lbl in xlbls]

fg, ax = plotData(area, 16, xlocs, xlbls)

plt.figtext(0.5, 0.98, u"Arctic Sea Ice Area (1,000,000km²)",
    ha='center', va="top", fontsize=12)

plt.figtext(0.98, 0.00, "Credit: The Cryosphere Today\n" + 
    "http://arctic.atmos.uiuc.edu/cryosphere/", \
    ha='right', va='bottom', fontsize=5)

showFigure(fg, "area-colored.png")

Out[107]:

NH Snow Cover¶

In [95]:

url = "http://climate.rutgers.edu/snowcover/files/wkcov.nhland.txt"
rawSnow = np.genfromtxt(url)
print rawSnow[-1:]

[[  2.01400000e+03   1.00000000e+00   4.57908620e+07]]

Snow comes as weekly data and includes a few gaps, so processing is a bit different.

In [96]:

url = "http://climate.rutgers.edu/snowcover/files/wkcov.nhland.txt"
rawSnow = np.genfromtxt(url)

snow = np.zeros((2015 - 1966) * 52); 
snow.fill(np.nan)
for w in rawSnow : snow[(w[0] - 1966) * 52 + w[1]] = w[2] / 1000000
snow = snow.reshape(2015 - 1966, 52).transpose()

print "shape: %s, min: %s, max: %s" % \
    (snow.shape, np.nanmin(snow), np.nanmax(snow))

## check
tmp = plot(rawSnow[:,2])

shape: (52, 49), min: 2.026924, max: 53.920493

In [108]:

xlbls = np.arange(1975, 2011, 5)
xlocs = [lbl - 1970 for lbl in xlbls]
ylbls = np.arange(4, 49, 4)
ylocs = ylbls -1

fg, ax = plotData(snow, 52, xlocs, xlbls, ylocs, ylbls)

plt.figtext(0.5, 0.98, u"NH Weekly Snow Cover (1,000,000km²)",
    ha='center', va="top", fontsize=12)

plt.figtext(0.98, 0.00, "Credit: Global Snow Lab\n" + 
    "http://climate.rutgers.edu/snowcover/index.php", \
    ha='right', va='bottom', fontsize=5)

showFigure(fg, "snow-colored.png")

Out[108]:

In [112]:

fg, axes = plt.subplots(figsize=(16, 9), dpi=72, nrows=3, ncols=1)
plt.subplots_adjust(left=0., right=1., top=1.0, bottom=0., wspace=0., hspace=0.) 

for im, ax in zip(["piomas-colored.png", "area-colored.png", "snow-colored.png"], axes) :
    ax.axes.get_yaxis().set_visible(False)
    ax.axes.get_xaxis().set_visible(False)
    ax.set_frame_on(False)
    ax.imshow(read_png(im))

There are comments.

EarthPy

Other articles

Near realtime data from Arctic ice mass balance buoys

Other articles

links

social