# Colormaps

Before we start, here are some useful links on some of the topics covered in this post.

Using Linear Segmented Colormaps

Scaling the indices in a cmap by a function

## The basics

The color dictionary in a colormap consists of a set of tuples for ‘red’, ‘green’, ‘blue’, each with three elements.  Each of the tuples has the form (x,y0,y1), which sets the RGB values that will be used for interpolation in the colorbar.  I find it more instructive to think about them as ( cm_index, below_interp_to, above_interp_from ):

•   ‘cm_index’ is the index in the colormap, from 0 to 1, that corresponds to the normalized range of the data values. That is, cm_index=0 in the colormap corresponds to the minimum data value plotted, and cm_index=1 in the colormap corresponds to the maximum data value that is plotted.  If you think about it like a colorbar, cm_index=0 is the bottom color, and cm_index=1 is the top color.
•   ‘below_interp_to’ is the R, G, or B value (from 0 to 1) that is used for interpolation below the ‘cm_index’ value in the tuple
•   ‘above_interp_from’ is the R, G, or B value (from 0 to 1) that is used for interpolation above the ‘cm_index’ value in the tuple

How to create a colormap:

from matplotlib import colors, cm, pyplot as plt

cdict = {
'red'  :  ((0., 0., 0.), (0.5, 0.25, 0.25), (1., 1., 1.)),
'green':  ((0., 1., 1.), (0.7, 0.0, 0.5), (1., 1., 1.)),
'blue' :  ((0., 1., 1.), (0.5, 0.0, 0.0), (1., 1., 1.))
}

#generate the colormap with 1024 interpolated values
my_cmap = colors.LinearSegmentedColormap('my_colormap', cdict,1024)



Accessing the values in an existing colormap.

cmap_old=cm.get_cmap('jet')

cmapvals=cmap_old._segmentdata
#--OR--#

print cmapvals
#  -->
#    {'blue': ((0.0, 0.5, 0.5),
#              (0.11, 1, 1),
#              (0.34, 1, 1),
#              (0.65, 0, 0),
#              (1, 0, 0)),
#     'green': ((0.0, 0, 0),
#              (0.125, 0, 0),
#              (0.375, 1, 1),
#              (0.64, 1, 1),
#              (0.91, 0, 0),
#              (1, 0, 0)),
#     'red': ((0.0, 0, 0),
#             (0.35, 0, 0),
#             (0.66, 1, 1),
#             (0.89, 1, 1),
#             (1, 0.5, 0.5))}


For examples of how this works, let’s look at some possible tuples in ‘blue’: (if this gets confusing, just jump ahead to the examples that follow)

• (0,0,0)  –> This tuple defines the beginning point, cm_index=0 (the bottom of the colorbar, for the lowest data value). Blue starts at B=0 here.
• (.25,0.1,0.1) –> This says that for the colormap value 0.25 (one-quarter of the way up the colorbar), B should be interpolated from the value in the previous tuple (B=0, from the one we talked about above) to B=0.1.  Furthermore, the interpolation value to use between cm_index=0.25 and that of the next tuple (0.5, coming up next) is also B=0.1.
• (0.5,0.8,1) –> This says that for cm_index=0.5 (the half-way point on the colorbar), the blue value should interpolate B values between the previous definition point (B=0.1 at cm_index=0.25 from before) and 0.8.  It also sets the interpolation point at 1.0 between cm_index=.5 and the next point.  The fact that we set different values 0.8 and 1.0 for bwlow_interp_to and above_interp_from means that you will see a break in the colormap at cm_index=0.5.  In other words, the colormap will not have a completely smooth color transition at the half-way point.  We’ll see an example of this later.
• (1,1,1) –> This tuple defines the end point: cm_index=1, the top of the colorbar. Interpolate between B=0.8 and B=1 for points between the last tuple and this one (cm_index=0.5 to 1).

Example: Make a simple colorbar of entirely blue values, starting from black. (Black is R=0,G=0,B=0.  White is R=1,G=1,B=1.  Pure Blue is R=0,G=0,B=1.)
To do this, we can set our ‘blue’ values to start at B=0 at ‘r’ = colorbar range = 0, then interpolate up to B=1 at ‘r’ = 1.  We don’t want red or green in this example, so we will set them both to be R=0 and G=0 at ‘r’=0 as well as at ‘r’=1.

