In one of the three chapter I submitted for this book, Prototype colourmaps for fault interpretation, I talk about building a widget for interactive generation of grayscale colourmaps with sigmoid lightness. The whole process is well described in the chapter and showcased in the accompanying GitHub repo.

This post is an opportunity for me to reflect on how revisiting old projects is very important. Indeed, I consider it an essential part of how I approach scientific computing, and a practical way to incorporate new insights, and changes (hopefully betterments) in your coding abilities.

In the fist version of the Jupyter notebook, submitted in 2017, all the calculations and all the plotting commands where packed inside a single monster function that was passed to ipywidgets.interact. Quite frankly, as time passed this approach seemed less and less Phytonic (aka mature) and no longer representative of my programming skills and style, and increased understanding of widget objects.

After a significant hiatus (2 years) I restructured the whole project in several ways:
– Converted Python 2 code to Python 3
– Created separate helper functions for each calculation and moved them to the top to improve on both clarity and reusability of the code.
– Improved and standardized function docstrings
– Optimized and reduced the number of parameters
– Switched from interact to interactive to enable access to the colormaparray in later cells (for printing, further plotting, and exporting).

In the first post of this series I argued that we should not build colormaps for azimuth (or phase) data by interpolating linearly between fully saturated hues in RGB or HSL space.

A first step towards the ideal colormap for azimuth data would be to interpolate between isoluminant colours instead. Kindlmann et al. (2002) published isoluminant RGB values for red, yellow, green, cyan, blue, and magenta based on a user study. The code in the next block show how to interpolate between those published colours to get 256-sample R, G, and B arrays (with magenta repeated at both ends), which can then be combined into a isoluminant colormap for azimuth data.

This is a good example in general of how to interpolate to a finer sampling one or more sequence of values using the interp function from the Numpy library. In line 04 we define 7 evenly spaced points between 0 and 255; this will be the sample coordinate for the r, g, and b colours created in lines 01-03. In line 05 we create the new coordinates we will be interpolating r, g, and b values at in lines 06-08 (all integers between 0 and 256). The full code will come in the Notebook accompanying the last post in this series.

This new colormap is used in the bottom map of the figure below, whereas in the top map we used a conventional HSV azimuth colormap (both maps show the dip azimuth calculated on the Penobscot horizon). The differences are subtle, but with the isoluminant colormap we are guaranteed there are no perceptual artifacts due to the random variations in lightness of the fully saturated HSV colors.

Another possible strategy to create a perceptual colormap for azimuth data would be to set lightness and chroma to constant values in LCH space and interpolate between hues. This is the Matlab colormap I previously created, and shown in Figure 4 of New Matlab isoluminant colormap for azimuth data. In the next post, I will show how to do this in Python.

In New Matlab isoluminant colormap for azimuth data I showcased a Matlab colormap that I believe is perceptually superior to the conventional, HSV-based colormaps for azimuth data, in that it does not superimposes on the data the color artifacts that plague all rainbows. However, it still has a limitation, which is that the main colours do not correspond exactly to the four compass directions N, E, W, and S.

My intention with this series is to go back to square one, deconstruct the conventional colormaps for azimuth, and build a new one that has all the desired properties of both perceptual linearity, and correct location of the main colors. All reproducible in Python.

If we wanted to build from scratch a colormap for azimuth (or phase) data the main tasks would be to generate a sequence of distinguishable colours at opposite quadrants, or compass directions (like 0 and 180 degrees, or N and S), and to wrap around the sequence with the same colour at the two ends.

But to do that, we should avoid interpolating linearly between fully saturated hues in RGB or HSL space.

To illustrate why, it is useful to look at the figure below. On the left is a hue circle with primary, secondary, and tertiary colours in a counter-clockwise sequence: red, rose, magenta, violet, blue, azure, cyan, aquamarine, electric green, chartreuse, yellow, and orange. The colour chips are placed at evenly spaced angular distances according to their hue (in radians).

Left, primary, secondary, and tertiary colour chips arranged using hue for angular distance; right, the same colour chips arranged using intensity for angular distance.

This looks familiar and seems like a natural ordering of colors, so we may be tempted in building a colormap, to just take that sequence, wrap it around at the red (or the magenta) and linearly interpolate to 256 colours to get a continuous colormap [1], and use it for azimuth data, which is how usually the conventional azimuth colormaps are built.

On the right side in the figure the chips have been rearranged according to their intensity on a counter-clockwise sequence from 0 to 255 with 0 at three hours; so, for example blue, which is the darkest colour with an intensity of 29, is close to the beginning of the sequence, and yellow, the brightest with an intensity of 225, is close to the end. Notice that the chips are no longer equidistant.

