Encode a (normalized) 1D-array as ligthness component of an HLS image (or any other color system supported by colorsys).

imSpec.py

#!/usr/bin/env python -*- coding: utf-8 -*-

"""
Source of inspiration: http://astro.u-strasbg.fr/~koppen/discharge/
"""

__version__ = "Time-stamp: <2012-08-13 13:12 ycopin@lyopc469>"
__author__ = "Yannick Copin <yannick.copin@laposte.net>"

import numpy as N
import matplotlib.pyplot as P

def imSpec(y, comp="hls:1", cmap='gist_rainbow_r', width=10):
    """Encode 1D-array *y* (values in the range [0,1]) into color
    component *comp*='abc:i', where *i* is the component (0,1,2) and
    *abc* is one of the color systems supported by standard module
    'colorsys':

    * 'hls': Hue, Luminance, Saturation
    * 'hsv': Hue, Saturation, Value
    * 'rgb': Red, Green, Blue components
    * 'yiq': Luminance, Chrominance (used by composite video signals)

    Components of interest are 'hls:1' (to mimick emission lines) and
    'hsv:2' (for transmission lines).
    """

    import matplotlib as M
    import colorsys

    # Decipher component 'abc:i'
    try:
        dest,i = comp.split(':')        # abc:i
        i = int(i) % 3
        if dest == 'rgb':               # Mimick rgb_to_rgb
            convForth = lambda r,g,b: (r,g,b)
            convBack  = lambda r,g,b: (r,g,b)
        else:
            convForth = getattr(colorsys, 'rgb_to_' + dest)
            convBack  = getattr(colorsys, dest + '_to_rgb')
    except (ValueError, AttributeError,) as err:
        raise ValueError("Cannot parse component '%s': %s" % (comp, err))
    
    # Get RGBA from cmap
    rgba = M.cm.get_cmap(cmap)(N.linspace(0,1,len(y)))

    # Convert to color system ABC 
    abc = N.array([ convForth(r,g,b) for r,g,b in rgba[:,:3] ])

    # Replace appropriate component with y, and bring back to RGB
    abc[:,i] = y
    rgb = N.array([ convBack(a,b,c) for a,b,c in abc ])

    # Resize (n,3) to (width,n,3)
    im = N.resize(rgb, (width,)+rgb.shape)

    return im

if __name__=='__main__':

    x = N.linspace(0,6*N.pi,512)
    y = 1/(1 + 50*N.sin(x)**2)     # Fabry-Pérot transmission in [0,1]

    dx = x[1]-x[0]
    extent = (x[0]-dx/2, x[-1]+dx/2, -0.5, 100.5)

    def showAll(x, y, show=("hls", "hsv", "rgb", "yiq")):

        for abc in show:
            title = abc.upper()
            fig,axs = P.subplots(4,1, sharex=True, subplot_kw={'aspect':'auto'})
            axs[0].plot(x,y)
            axs[0].set(title=title, ylabel="Intensity")
            for i in range(3):
                comp = ':'.join((abc,str(i)))
                axs[i+1].imshow(imSpec(y, comp=comp, width=100), extent=extent)
                axs[i+1].set_ylabel(abc[i].upper())

    # Mimick emission and tranmission line spectra
    fig,(ax0,ax1,ax2) = P.subplots(3,1, sharex=True,
                                   subplot_kw={'aspect':'auto'})
    ax0.plot(x,  y, label='Emission')
    ax0.plot(x,1-y, label='Transmission')
    ax0.set_ylabel("Intensity")    
    ax0.legend(prop=dict(size='small'), ncol=2)

    ax1.imshow(imSpec(y, comp='hls:1', width=100), extent=extent)
    ax1.set(ylabel="Emission", yticklabels=[])

    ax2.imshow(imSpec(1-y, comp='hsv:2', width=100), extent=extent)
    ax2.set(ylabel="Transmission", yticklabels=[])

    P.show()