cdict_blue = {
'blue'  : ((0,0,0),  (1,1,1)),
'green' : ((0,0,0),  (1,0,0)),
'red'   : ((0,0,0),  (1,0,0))
}

cmap_blue=colors.LinearSegmentedColormap('cmap_blue',cdict_blue,1024)


Whenever the R,G,B values all equal the same value, you have a shade of gray.  Thus, we can make a simple grayscale with:

cdict_gray = {
'blue'  : ((0,0,0),  (1,1,1)),
'green' : ((0,0,0),  (1,1,1)),
'red'   : ((0,0,0),  (1,1,1))
}

cmap_gray=colors.LinearSegmentedColormap('cmap_gray',cdict_gray,1024)


Now, let’s make a special cmap that uses entirely red values (starting with black at the bottom).  But now, let’s give it a sharp break at 80%, at which point it will transition from blue to pure white.  This is useful if you want to demonstrate that some of your data goes above a threshold (80% here), and still show the variation above that threshold in another color (blue) without simply blanking it at a uniform value.

cdict_redtoblue = {
'blue'  : ((0,0,0), (0.8,0,1), (1,1,1)),
'green' : ((0,0,0), (0.8,0,0), (1,1,1)),
'red'   : ((0,0,0), (0.8,1,0), (1,1,1))
}

cmap_redtoblue=colors.LinearSegmentedColormap('cmap_redtoblue',cdict_redtoblue,1024)


Show Basic Colormaps Plot
#Plotting the example colormaps
import numpy as np
from matplotlib import colors, cm, pyplot as plt

zvals=[np.arange(1,10,.1)]

fig=plt.figure(1)
plt.imshow(zvals,cmap=cmap_blue)
plt.title(r'cmap_blue')
ax1.axes.get_xaxis().set_visible(False); ax1.axes.get_yaxis().set_visible(False)

plt.imshow(zvals,cmap=cmap_gray)
plt.title(r'cmap_gray')
ax2.axes.get_xaxis().set_visible(False); ax2.axes.get_yaxis().set_visible(False)

plt.imshow(zvals,cmap=cmap_redtoblue)
plt.title(r'cmap_redtoblue')
ax3.axes.get_yaxis().set_visible(False)

plt.savefig('cmap_examples.png',bbox_inches='tight')
plt.show()
plt.clf()


If you want to register the new cmap for permanent use on your machine:
cm.register_cmap(name=’my_cmap_name’, cmap=cmap_you_created)

## Scaling your data (rather than scaling the colormap).

It is often easier and more desirable to scale your data before plotting.  The colorbar scale will remain linear.  Simply apply some type of transformation to your data array (such as logdata=numpy.log10(data) ), then plot as you normally would.  Here is an example.

Scaling Data Code
import numpy as np

y,x=np.mgrid[0:1,0:10:0.1]
z=x**(7)
zscale=z**(1/7.)
#For z=x^beta: beta=logz/logx.  Take the last data point for a quick & dirty estimate of beta.
#beta=np.log10(z[-1])/np.log10(x[-1])

plt.rc('text',usetex=True)

fig=plt.figure(1)
plt.imshow(z,cmap='jet')
plt.title('Original Data: z = x$^{\\beta}$')
ax1.axes.get_xaxis().set_visible(False); ax1.axes.get_yaxis().set_visible(False)

plt.imshow(zscale,cmap='jet')
plt.title('Scaled Data: $\\textrm{z}^{1/ \\beta} = {\\left( \\textrm{x}^{\\beta} \\right) }^{1/\\beta}$')
ax2.axes.get_xaxis().set_visible(False); ax2.axes.get_yaxis().set_visible(False)

plt.savefig('cmap_powerscaledata.png',bbox_inches='tight')
plt.show()
plt.clf()


## Scaling the colormap

If you are dealing with log scale, you can simply set the normalization to ‘LogNorm’ within the plot command. For instance:
pyplot.imshow(data,cmap=’jet’,norm=matplotlib.colors.LogNorm(vmin=data.min(), vmax=data.max()) )
See http://matplotlib.org/examples/pylab_examples/pcolor_log.html for an example.