The most striking is that the blue and the yellow chips are more separated than the other chips, and for this reason blue and yellow features seem to stand out a lot more in a map when using this color sequence, which can be both distracting and confusing. A good example is Figure 3 in New Matlab isoluminant colormap for azimuth data.

Also, yellow and red, being two chips apart in the left circle in the figure above, are used to colour azimuths 60 degrees apart, and so do cyan and green. However, if we look at the right circle, we realize that the yellow and red chips are much further apart than the cyan and green chips [2] in the perceptual dimension of intensity; therefore, features colored in yellow and red could be perceived as much further apart (in azimuth) than cyan and green.

These differences may be subtle, but in my opinion they become important when dip azimuth is combined with other attributes, perhaps using a 3D colormap, and the resulting map is used for detailed structural interpretation. There is a really good example of this type of 3D colormap in Chopra and Marfurt (2007), where dip azimuth is rendered with hue modulation, dip magnitude with saturation modulation, and coherence with lightness modulation.

A code snippet with the main Python commands to generate the two polar scatterplots in the figure is listed, and explained below. The full code can be found in this Jupiter Notebook.

01 import matplotlib.colors as clr
02 keys=['red', '#FF007F', 'magenta', '#7F00FF', 'blue', '#0080FF','cyan', '#00FF80',
'#00FF00', '#7FFF00', 'yellow', '#FF7F00']
03 my_cmap = clr.ListedColormap(keys)
04 x = np.arange(12)
05 color = my_cmap(x)
06 n = 12
07 theta = 2*np.pi*(np.linspace(0,1,13))
08 r = np.ones(13)*2.5
09 area = 200*r**2 # size of color chips
10 c = plt.scatter(theta, r, c=color, s=area)
11 theta_i = 2*np.pi*(sorted_intensity/255.0)
12 colors = my_sorted_cmap(np.arange(12))
13 c = plt.scatter(theta_i, r, c=colors, s=area)

In line 01 we import the Colors module from the Matplotlib library, then line 02 creates the desired sequence of colours (red, rose, magenta, violet, blue, azure, cyan, aquamarine, electric green, chartreuse, yellow, and orange) using either the name or Hex code, and line 03 generates the colormap. Then we use lines 04 and 05 to assign colours to the chips in the first scatterplot (left), and lines 06, 07, and 09 to specify the number of chips, the angular distances between chips, and the area of the chips, respectively. Line 10 generates the plot. The modifications in lines 11-14 will result in the scatterplot on the right side in the figure (the sorted intensity is calculated in much the same way as in my Geophysical tutorial – How to evaluate and compare colormaps in Python).

[1] Or, perhaps, just create 12 discrete colour classes to group azimuth values in bins of pi/6 (30 degrees) each, and wrap around again at the magenta, to generate a discrete colormap.

[2] The green chip is almost completely covered by the orange chip.

I was curious to see how all these colormaps fared, but my expectation was that Jet would sink to the bottom. I was really surprised to see it came on top, one vote ahead of the linear lightness rainbow (21 and 20 votes out of 62, respectively). The modified heated body followed with 11 votes.

My surprise comes from the fact that Jet carries perceptual artifacts within the progression of colours (see for example this post). One way to demonstrate these artifacts is to convert the 2D map into a 3D surface where again we use Jet to colour amplitude values, but we use the intensities from the 2D map for the elevation. This can be done for example using the Interactive 3D Surface Plot plugin for ImageJ (as in my previous post Lending you a hand with image processing – introduction to ImageJ). The resulting surface is shown in Figure 1. This is almost exactly what your brain would do when you look at the 2D map colored with Jet in the previous post.

Figure 1

In Figure 2 the same data is now displayed as a surface where amplitude values were used for the elevation, with a very light sun shading to help a bit with the perception of relief, but no colormap at all. to When comparing Figure 1 with Figure 2 one of the artifacts is immediately recognized: the highest values in Figure 2, which honours the data, become a relative low in Figure 1. This is because red has lower intensity than yellow and therefore data colored in red in 2D are plotted at a lower elevation than data colored in yellow, even though the amplitudes of the latter were lowest.

Figure 2

For these reasons, I did not expect Jet to be the top pick. On the other hand, I think Jet is perhaps favoured because with consistent use, our brain, learns in part to accommodate for these non-perceptual artifacts in 2D maps, and because it has at least two regions of higher contrast (higher magnitude gradient) than other colormaps. Unfortunately, as I wrote in a recently published tutorial, these regions are randomly placed in the colormap, and the gradients are variable, so we gain on contrast but lose on faithfulness in representing the data structure.

