# Spectral lightness rainbow colormap

#### Spectral lightness rainbow

Quick post to share my replica of Ethan Montag ‘s Spectral lightness colormap from this paper. My version has a linear Lightness profile (Figure 1) that increases monotonically from 0 (black) to 100 (white) while Hue cycles from 180 to 360 degrees and then wraps around and continues from 0 to 270.

Figure 1

Figure 2

#### Matlab code

To run this code you will need Colorspace, a free function from Matlab File Exchange, for the color space transformations.

```%% generate Chroma, Lightness, Hue components
h2=linspace(0,270,60*2.56);
h1=linspace(180,360,103);
h=horzcat(h1,h2);
c=ones(1,256)*25;
l=linspace(0,100,256);```
```%% merge together
lch = vertcat(l,c,h)';```
```%% convert to RGB
rgb = colorspace('LCH->rgb',lch);```
```%% Plot color-coded lightness profile
figure;
h=colormapline(1:1:256,lch(:,1),[],rgb);
set(h,'linewidth',2.8);
title('Montag Spectral L* lightness plot ','Color','k','FontSize',12,'FontWeight','demi');```
```%%  Pyramid test data
PY=zeros(241,241);
for i = 1:241
temp=(0:1:i-1);
PY(i,1:i)=temp(:);
end
test=PY.';
test=test-1;
test(test==-1)=0;
test1=test([2:1:end,1],:);
PY1=PY+test1;
PY2=fliplr(PY1);
PY3= flipud(PY1);
PY4=fliplr(PY3);
GIZA=[PY1,PY2;PY3,PY4].*2;
x=linspace(1,756,size(GIZA,1));
y=x;
[X,Y]=meshgrid(x,y);
clear i test test1 PY PY1 PY2 PY3 PY4 temp;

%% display Pyramid surface
fig1=figure;
surf(X,Y,GIZA,GIZA,'FaceColor','interp','EdgeColor','none');
view(-35,70);
colormap(rgb);
axis off;
grid off;
colorbar;
% set(fig1,'Position',[720 400 950 708]);
% set(fig1,'OuterPosition',[716 396 958 790]);
title('Montag Spectral L*','Color','k','FontSize',12,'FontWeight','demi');```

#### Aknowledgements

The coloured lightness profiles were made using the Colormapline submission from the Matlab File Exchange.

# Colormap compromise

At the end the series The rainbow is dead…long live the rainbow! I introduced CubicYF, a rainbow-like color palette with a slightly compressive, monotonic increasing lightness profile. The red color is missing from the palette because green and red at people with deuteranopia, the most common type of color vision deficiency, confuse red and green, so I took one out. An alternative solution to taking the red out is to have a yellow-like color with lower lightness followed by a red-like color of higher lightness, like in the LinearL palette described The rainbow is dead…long live the rainbow! – Part 5. A third solution, which is a compromise in perceptual terms, is CubicL (sometimes I called it cube1), the color palette I used with the gravity data from my degree thesis in geology in the series Visualization tips for geoscientists. An example map from the series in which I used CubicL is in Figure 1 below:

Figure 1

To generate CubicL  I used a compressive function for the increase in lightness very much like the one for CubicYF, except in this case I allowed lightness to decrease from about 90 in the yellow portion of the colormap to about 80 to get to the red. As a result, there is an inversion in the lightness trend, indicated by the arrow in the bottom row in Figure 2. However, I believe this is an acceptable compromise, and one that is pleasing for people accustomed to rainbow (it has red), because the inversion is smooth, with a gentle downward-facing concavity. I also believe this colormap is still a perceptual improvement over the rainbow or spectrum so commonly used in seismic visualization, which typically have at least 3 sharp lightness inversions (indicated by arrows in the top row of Figure 2 below).

Figure 2

CubicL is one of the colormaps available in my Matlab File Exchange function Perceptually Improved Colormaps.

# Convert color palettes to python matplotlib colormaps

This is a quick post to show you how to import my perceptual color palettes – or any other color palette – into Python and convert them into Matplotlib colormaps. We will use as an example the CIE Lab linear L* palette, which was my adaptation to Matlab of the luminance controlled colormap by Kindlmann et al.

#### Introduction

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.

`LinL = np.loadtxt('Linear_L_0-1.txt')`

A quick QC of the data:

`LinL[:3]`

gives us:

```array
([[ 0.0143,  0.0143,  0.0143],
[ 0.0404,  0.0125,  0.0325],
[ 0.0516,  0.0154,  0.0443]])```

looking good!

#### Code snippets – creating colormap

Now we want to go from an array of values for blue in this format:

`b1=[0.0000,0.1670,0.3330,0.5000,0.6670,0.8330,1.0000]`

to a list that has this format:

```[(0.0000,1.0000,1.0000),
(0.1670,1.0000,1.0000),
(0.3330,1.0000,1.0000),
(0.5000,0.0000,0.0000),
(0.6670,0.0000,0.0000),
(0.8330,0.0000,0.0000),
(1.0000,0.0000,0.0000)]```

# to a dictionary entry that has this format:

```{
'blue': [
(0.0000,1.0000,1.0000),
(0.1670,1.0000,1.0000),
(0.3330,1.0000,1.0000),
(0.5000,0.0000,0.0000),
(0.6670,0.0000,0.0000),
(0.8330,0.0000,0.0000),
(1.0000,0.0000,0.0000)],
...
...
}```

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

```

#### Code snippets – testing colormap

That’s it, with the next 3 lines we are done:

```my_cmap = matplotlib.colors.LinearSegmentedColormap('my_colormap',LinearL)
pcolor(rand(10,10),cmap=my_cmap)
colorbar()```

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:

`my_cmap = matplotlib.colors.LinearSegmentedColormap('my_colormap',LinearL,256)`

Have fun!

#### Related Posts (MyCarta)

The rainbow is dead…long live the rainbow! – Part 5 – CIE Lab linear L* rainbow

Visualize Mt S Helens with Python and a custom color palette

#### A bunch of links on colors in Python

http://bl.ocks.org/mbostock/5577023

http://matplotlib.org/api/cm_api.html

http://matplotlib.org/api/colors_api.html

http://matplotlib.org/examples/color/colormaps_reference.html

http://matplotlib.org/api/colors_api.html#matplotlib.colors.ListedColormap

http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps

http://stackoverflow.com/questions/21094288/convert-list-of-rgb-codes-to-matplotlib-colormap

http://stackoverflow.com/questions/16834861/create-own-colormap-using-matplotlib-and-plot-color-scale

http://stackoverflow.com/questions/11647261/create-a-colormap-with-white-centered-around-zero

# Perceptual rainbow palette – the method

With this post I would like to introduce my new, perceptually balanced rainbow color palette. I used the palette for the first time in How to assess a colourmap, an essay I wrote for 52 Things You Should Know About Geophysics, edited by Matt Hall and Evan Bianco of Agile Geoscience.

In my essay I started with the analysis of the spectrum color palette, the default  in some seismic interpretation softwares, using my Lightness L* profile plot and Great Pyramid of Giza test surface (see this post for background on the tests and to download the Matlab code). The profile and the pyramid are shown in the top left image and top right image in Figure 1, from the essay.

Figure 1

In the plot the value of L* varies with the color of each sample in the spectrum, and the line is colored accordingly. This erratic profile highlights several issues with spectrum: firstly, the change in lightness is not monotonic. For example it increases from black (L*=0) to magenta [M] then drops from magenta to blue [B], then increases again and so on. This is troublesome if spectrum is used to map elevation because it will interfere with the correct perception of relief, particularly if shading is added. Additionally, the curve gradient changes many times, indicating a nonuniform perceptual distance between samples. There are also plateaus of nearly flat L*, creating bands of constant color (a small one at the blue, and a large one at the green [G]).

The Great Pyramid has monotonically increasing elevation (in feet – easier to code) so there should be no discontinuities in the surface if the color palette is perceptual. However, clearly using the spectrum we have introduced many artificial discontinuities that are not present in the data. For the bottom row in FIgure 1 I used my new color palette, which has a nice, monotonic, compressive Lightness profile (bottom left). Using this palette the pyramid surface (bottom right) is smoothly colored, without any perceptual artifact.

This is how I created the palette: I started with RGB triplets for magenta, blue, cyan, green, and yellow (no red), which I converted to L*a*b* triplets using Colorspace transformations, a Matlab function available on the Matlab File Exchange. I modified the new L* values by fitting them to an approximately cube law L* function (this is consistent with Stevens’ power law of perception), and adjusted a* and b* values using Lab charts like the one in Figure 2 (from CIELab Color Space by Gernot Hoffmann, Department of Mechanical Engineering, University of Emden)  to get 5 colors moving up the L* axis along an imaginary spiral (I actually used tracing paper). Then I interpolated to 256 samples using the same ~cube law, and finally reconverted to RGB [1].

Figure 2

There was quite a bit of trial and error involved, but I am very happy with the results. In the animations below I compare the spectrum and the new palette, which I call cubeYF, as seen in CIELab color space. I generated these animations with the method described in this post, using the 3D color inspector plugin in ImageJ:

I also added Matlab’s default Jet rainbow – a reminder that defaults may be a necessity, but in many instances not the ideal choice:

OK, the new palette looks promising, insofar as modelling is concerned. But how would it fare using some real data? To answer this question I used a residual gravity map from my unpublished thesis in Geology at the University of Rome. I introduced this map and discussed the geological context and objectives of the geophysical study in a previous post, so please refer to that if you are curious about it. In this post I will go straight to the comparison of the color palettes; if you are unfamiliar with gravity data, try to imagine negative residuals as elevation below sea level, and positive residuals as elevation above seal level – you won’t miss out on anything.

In Figures 3 to 6 I colored the data using the above three color palettes, and grayscale as benchmark. I generated these figures using Matlab code I shared in my post Visualization tips for geoscientists: Matlab, and I presented three of them (grayscale, Spectrum, and cubeYF) at the 2012 convention of the Canadian Society of Exploration Geophysicists in Calgary (the extended abstract, which I co-authored with Steve Lynch of 3rd Science, is available here).

In Figure 3, the benchmark for the following figures, I use grayscale to represent the data, assigning increasing intensity from most negative gravity residuals in black to most positive residuals in white (as labeled next to the colorbar). Then, I used terrain slope to create shading: the higher the slope, the darker the shading that is assigned, which results in a pseudo-3D display that is very effective (please refer to Visualization tips for geoscientists: Surfer, for an explanation of the method, and Visualization tips for geoscientists: Matlab for code).

Figure 3 – Grayscale benchmark

In Figure 4 I color the pseudo-3D surface with the cubeYF rainbow. Using this color palette instead of grayscale allows viewers to appreciate smaller changes, more quickly assess differences, or conversely identify areas of similar anomaly, while at the same time preserving the peudo-3D effect. Now compare Figure 4 with Figure 5, where we use the spectrum to color the surface: this palette introduces several artefacts (sharp edges and bands of constant hue) which confuse the display and interfere with the perception of pseudo-relief, all but eliminating the effect.  For Figure 6 I used Matlab’s default Jet color palette, which is better that the spectrum, and yet the relief effect is somewhat lost (due mainly to a sharp yellow edge and cyan band).

Figure 4 – cube YF rainbow

Figure 5 – Industry spectrum

Figure 6 – Matlab Jet

It looks like both spectrum and jet are poor choices when used for color representation of a surface, with the new color palette a far superior alternative. In the CSEG convention paper mentioned above (available here) Steve and I went further by showing that the spectrum not only has these perceptual artifacts and edges, but it is also very confusing for viewers with deficient color vision, a condition that occurs in about 8% of Caucasian males. We did that using computer software [2] to simulate how viewers with two types of deficient color vision, Deuteranopia and Tritanopia, would see the two colored surfaces, and we compare the results. In other words, we are now able to see the images as they would see them. Please refer to the paper for a full discussion on these simulation.

In here, I show in Figures  7 to 9 the Deuteranope simulations for cubeYF, spectrum, and jet, respectively. In all three simulations the hue discrimination has decreased, but while the spectrum and jet are now even more confusing, the cubeYF has preserved the relief effect.

Deuteranope Simulation of cube YF

Deuteranope Simulation of Industry spectrum

Deuteranope Simulation of Matlab Jet

That’s it for today. In my next post, to be published very shortly, you will get the palette, and a lot more.

#### References

A more perceptual color palette for structure maps, CSEG/CSPG 2012 convention, Calgary

How to assess a colourmap, in 52 Things You Should Know About Geophysics

#### Notes

[1] An alternative to the method I used would be to start directly in CIELab color space, and use a some kind of spiral *L lightness profile programmatically.  For example:

– Using 3D helical curves from: http://www.mathworks.com/matlabcentral/fileexchange/25177-3d-curves

– Using Archimedes spiral

– Expanding on code by Steve Eddins at Mathworks (A path through L*a*b* color space) in this article , one could create a spiral cube lightness with something like:

```%% this creates best-fit pure power law function
%  Inspired by wikipedia - http://en.wikipedia.org/wiki/Lightness
l2=linspace(1,power(100,0.42),256);
L2=(power(l2,1/0.42))';