If you want to scale the colormap indices to a function manually, there are a few resources that are already available online.  The scipy documentation has a helpful piece of code to accomplish this.  Below is a slightly modified version of the cmap_xmap function available at http://www.scipy.org/Cookbook/Matplotlib/ColormapTransformations

import matplotlib
def cmap_xmap(function,cmap):
""" Applies function, on the indices of colormap cmap. Beware, function
should map the [0, 1] segment to itself, or you are in for surprises.

"""
cdict = cmap._segmentdata.copy()
function_to_map = (lambda x : (function(x[0]), x[1], x[2]))
for key in ('red','green','blue'):
cdict[key] = map(function_to_map, cdict[key])
cdict[key].sort()
assert (cdict[key][0]<;0 or cdict[key][-1]>;1), "Resulting indices extend out of the [0, 1] segment."
return matplotlib.colors.LinearSegmentedColormap('colormap',cdict,1024)

#Example usage:
pscale=2.
cmap_old=cm.jet
cmap_new=cmap_xmap(lambda funcvar: funcvar**pscale,cmap_old)


If your function results in cmap indices over 1, there is an easy way to re-normalize to the 0:1 scale.  Use matplotlib.colors.BoundaryNorm() with appropriate boundaries based on the values you got from cmap_xmap

cmap_new=cmap_xmap(lambda funcvar: np.log10(funcvar),cmap_old)
bounds=[0,8,10]
scalednorm = colors.BoundaryNorm(bounds, cmap_new.N)

plt.imshow(...,cmap=...,norm=scalednorm)


Here is an example of colormaps with different power scales, making use of cmap_xmap.

z=[np.arange(1,10,.1)]

cmap_4=cmap_xmap(lambda funcvar: funcvar**4.,cm.spectral)
cmap_1_4=cmap_xmap(lambda funcvar: funcvar**(1/4.),cm.spectral)  #Note the dot after the 4!

plt.rc('font',family='serif')
fig=plt.figure(1)
plt.imshow(z,cmap='spectral')
plt.title("Original cmap 'spectral'")
ax1.axes.get_xaxis().set_visible(False); ax1.axes.get_yaxis().set_visible(False)

plt.imshow(z,cmap=cmap_4)
plt.title("cmap 'spectral' scaled to the power p=4")
ax2.axes.get_xaxis().set_visible(False); ax2.axes.get_yaxis().set_visible(False)

plt.imshow(z,cmap=cmap_1_4)
plt.title("cmap 'spectral' scaled to the power p=1/4")
ax3.axes.get_xaxis().set_visible(False); ax3.axes.get_yaxis().set_visible(False)

plt.savefig('cmap_powerscale.png',bbox_inches='tight')
plt.show()
plt.clf()


## Multiple Scaling

Sometimes normal colormap scaling tricks just don’t work.  It is difficult to highlight variations in faint features and bright features at the same time within a single data set.  In these cases, you would want the colormap to have high range at the low end and again at the high end, with little range in the middle.  I’ve created one simple solution below based on the code for cmap_xmap in the SciPy documentation.  Of course there are many ways to accomplish this (in better ways, I’m quite sure), but this is simply a tutorial.

We are scaling the colormap by two different powers (pscale1 and pscale2), and manually defining the point where we want it to make the transition (breakpoint).  We need to define a function that will scale the cmap indices by one power before the break point and then by the other power after. We must also take care this time to make sure that there is no overlap between the two powers – let’s say we decide to break at colormap index 0.5.  If we choose our scaling powers as 1/2. and 3 (–> x^(1/2.) and x^3), then colormap index of 0.4 would be scaled to 0.6324… (which is greater than the breakpoint of 0.5) while cmap index 0.7 would remap to 0.343 (below the breakpoint).  This would create a mix-up in the colorscale, and the resulting cbar would not follow the same color progression as the original.  One simple fix would be to set those values that would pass the breakpoint to be just the breakpoint.  In other words, 0.6324 –> 0.5 and 0.343 –> 0.5.  Of course this isn’t ideal since it creates a hard break, but it’s a start.

