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))
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__
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])
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))
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:]
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))
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:]
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])
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))
Comments !