%% this makes cielab real cube function spiral
theta = linspace(0.6*pi, 2*pi, 256).';
Lab1 = [L2, a, b]; RGB_realcube=colorspace('RGB<-Lab',(Lab1));```

[2] The simulations are created using ImageJ, an open source image manipulation program, and the Vischeck plug-in. I later discovered Dichromacy, anther ImageJ plug-in for these simulations, which has the advantage of being an open source plugin. They can also be performed on the fly (no upload needed) using the online tool Color Oracle.

# The rainbow is dead…long live the rainbow! – Perceptual palettes, part 5 – CIE Lab linear L* rainbow

#### Some great examples

After my previous post in this series there was a great discussion on perceptual color palettes with some members of the Worldwide Geophysicists group on LinkedIn. Ian MacLeod shared some really good examples, and uploaded it in here.

#### HSL linear L rainbow palette

Today I’d like to share a color palette that I really like:

It is one of the palettes introduced in a paper by Kindlmann et al. [1]. The authors created their palettes with a technique they call luminance controlled interpolation. They explain it in this online presentation. However they used different palettes (their isoluminant rainbow, and their heated body) so if you find it confusing I recommend you look at the paper first. Indeed, this is a good read if you are interested in colormap generation techniques; it is one of the papers that encouraged me to develop the methodology for my cube law rainbow, which I will introduce in an upcoming post.

This is how I understand their method to create the palette: they mapped six pure-hue rainbow colors (magenta, blue, cyan, green, yellow, and red) in HSL space, and adjusted the Luminance by changing the HSL Lightness value to ‘match’ that of six control points evenly spaced along the gray scale palette. After that, they interpolated linearly along the L axis between 0 and 1 using the equation presented in the paper.

#### CIE Lab linear L* rainbow palette

For this post I will try to create a similar palette. In fact, initially I was thinking of just replicating it, so I imported the palette as a screen capture image into Matlab, reduced it to a 256×3 RGB colormap matrix, and converted RGB values to Lab to check its linearity in lightness. Below I am showing the lightness profile, colored by value of L*, and the Great Pyramid of Giza – my usual test surface –  also colored by L* (notice I changed the X axis of both L* plots from sample number to Pyramid elevation to facilitate comparison of the two figures).

Clearly, although the original palette was constructed to be perceptually linear, it is not linear following my import. Notice in particular the notch in the profile in the blue area, at approximately 100 m elevation. This artifact is also visible as a flat-looking blue band in the pyramid.

I have to confess I am not too sure why the palette has this peculiar lightness profile. I suspect this may be because their palette is by construction device dependent (see the paper) so that when I took the screen capture on my monitor I introduced the artifacts.

The only way to know for sure would be to use their software to create the palette, or alternatively write the equation from the paper into Matlab code and create a palette calibrated on my monitor, then compare it to the screen captured one. Perhaps one day I will find the time to do it but having developed my own method to create a perceptual palette my interest in this one became just practical: I wanted to get on with it and use it.

#### Fixing and testing the palette

Regardless of what the cause might be for this nonlinear L* profile, I decide to fix it and I did it by simply replacing the original profile with a new one, linearly changing between 0.0 and 1.0. Below I am showing the L* plot for this adjusted palette, and the Great Pyramid of Giza, both again colored by value of L*.

The pyramid with the adjusted palette seems better: the blue band is gone, and it looks great. I am ready to try it on a more complex surface. For that I have chosen the digital elevation data for South America available online through the Global Land One-km Base Elevation Project at the National Geophysical Data Center. To load and display the data in Matlab I used the first code snippet in Steve Eddin’s post on the US continental divide  (modified for South America data tiles). Below is the data mapped using the adjusted palette. I really like the result: it’s smooth and it looks right.

But how do I know, really? I mean, once I move away from my perfectly flat pyramid surface, how do I know what to expect, or not expect? In other words, how would I know if an edge I see on the map above is an artifact, or worse, that the palette is not obscuring real edges?

In some cases the answer is simple. Let’s take a look at the four versions of the map in my last figure. The first on the left was generated using th ROYGBIV palette I described in this post. It would be obvious to me, even if I never looked at the L* profile, that the blue areas are darker than the purple areas, giving the map a sort of inverted image look.

But how about the second map from the left? For this I used the default rainbow from a popular mapping program. This does not look too bad at first sight. Yes, the yellow is perceived as a bright, sharp edge, and we now know why that is, but other than that it would be hard to tell if there are artifacts. After a second look the whole area away from the Andes is a bit too uniform.

A good way to assess these maps is to use grayscale, which we know is a good perceptual option, as a benchmark. This is the last map on the right. The third map of South America was coloured using my adjusted linear L* palette. This maps looks more similar to our grayscale benchmark. Comparison of the colorbars will also help: the third and fourth are very similar and both look perceptually linear, whereas the third does show flatness in the blue and green areas.

Let me know what you think of these examples. And as usual, you are welcome to use the palette in your work. You can download it here.

#### UPDATE

With my following post, Comparing color palettes, I introduced my new method to compare palettes with ImageJ and the 3D color inspector plugin. Here below are the recorded 3D animations of the initial and adjusted palettes respectively. In 3D it is easier to see there is an area of flat L* between the dark purple and dark blue in the initial color palette. The adjusted color palette instead monotonically spirals upwards.

#### References

[1] Kindlmann, G. Reinhard, E. and Creem, S., 2002, Face-based Luminance Matching for Perceptual Colormap Generation, IEEE – Proceedings of the conference on Visualization ’02

#### Related posts (MyCarta)

The rainbow is dead…long live the rainbow! – the full series

What is a colour space? reblogged from Colour Chat

Color Use Guidelines for Mapping and Visualization

A rainbow for everyone

Is Indigo really a colour of the rainbow?

Why is the hue circle circular at all?

A good divergent color palette for Matlab

#### Related topics (external)

Color in scientific visualization

The dangers of default disdain

#### Color tools

How to avoid equidistant HSV colors

Colormap tool

Color Oracle – color vision deficiency simulation – stand alone (Window, Mac and Linux)

Dichromacy –  color vision deficiency simulation – open source plugin for ImageJ

Vischeck – color vision deficiency simulation – plugin for ImageJ and Photoshop (Windows and Linux)

#### For teachers

NASA’s teaching resources for grades 6-9: What’s the Frequency, Roy G. Biv?

#### ImageJ and 3D Color inspector plugin

http://rsbweb.nih.gov/ij/docs/concepts.html

http://rsb.info.nih.gov/ij/plugins/color-inspector.html

# Visualization tips for geoscientists: Matlab, part III

#### Introduction

Last weekend I had a few hours to play with but needed a short break from writing about color palettes, so I decided to go back and finish up (for now) this series on geoscience visualization in Matlab. In the first post of the series I expanded on work by Steve Eddins at Mathworks on overlaying images using influence maps and demonstrated how it could be used to enhance the display of a single geophysical dataset.

#### Using transparency to display multiple data sets an example

At the end of the second post I promised I would go back and show an example of using transparency and influence maps for other tasks, like overlaying of different attributes. Here’s my favorite example in Figure 1. The image is a map in pastel colors of the Bouguer Gravity anomaly for the Southern Tuscany region of Italy, with three other layers superimposed using the techniques above mentioned.

It is beyond the objectives of this post to discuss at length about gravity exploration methods or to attempt a full interpretation of the map. I will go back to it at a later time as I am planning a full series on gravity exploration using this data set, but if you are burning to read more about gravity interpretation please check these excellent notes by Martin Unsworth, Professor of Physics at the Earth and Atmospheric Sciences department, University of Alberta, and note 4 at the end of this post. Otherwise, and for now, suffice it to say that warm colors (green to yellow to red) in the Bouguer gravity map indicate, relatively speaking, excess mass in the subsurface and blue and purple indicate deficit of mass in the subsurface.

The black and grey lines are lineaments extracted from derivatives of the Bouguer gravity data using two different methods [1]. The semitransparent, white-filled polygons show the location of some of the  basement outcrops (the densest rocks in this area).

Lineaments extracted from gravity data can correspond to contacts between geological bodies of different density, so a correlation can be expected between basement outcrops and some of the lineaments, as they are often placed in lateral contact with much lesser dense rocks. This is often exploited in mineral exploration in areas such as this where mineralization occurs at or in the vicinity of this contacts. As an example, I show in Figure 2 the occurrences (AGIP – RIMIN, unpublished industry report, 1989) of silicization (circles) and antimony deposits (triangles), superimposed on the distance from one of the set of lineaments (warm colors indicate higher distance) from Figure 1.

The fact that different methods give systematically shifted results is a known fact, often due the trade-off between resolution and stability, whereby the more stable methods are less affected by noise, but often produce smoother edges over deeper contacts, and their maxima may not correspond. This is in addition to the inherent ambiguity of gravity data, which cannot, by themselves, be interpreted uniquely. To establish which method might be more correct in this case (none is a silver bullet) I tried to calibrate the results using basement outcrops (i.e. does either method more closely match the outcrop edges?). Having done that, I would have more confidence in making inferences on possible other contacts in the subsurface suggested by lineament. I would say the black lines do a better overall job in the East, the gray perhaps in the West. So perhaps I’m stuck? I will get back to this during my gravity series.

Figure 1

Figure 2

#### Matlab code

As usual I am happy to share the code I used to make the combined map of Figure 1. Since the data I use is in part from my unpublished thesis in Geology and in part from Michele di Filippo at the University of Rome, I am not able to share it, and you will have to use your own data, but the Matlab code is simply adapted. The code snippet below assume you have a geophysical surface already imported in the workspace and stored in a variable called “dataI”, as well as the outcrops in a variable called “basement”, and the lineaments in “lnmnt1” and “lnmnt2”. It also uses my cube1 color palette.

```
% part 1 - map gravity data
figure; imagesc(XI,YI,dataI); colormap(cube1); hold on;
%
% part 2 - dealing with basement overlay
white=cat(3, ones(size(basement)), ones(size(basement)),...
ones(size(basement)));
ttt=imagesc(Xb,Yb,white); % plots white layer for basement
%
% part 3 - dealing with lineaments overlays
black=cat(3, zeros(size(lnmnt1)), zeros(size(lnmnt1)),...
zeros(size(lnmnt1)));
grey=black+0.4;
basement_msk=basement.*0.6;
kkk=imagesc(XI,YI,black); % plots black layer for lineament 1
sss=imagesc(XI,YI,gray); % plots gray layer for lineament 2
hold off
%
% part 4 - set influence maps
set(ttt, 'AlphaData', basement_msk); % influence map for basement
set(kkk, 'AlphaData', lnmnt1); % influence map for linement 1
set(sss, 'AlphaData', lnmnt2); % influence map for linement 2
%
% making it pretty
axis equal
axis tight
axis off
set(gca,'YDir','normal');
set(gcf,'Position',[180 150 950 708]);
set(gcf,'OuterPosition',[176 146 958 790]);```

#### Matlab code, explained

OK, let’s break it down starting from scratch. I want first to create a figure and display the gravity data, then hold it so I can overlay the other layers on top of it. I do this with these two commands:

figure;imagesc(XI,YI,dataI);

hold on;

The layer I want to overlay first is the one showing the basement outcrops. I make a white basement layer covering the full extent of the map, which is shown in Figure 3, below.

Figure 3

I create it and plot it with the commands:

white=cat(3, ones(size(basement)), ones(size(basement)), ones(size(basement)));

ttt=imagesc(Xb,Yb,white);

The handle  ttt is to be used in combination with the basement influence map to produce the partly transparent basement overlay: remember that I wanted to display the outcrops in white color, but only partially opaque so the colored gravity map can still be (slightly) seen underneath. I make the influence map, shown in Figure 4, with the command:

basement_msk=basement.*0.6;

Since the original binary variable “basement” had values of 1 for the outcrops and 0 elsewhere, whit the command above I assign an opacity of 0.6 to the outcrops, which will be applied when the next command, below, is run, achieving the desired result.

set(ttt, ‘AlphaData’, basement_msk); % uses basement influence map

Figure 4

For the lineaments I do things in a similar way, except that I want those plotted with full opacity since they are only 1 pixel wide.

As an example I am showing in Figure 5 the black layer lineament 1 and in Figure 6 the influence map, which has values of 1 (full opacity) for the lineament and 0 (full transparency) for everywhere else.

Figure 5

Figure 6

Now a few extra lines to make things pretty, and this is what I get, shown below in Figure 7: not what I expected!

Figure 7

The problem is in these two commands:

white=cat(3, ones(size(basement)), ones(size(basement)), ones(size(basement)));

ttt=imagesc(Xb,Yb,white);

I am calling the layer white but really all I am telling Matlab is to create a layer with maximum intensity (1). But the preceding colormap(cube1) command assigned a salmon-red color to the maximum intensity in the figure, and so that is what you get for the basement overlay.

Again, to get the result I wanted, I had to come up with a trick like in the second post examples. This is the trick:

I create a new color palette with this command:

cube1edit=cube1; cube1edit(256,:)=1;

The new color palette has last RGB triplet actually defined as white, not salmon-red.

Then I replace this line:

figure; imagesc(XI,YI,dataI); colormap(cube1); hold on;

with the new line:

figure; imagesc(XI,YI,dataI, [15 45]); colormap (cube1edit); hold on;

The highest value in dataI is around 43. By spreading the color range from [15 43] to [15 45], therefore exceeding max(dataI) I ensure that white is used for the basement overlay but not in any part of the map where gravity is highest but there is no basement outcrop. In other words, white is assigned in the palette but reserved to the overlay.

Please let me know if that was clear. If it isn’t I will try to describe it better.

#### Notes

[1] One method is the total horizontal derivative. The other method is the hyperbolic tilt angle – using Matlab code by Cooper and Cowan (reference). This is how I produced the two overlays:  first I calculated the total horizontal derivative and the tilt angle, then I found the maxima to use as the overlay layers. This is similar to Figure 3e in Cooper and Cowan, but I refined my maxima result by reducing them to 1-pixel-wide lines (using a thinning algorithm).

#### Reference

Cooper, G.R.J., and Cowan, D.R. (2006) – Enhancing potential field data using filters based on the local phase  Computers & Geosciences 32 (2006) 1585–1591

#### Related posts (MyCarta)

Visualization tips for geoscientists: Surfer

Visualization tips for geoscientists: Matlab

Visualization tips for geoscientists: Matlab, part II

Image Processing Tips for Geoscientists – part 1

# Visualization tips for geoscientists – Matlab

#### Introduction

In my last post I described how to create a powerful, nondirectional shading for a geophysical surface using the slope of the data to assign the shading intensity (i.e. areas of greater slope are assigned darker shading). Today I will show hot to create a similar effect in Matlab.

Since the data set I use is from my unpublished thesis in Geology, I am not able to share it, and you will have to use your own data, but the Matlab code is simply adapted. The code snippets below assume you have a geophysical surface already imported in the workspace and stored in a variable called “data”, as well as the derivative in a variable called “data_slope”.

#### Method 1 – with a slope mask and transparency

Some time ago I read this interesting Image Processing blog post by Steve Eddins at Mathworks on overlaying images using transparency. I encourage readers to take a look at this and other posts by Steve, he’s great! That particular blog post gave me the idea to use transparency and the slope to create my favorite shading in Matlab.

In addition to the code below you will need normalise.m from Peter Kovesi‘s website, and to import the color palette cube1.

```%% alpha transparency code snippet
black = cat(3, zeros(size(data)), zeros(size(data)), ...
zeros(size(data)));             % make a truecolor all-black image
gray=black+0.2;                     % make a truecolor all-gray image
alphaI=normalise(data_slope);       % create transparency weight matrix
% using data_slope

imagesc(data);colormap(cube1);      % display data
hold on
h = imagesc(gray);                  % overlay gray image on data
hold off
set(h, 'AlphaData', alphaI);        % set transparency of gray layer using
axis equal;                         % weight matrix
axis tight;
axis off;```

And here is the result in Figure 1 below – not bad!

# Visualization tips for geoscientists: Surfer

I have been thinking for a while about writing on visualization of geophysical data. I finally got to it, and I am now pleased  to show you a technique I use often.  This tutorial has shaped up into 2 independent posts: in the first post I will show how to implement the technique with Surfer, in the second one with Matlab (you will need access to a license of Surfer 8.08 or later, and Matlab 2007a or later to replicate the work done in the tutorial).

I will illustrate the technique using gravity data since it is the data I developed it for. In an upcoming series of gravity exploration tutorials I will discuss in depth the acquisition, processing, enhancement, and interpretation of gravity data (see [1] and [4]). For now, suffice it to say that gravity prospecting is useful in areas where rocks with different density are laterally in contact, either stratigraphic or tectonic, producing a measurable local variation of the gravitational field. This was the case for the study area (in the Monti Romani of Southern Tuscany) from my thesis in Geology at the University of Rome [2].

In this part of the Apennine belt, a Paleozoic metamorphic basement (density ~2.7 g/cm3) is overlain by a thick sequence of clastic near-shore units of the Triassic-Oligocene Tuscany Nappe (density ~2.3 g/cm3). The Tuscan Nappe is in turn covered by the Cretaceous-Eocene flish units of the Liguride Complex (density ~2.1 g/cm3).

During the deformation of the Apennines, NE verging compressive thrusts caused doubling of the basement. The tectonic setting was later complicated by tensional block faulting with formation of horst-graben structures generally extend along NW-SE and N-S trends which were further disrupted by later and still active NE-SW normal faulting (see [2], and reference therein, for example [3]).

This complex tectonic history placed the basement in lateral contact with the less dense rocks of the younger formations and this is reflected in the residual anomaly map [4] of Figure 1. Roughly speaking, there is a high in the SE quadrant of ~3.0 mgal corresponding to the location of the largest basement outcrop, an NW-SE elongated high of ~0.5 mgal in the centre bound by lows on both the SW and NE (~-6.0 and ~-5.0 mgal, respectively), and finally a local high in the N.W. quadrant of ~-0.5 mGal. From this we can infer that in this area can infer that the systems of normal faults caused differential sinking of the top of basement in different blocks leaving an isolated high in the middle, which is consistent with the described tectonic history [2]. Notice that grayscale representation is smoothly varying, reflecting (and honoring) the structure inherent in the data. It does not allow good visual discrimination and comparison of differences, but from the interpretation standpoint I recommend to always start out with it: once a first impression is formed it is difficult to go back. There is time later to start changing the dispaly.

Figure 1 – Grayscale residual anomalies in milligals. This version of the map was generate using the IMAGE MAP option in Surfer.

OK, now that we formed a first impression, what can we try to improve on this display? The first thing we can do to increase the perceptual contrast is to add color, as I have done in Figure 2. This is an improvement, now we are able to appreciate smaller changes, quickly assess differences, or conversely identify areas of similar anomaly. Adding the third dimension and perspective is a further improvement, as seen in figure 3. But there’s still something missing. Even though we’ve added color, relief, and perspective, the map looks a bit “flat”.

Figure 2 – Colored residual anomalies in milligals. This version of the map was generate using the IMAGE MAP option in Surfer.
Figure 3 – Colored 3D residual anomaly map in milligals. This version of the map was generate using the SURFACE MAP option in Surfer.

Adding contours is a good option to further bring out details in the data, and I like the flexibility of contours in Surfer. For example, for Figure 4 I assigned (in Contour Properties, Levels tab) a dashed line style to negative residual contours, and a solid line style to positive residual contours, with a thicker line for the zero contour. This can be done by modifying the style for each level individually, or by creating two separate contours, one for the positive data, one for the negative data, which is handy when several contour levels are present. The one drawback of using contours this way is that it is redundant. We used 3 weapons  – color, relief, and contours – to dispaly one dataset, and to characterize just one property, the shape of gravity anomaly. In geoscience it is often necessary, and desireable to show multiple (relevant) datasets in one view, so this is a bit of a waste. I would rather spare the contours, for example, to overlay and compare anomalous concentrations in gold pathfinder elements on this gravity anomaly map (one of the objectives of the study, being the Monti Romani an area of active gold exploration).

Figure 4 – Colored 3D residual anomaly map in milligals. Contours were added with the the CONTOUR MAP option in Surfer.  Figure 5 – Colored 3D residual anomaly map in milligals with lighting (3D Surface Properties menu). Illumination is generated by a point source with -135 deg azimuth and 60 deg elevation, plus an additional 80% gray ambient light, a 30% gray diffuse light, and a 10% gray specular light.

The alternative to contours is the use of illumination, or lighting, which I used in Figure 5. Lighting is doing a really good job: now we can recognize there is a high frequency texture in the data and we see some features both in the highs and lows. But there’s a catch:  we are now introducing perceptual artifacts, in the form of bright white highlight, which is obscuring some of the details where the surface is orthogonal to the point source light.

There is a way to illuminate the surface without introducing artifact – and that is really wanted to show you with this tutorial – which is to use a derivative of the data to assign the shading intensity (areas of greater gradient were assigned darker shading) [5]. In this case  I choose the terrain slope, which is the slope in the direction of steepest gradient at any point in the data (calculated in a running window). The result is a very powerful shading. Here is how you can do it in Surfer:

1) CREATE TERRAIN SLOPE GRID (let’s call this A): go to GRID > CALCULUS > TERRAIN SLOPE

Result is shown in Figure 6 below:

Figure 6 – Terrain slope of residual anomaly. Black for low gradients, white for high gradients. Displayed using IMAGE MAP option.

2) CREATE COMPLEMENT OF TERRAIN SLOPE AND NORMALIZE TO [1 0] RANGE (to assign darker shading to areas of greater slope. This is done with 3 operations:

i)     GRID > MATH> B=A – min(A)

where min(A) is the minimum value, which you can read off the grid info (for example you would double click on the map above to open the Map Properties and there’s an info button next to the Input File field) .

ii)    GRID > MATH> C=B /max(B)

iii)   GRID > MATH> D= 1-C

Result is shown in Figure 7 below. This looks really good, see how now the data seems almost 3D? It would work very well just as it is. However, I do like color, so I added it back in Figure 8. This is done by draping the grayscale terrain slope complement IMAGE MAP as an overlay over the colored residual anomaly SURFACE MAP, and setting the Color Modulation to BLEND in the 3D Surface Properties in the Overlay tab. I really do like this display in Figure 8, I think it is terrific. Let me know if you like it best too.

Finally, in Figure 9, I added a contour of the anomaly in the Gold Pathfiners, to reiterate the point I made above that contours are best spared for a second dataset.

In my next post I will show you how to do all of the above programmatically in Matlab (and share the code). Meanwhile, comments, suggestions, requests are welcome. Have fun mapping and visualizing!

Figure 7 – Complement of the terrain slope. White for low gradients, black for high gradients. Displayed using IMAGE MAP option.

Figure 8 – Complement of the terrain slope with color added back.
Figure 9 – Complement of the terrain slope with color added back and contour overlay of gold pathfinders in stream sediments.

#### GOODIES

Did you lie the colormap? In future series on perceptually balanced colormaps I will tell you how I created it. For now, if you’d like to try it on your data you can download it here:

cube1 – generic format with RGB triplets;

Cube1_Surfer – this is preformatted for Surfer with 100 RGB triplets and header line. Dowload the .doc file, open and save as plain text, then change extension to .clr;

Cube1_Surfer_inverse – the ability to flip color palette is not implemented in Surfer (at least not in version 8) so I am including the flipped version of above. Again, dowload the .doc file, open and save as plain text, then change extension to .clr.

#### RELATED POSTS (MyCarta)

Visualization tips for geoscientists: Matlab

Visualization tips for geoscientists: Matlab, part III

Image Processing Tips for Geoscientists – part 1

Compare lists from text files using Matlab – an application for resource exploration

#### RELATED CONTENT (External)

Basement structure in central and southern Alberta: insights from gravity and magnetic maps

Making seismic data come alive

Visual Crossplotting

Mapping multiple attributes to three- and four-component color models — A tutorial

Enhancing Fault Visibility Using Bump Mapped Seismic Attributes

#### ACKNOWLEDGEMENTS

I would like to thank Michele di Filippo at the Department of Earth Science, University of Rome La Sapienza, to whom I owe a great deal. Michele, my first mentor and a friend, taught me everything I know about the planning and implementation of a geophysical field campaign. In the process I also learned from him a great deal about geology, mapping, Surfer, and problem solving. Michele will make a contribution to the gravity exploration series.

###### NOTES

[1] If you would like to learn more about gravity data interpretation please check these excellent notes by Martin Unsworth Unsworth, Professor of Physics at the Earth and Atmospheric Sciences department, University of Alberta.

[2] Niccoli, M., 2000:  Gravity, magnetic, and geologic exploration in the Monti Romani of Southern Tuscany, unpublished field and research thesis, Department of Earth Science, University of Rome La Sapienza.

[3] Moretti A., Meletti C., Ottria G. (1990) – Studio stratigrafico e strutturale dei Monti Romani (GR-VT) – 1: dal Paleozoico all’Orogenesi Alpidica. Boll. Soc. Geol. It., 109, 557-581. In Italian.

[4] Typically reduction of the raw data is necessary before any interpretation can be attempted. The result of this process of reduction is a Bouguer anomaly map, which is conceptually equivalent to what we would measure if we stripped away everything above sea level, therefore observing the distribution of rock densities below a regular surface. It is standard practice to also detrend the Bouguer anomaly to separate the influence of basin or crustal scale effects, from local effects, as either one or the other is often the target of the survey. The result of this procedure is typically called Residuals anomaly and often shows subtler details that were not apparent due to the regional gradients. Reduction to rsiduals makes it easier to qualitatively separate mass excesses from mass deficits. For a more detailed review of gravity exploration method check again the notes in [1] and refer to this article on the CSEG Recorder and reference therein.

[5] Speaking in general, 3D maps without lighting often have a flat appearance, which is why light sources are added. The traditional choice is to use single or multiple directional light sources, but the result is that only linear features orthogonal to those orientations will be highlighted. This is useful when interpreting for lineaments or faults (when present), but not in all circumstances, and requires a lot of experimenting. in other cases, like this one , directional lighting introduces a bright highlight, which obscures some detail. A more generalist, and in my view more effective alternative, is to use information derived from the data itself for the shading. One way to do that is to use a high pass filtered version of the data. i will show you how to do that in matlab in the next tutorial. Another solution, which I favored in this example, is to use a first derivative of the data.