In this simplisitic solution, we make a new function called pscalefn() that takes care of all that.

#Multi-power scaler for colormaps
def cmap_powerscaler(cmap_in,breakpoint=.5,pscale1=2,pscale2=(1/2.)):
#*args are breakpoint, pscale1, pscale2
#defaults are 0.5, 2, 1/2.
try: matplotlib
except: import matplotlib

cdict = cmap_in._segmentdata.copy()
function_to_map = (lambda x : (pscalefn(x[0],breakpoint,pscale1,pscale2), x[1], x[2]))
for key in ('red','green','blue'):
cdict[key] = map(function_to_map, cdict[key])
#cdict[key].sort()
assert (cdict[key][0]<0 or cdict[key][-1]>1), "Resulting indices extend out of the [0, 1] segment."
return matplotlib.colors.LinearSegmentedColormap('colormap',cdict,1024)

def pscalefn(cval,breakpoint,pscale1,pscale2):
if cval<=breakpoint:
if cval**pscale1<=breakpoint: return cval**pscale1
else: return breakpoint
else:
if cval**pscale2<=breakpoint: return breakpoint
else: return cval**pscale2


Sample usage

z=[np.arange(1,10,.1)]

cmap_3to1_3=cmap_powerscaler(cm.spectral,breakpoint=.5,pscale1=3.,pscale2=(1/3.))
cmap_1_3to3=cmap_powerscaler(cm.spectral,breakpoint=.7,pscale1=(1/3.),pscale2=3.)

Multiple Scales Code
z=[np.arange(1,10,.1)]

cmap_3to1_3=cmap_powerscaler(cm.spectral,breakpoint=.5,pscale1=3.,pscale2=(1/3.)) #Note the dot!
cmap_1_3to3=cmap_powerscaler(cm.spectral,breakpoint=.7,pscale1=(1/3.),pscale2=3.)

plt.rc('font',family='serif')
fig=plt.figure(1)
plt.imshow(z,cmap='spectral')
plt.title("Original cmap 'spectral'",size=12)
ax1.axes.get_xaxis().set_visible(False); ax1.axes.get_yaxis().set_visible(False)

plt.imshow(z,cmap=cmap_3to1_3)
plt.title("Scaled to the power p=3 up to cm_index=0.5, and to p=1/3. after",size=12)
ax2.axes.get_xaxis().set_visible(False); ax2.axes.get_yaxis().set_visible(False)

plt.imshow(z,cmap=cmap_1_3to3)
plt.title("Scaled to the power p=1/3. up to cm_index=0.7, and to p=3 after",size=12)
ax3.axes.get_xaxis().set_visible(False); ax3.axes.get_yaxis().set_visible(False)

plt.savefig('cmap_powerscale_multiscale.png',bbox_inches='tight')
plt.show()
plt.clf()


The most useful way to see the different scales is to apply them to some real data.  Here we will again use the image of M101 as used extensively in the Kapteyn package tutorials.  http://www.astro.rug.nl/software/kapteyn/

The examples below show the following settings:

•      Original matplotlib.cm.spectral colormap with no scaling done (linear)
•      Scaled as x6, but with a break at cmap index 0.2
•      Transition from x3 to x2 with the break at cmap index 0.2
•      Transition from x(1.2) to x(1/1.5) with the break at cmap index 0.5
Galaxy Plots Code
import pyfits
import pywcsgrid2
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

m101dat=pyfits.getdata('m101.fits')

cmap2=cmap_powerscaler(cm.spectral,breakpoint=.2,pscale1=6,pscale2=6)
cmap3=cmap_powerscaler(cm.spectral,breakpoint=.2,pscale1=3,pscale2=2)
cmap4=cmap_powerscaler(cm.spectral,breakpoint=.5,pscale1=1.2,pscale2=(1/1.5)) #Note the dot!

fig=plt.figure(2)
ax1.set_ticklabel_type("delta",nbins=6)
ax1.set_xlabel("$\Delta$RA (J2000)",size=8)
ax1.set_ylabel('$\Delta$DEC (J2000)',size=8)
ax1.axis[:].major_ticklabels.set_fontsize(10)
plt.imshow(m101dat,origin='lower',interpolation='nearest',cmap='spectral')
plt.title("Original cmap 'spectral'\n(linear)",size=8)
cax1 = inset_axes(ax1,
width="5%", # width = 5% of parent_bbox width
height="100%", # height : 100%
loc=2,
bbox_to_anchor=(1.05, 0, 1, 1),
bbox_transform=ax1.transAxes,
)
cbar1=plt.colorbar(cax=cax1,orientation='vertical',cmap='spectral')
cbar1.set_label('Intensity (Some Units)',fontsize=8)
cbar1.ax.tick_params(labelsize=8)

ax2.set_ticklabel_type("delta",nbins=6)
ax2.axis[:].major_ticklabels.set_fontsize(10)
ax2.set_xlabel(""); ax2.set_ylabel("")
plt.imshow(m101dat,origin='lower',interpolation='nearest',cmap=cmap2)
plt.title("Map scaled to the power p=6,\nwith a break at 0.2",size=8)
cax2 = inset_axes(ax2,width="5%",height="100%",loc=2,bbox_to_anchor=(1.05, 0, 1, 1),bbox_transform=ax2.transAxes,borderpad=0)
cbar2=plt.colorbar(cax=cax2,orientation='vertical',cmap=cmap2)
cbar2.set_label('Intensity (Some Units)',fontsize=8)
cbar2.ax.tick_params(labelsize=8)

ax3.set_ticklabel_type("delta",nbins=6)
ax3.axis[:].major_ticklabels.set_fontsize(10)
ax3.set_xlabel(""); ax3.set_ylabel("")
plt.imshow(m101dat,origin='lower',interpolation='nearest',cmap=cmap3)
plt.title("pscale1=3, pscale2=2,\nbreakpoint=.2",size=8)
cax3 = inset_axes(ax3,width="5%",height="100%",loc=2,bbox_to_anchor=(1.05, 0, 1, 1),bbox_transform=ax3.transAxes,borderpad=0)
cbar3=plt.colorbar(cax=cax3,orientation='vertical',cmap=cmap3)
cbar3.set_label('Intensity (Some Units)',fontsize=8)
cbar3.ax.tick_params(labelsize=8)

ax4.set_ticklabel_type("delta",nbins=6)
ax4.axis[:].major_ticklabels.set_fontsize(10)
ax4.set_xlabel(""); ax4.set_ylabel("")
plt.imshow(m101dat,origin='lower',interpolation='nearest',cmap=cmap4)
plt.title("pscale1=1.2, pscale2=(1/1.5),\nbreakpoint=.5",size=8)
cax4 = inset_axes(ax4,width="5%",height="100%",loc=2,bbox_to_anchor=(1.05, 0, 1, 1),bbox_transform=ax4.transAxes,borderpad=0)
cbar4=plt.colorbar(cax=cax4,orientation='vertical',cmap=cmap4)
cbar4.set_label('Intensity (Some Units)',fontsize=8)
cbar4.ax.tick_params(labelsize=8)

plt.suptitle('M101 with Power-Scaled Colormaps',size=16,x=.5,y=.999)
plt.savefig('cmap_powerscale_images.png',dpi=200,bbox_inches='tight')
plt.show()
plt.clf()


Update: I decided to add another of my favorite colormaps.  This is similar to the STD GAMMA-II color table in IDL.  I find it really nice for astronomical images. It goes from black to blue to red to yellow to white.

cdict_coolheat={
'red'  :  ((0., 0., 0.), (0.25,0.,0.), (0.5,1.,1.), (0.75,1.0,1.0),  (1., 1., 1.)),
'green':  ((0., 0., 0.), (0.25,0.,0.), (0.5,0.,0.), (0.75,1.0,1.0),  (1., 1., 1.)),
'blue' :  ((0., 0., 0.), (0.25,1.,1.), (0.5,0.,0.), (0.75,0.0,0.0),  (1., 1., 1.))
}

coolheat = colors.LinearSegmentedColormap('coolheat', cdict_coolheat,1024)


Here it is in action with the M101 image.  On the left is the simple cmap, on the right I’ve stretched the low end to separate the background from the low-level emission.

Code for plotting M101 with coolheat colormap
import numpy as np
from matplotlib import colors, pyplot as plt
import matplotlib
import pyfits
import pywcsgrid2
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

#Defining the coolheat colormap manually
cdict_coolheat={
'red'  :  ((0., 0., 0.), (0.25,0.,0.), (0.5,1.,1.), (0.75,1.0,1.0),  (1., 1., 1.)),
'green':  ((0., 0., 0.), (0.25,0.,0.), (0.5,0.,0.), (0.75,1.0,1.0),  (1., 1., 1.)),
'blue' :  ((0., 0., 0.), (0.25,1.,1.), (0.5,0.,0.), (0.75,0.0,0.0),  (1., 1., 1.))
}

coolheat = colors.LinearSegmentedColormap('coolheat', cdict_coolheat,1024)

### Stretching the colormap manually
cdict_coolheatstretch={
'red'  :  ((0., 0., 0.), (0.12,0.0,0.0), (0.35,0.,0.), (0.55,1.,1.), (0.8,1.0,1.0),  (1., 1., 1.)),
'green':  ((0., 0., 0.), (0.12,0.0,0.0), (0.35,0.,0.), (0.55,0.,0.), (0.8,1.0,1.0),  (1., 1., 1.)),
'blue' :  ((0., 0., 0.), (0.12,0.1,0.1), (0.35,1.,1.), (0.55,0.,0.), (0.8,0.0,0.0),  (1., 1., 1.))
}
coolheat_stretch = colors.LinearSegmentedColormap('coolheatstretch', cdict_coolheatstretch,1024)

### Alternatively, we can stretch the cmap using the cmap_xmap (from
### http://www.scipy.org/Cookbook/Matplotlib/ColormapTransformations).
### Here, use power=0.8
#
# def cmap_xmap(function,cmap):
#     cdict = cmap._segmentdata.copy()
#     function_to_map = (lambda x : (function(x[0]), x[1], x[2]))
#     for key in ('red','green','blue'):
#         cdict[key] = map(function_to_map, cdict[key])
#         cdict[key].sort()
#         assert (cdict[key][0]<0 or cdict[key][-1]>1), "Resulting indices extend out of the [0, 1] segment.";
#     return matplotlib.colors.LinearSegmentedColormap('colormap',cdict,1024)
#
# coolheat_stretch=cmap_xmap(lambda funcvar: funcvar**(0.5),coolheat)

plt.rc('font',family='serif')
fig=plt.figure(3)
ax1.set_ticklabel_type("delta",nbins=6)
ax1.set_xlabel("$\Delta$RA (J2000)",size=8)
ax1.set_ylabel('$\Delta$DEC (J2000)',size=8)
ax1.axis[:].major_ticklabels.set_fontsize(10)
plt.imshow(m101dat,origin='lower',interpolation='nearest',cmap=coolheat)
plt.title('Basic',size=10)
cax1 = inset_axes(ax1,
width="5%", # width = 5% of parent_bbox width
height="100%", # height : 100%
loc=2,
bbox_to_anchor=(1.02, 0, 1, 1),
bbox_transform=ax1.transAxes,
)
cbar1=plt.colorbar(cax=cax1,orientation='vertical',cmap=coolheat)
cbar1.set_label('Intensity (Some Units)',fontsize=8)
cbar1.ax.tick_params(labelsize=8)

ax2.set_ticklabel_type('delta',nbins=6)
ax2.axis[:].major_ticklabels.set_fontsize(10)
ax2.set_xlabel(""); ax2.set_ylabel("")
plt.imshow(m101dat,origin='lower',interpolation='nearest',cmap=coolheat_stretch)
plt.title('Cmap scaled manually',size=10)
cax2 = inset_axes(ax2,width="5%",height="100%",loc=2,bbox_to_anchor=(1.02, 0, 1, 1),bbox_transform=ax2.transAxes,borderpad=0)
cbar2=plt.colorbar(cax=cax2,orientation='vertical',cmap=coolheat_stretch)
cbar2.set_label('Intensity (Some Units)',fontsize=8)
cbar2.ax.tick_params(labelsize=8)