Matt Hall wrote a great comment following the previous post, really making an argument for switching between multiple colormaps in the interpretation stage to explore and highlight features in both the signal and the noise in the data, and that perhaps no single colormap is best overall. I agree 100% on almost everything Matt said, except perhaps on the best overall: looking at the 2D maps, at least with this dataset, I feel the heated body could be the best overall colormap, even if marginally. In Figure 3, Figure 4, Figure 5, and Figure 6 I show the 3D displays obtained by converting the 2D grayscale, linear lightness rainbow, modified heated body, and cube llightness rainbow, respectively. Looking at the 3D displays altogether gives me a confirmation of that feeling.

I recently added to my Matlab File Exchange function, Perceptually improved colormaps, a colormap for periodic data like azimuth or phase. I am going to briefly showcase it using data from my degree thesis in geology, which I used before, for example in Visualization tips for geoscientists – Matlab. Figure 1, from that post, shows residual gravity anomalies in milligals.

Figure 1

Often we’re interested in characterizing these anomalies by calculating the direction of maximum dip at each point on the surface, and for that direction display the azimuth, or dip azimuth. I’ve done this for the surface of residual anomalies from Figure 1 and displayed the azimuth in Figure 2. Azimuth from 0 to 360 degrees are color-coded using Jet, Matlab’s standard colormap (until recently). Typically I do not trust azimuth values when the dip is close to zero because it is often contaminated by noise so I would use shading to de-saturate the colors where dip has the lowest values, but for ease of discussion I haven’t done so in this case.

Figure 2. Azimuth values color-coded with Jet.

There are two problems with Figure 2. First, the well-known problems with the jet colormap. For example, blue is too dark and blue areas appear as bands of constant colour. Yellow is much lighter than any other colour so we see artificial yellow edges that are not really present in the data. But there is an additional issue in Figure 2 because azimuths close in value to 0 and 360 degrees are colored with blue and red, respectively, instead of a single color as they should, causing an additional artificial edge.

In Figure 3 I recolored the map using a colormap that replicates those used in many geophysical software tools to display azimuth or phase data. This is better because it wraps around at 360 degrees but the perceptual issues are unresolved: in this case red, yellow and blue all appear as sharp perceptual edges.

Figure 3. Azimuth values color-coded with generic azimuth colormap.

Figure 4. Azimuth values color-coded with isoluminant azimuth colormap.

In Figure 4 I used my new colormap, called isoAZ (for isoluminant azimuth). This colormap is much better because not only does it wraps around at 360 degrees, but also lightness is held constant for all colors, which eliminates the perceptual anomalies. All the artificial yellow, red, and blue edges are gone, only real edges are left. This can be more easily appreciated in the figure below: if you hover with your mouse over it you are able to switch back and forth between Figure 3 and Figure 4.

From an interpretation point of view, azimuths 180 degrees apart are of opposing colours, which is ideal for dip azimuth data because it allows us to easily recognize folds where dips of opposite direction are juxtaposed at an edge. One example is the sharp edge in the northwest quadrant of Figure 4, where magenta is juxtaposed to green. If you look at Figure 1 you see that there’s a relative high in this area (the edge in Figure 4) with dips of opposite direction on either side (East and West, or 0 and 360 degrees).

The colormap was created in the Lightness-Chroma-Hue color space, a polar transform of the Lab color space, where lightness is the vertical axis and at each value of lightness, chroma is the radial coordinate and hue the polar angle. One limitation of this approach is that due to theirregular shape of the color gamut section at each lightness value, we can never exceed chroma values of about 38-40 (at lightness = 65 in Matlab; in Python, with extensive trial and error, I have not been able to go past 36 using the Scikit-image Color module), which make the resulting colors pale, pastely.

it creates For those that want to experiment with it further, I used just a few lines of code similar to the ones below:

radius = 38; % chroma
theta = linspace(0, 2*pi, 256)'; % hue
a = radius * cos(theta);
b = radius * sin(theta);
L = (ones(1, 256)*65)'; % lightness
Lab = [L, a, b];
RGB=colorspace('RGB<-Lab',Lab(end:-1:1,:));

This code is a modification from an example by Steve Eddins on a post on his Matlab Central blog. In Steve’s example the colormap cycles through the hues as lightness increases monotonically (which by the way is an excellent way to generate a perceptual rainbow). In this case lightness is kept constant and hue cycles through the entire 360 degrees and wraps around. Also, instead of using the Image Processing Toolbox, I used Colorspace, a free function from Matlab File Exchange, for the color space transformations.

For data like fracture orientation where azimuths 180 degrees apart are equivalent it is better to stack two of these isoluminant colormaps in a row. In this way we place opposing colors 90 degrees apart, whereas color 180 degrees apart are the same. You can do it using Matlab commands repmat or vertcat, as below:

radius = 38; % chroma
theta = linspace(0, 2*pi, 128)'; % hue
a = radius * cos(theta);
b = radius * sin(theta);
L = (ones(1, 128)*65)'; % lightness
Lab = [L, a, b];
rgb=colorspace('RGB<-Lab',Lab(end:-1:1,:));
RGB=vertcat(rgb,rgb);

Since then Giuliano has been kind enough to provide me with the data for one of his spectrograms, so I am resuming the discussion. Below here is a set of 5 figures generated in Matlab from the same data using different colormaps. With this post I’d like to get readers involved and ask to cast your vote for the colormap you prefer, and even drop a line in the comments section to tell us the reason for your preference.

In the second post I’ll show the data displayed with the same 5 colormaps but using a different type of visualization, which will reveal what our brain is doing with the colours (without our full knowledge and consent), and then I will ask again to vote for your favourite.

Evan Bianco of Agile Geoscience wrote a wonderful post on how to use python to import, manipulate, and display digital elevation data for Mt St Helens before and after the infamous 1980 eruption. He also calculated the difference between the two surfaces to calculate the volume that was lost because of the eruption to further showcase Python’s capabilities. I encourage readers to go through the extended version of the exercise by downloading his iPython Notebook and the two data files here and here.

I particularly like Evan’s final visualization (consisting of stacked before eruption, difference, and after eruption surfaces) which he created in Mayavi, a 3D data visualization module for Python. So much so that I am going to piggy back on his work, and show how to import a custom palette in Mayavi, and use it to color one of the surfaces.

Python Code

This first code block imports the linear Lightness palette. Please refer to my last post for instructions on where to download the file from.

import numpy as np
# load 256 RGB triplets in 0-1 range:
LinL = np.loadtxt('Linear_L_0-1.txt')
# create r, g, and b 1D arrays:
r=LinL[:,0]
g=LinL[:,1]
b=LinL[:,2]
# create R,G,B, and ALPHA 256*4 array in 0-255 range:
r255=np.array([floor(255*x) for x in r],dtype=np.int)
g255=np.array([floor(255*x) for x in g],dtype=np.int)
b255=np.array([floor(255*x) for x in b],dtype=np.int)
a255=np.ones((256), dtype=np.int); a255 *= 255;
RGBA255=zip(r255,g255,b255,a255)

This code block imports the palette into Mayavi and uses it to color the Mt St Helens after the eruption surface. You will need to have run part of Evan’s code to get the data.

from mayavi import mlab
# create a figure with white background
mlab.figure(bgcolor=(1, 1, 1))
# create surface and passes it to variable surf
surf=mlab.surf(after, warp_scale=0.2)
# import palette
surf.module_manager.scalar_lut_manager.lut.table = RGBA255
# push updates to the figure
mlab.draw()
mlab.show()

You can use the code snippets in here or download the iPython notebook from here (*** please see update ***). You will need NumPy in addition to Matplotlib.

*** UPDATE ***

Fellow Pythonista Matt Hall of Agile Geoscience extended this work – see comment below to include more flexible import of the data and formatting routines, and code to reverse the colormap. Please use Matt’s expanded iPython notebook.

Preliminaries

First of all, get the color palettes in plain ASCII format rom this page. Download the .odt file for the RGB range 0-1 colors, change the file extension to .zip, and unzip it. Next, open the file Linear_L_0-1, which contains comma separated values, replace commas with tabs, and save as .txt file. That’s it, the rest is in Python.

Code snippets – importing data

Let’s bring in numpy and matplotlib:

%pylab inline
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

Now we load the data. We get RGB triplets from tabs delimited text file. Values are expected in the 0-1 range, not 0-255 range.

which is the format required to make the colormap using matplotlib.colors. Here’s the code:

b3=LinL[:,2] # value of blue at sample n
b2=LinL[:,2] # value of blue at sample n
b1=linspace(0,1,len(b2)) # position of sample n - ranges from 0 to 1
# setting up columns for list
g3=LinL[:,1]
g2=LinL[:,1]
g1=linspace(0,1,len(g2))
r3=LinL[:,0]
r2=LinL[:,0]
r1=linspace(0,1,len(r2))
# creating list
R=zip(r1,r2,r3)
G=zip(g1,g2,g3)
B=zip(b1,b2,b3)
# transposing list
RGB=zip(R,G,B)
rgb=zip(*RGB)
# print rgb
# creating dictionary
k=['red', 'green', 'blue']
LinearL=dict(zip(k,rgb)) # makes a dictionary from 2 lists

which gives you this result: Notice that if the original file hadn’t had 256 samples, and you wanted a 256 sample color palette, you’d use the line below instead: