Issue creating map shading in Matplotlib imshow by setting opacity to data gradient

I am trying to add shading to a map of some data by calculating the gradient of the data and using it to set alpha values.

I start by loading my data (unfortunately I cannot share the data as it is being used in a number of manuscripts in preparation. EDIT - December, 2020: the published paper is available with open access on the Society of Exploration Geophysicists library, and the data is available with an accompanying Jupyter Notebook):

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from pylab import imread, imshow, gray, mean
import matplotlib.colors as cl
%matplotlib inline 

data = np.loadtxt('data.txt')
plt.imshow(data, cmap='cubehelix')
plt.show()

gets me a plot of the data:

Then I calculate the total horizontal gradient and normalize it to use for shading:

dx,dy = np.gradient(data, 1, 1)
tdx=np.sqrt(dx*dx + dy*dy)
tdx_n=(tdx-tdx.min())/(tdx.max()-tdx.min())
tdx_n=1-tdx_n

which looks as I expected:

plt.imshow(tdx_n[4:-3,4:-3],  cmap='bone')
plt.show()

To create the shading effect I thought I would get the colour from the plot of the data, then replace the opacity with the gradient so as to have dark shading proportional to the gradient, like this:

img_array = plt.get_cmap('cubehelix')(data[4:-3,4:-3])
img_array[..., 3] = (tdx_n[4:-3,4:-3]) 
plt.imshow(img_array)
plt.show()

But the result is not what I expected:

This is what I am looking for (created in Matlab, colormap is different):

Any suggestion as to how I may modify my code?

UPDATED

With Ran Novitsky's method, using the code suggested in the answer by titusjan, I get this result:

which gives the effect I was looking for. In terms of shading though I do like titusjan's own suggestion of using HSV, which gives this result: .

However, I could not get the colormap to be cubehelix, even though I called for it:

from matplotlib.colors import rgb_to_hsv, hsv_to_rgb
hsv = rgb_to_hsv(img_array[:, :, :3])
hsv[:, :, 2] = tdx_n
rgb = hsv_to_rgb(hsv)
plt.imshow(rgb[4:-3,4:-3], cmap='cubehelix')
plt.show()

First of all, Matplotlib includes a hill shading implementation. This calculates the intensity by comparing the gradient with a light source at a certain angle. So it's not exactly what you are implementing, but close and may even give better results.

Ran Novitsky has made another hill shading implementation that differs from Matplotlib in the way how the color and intensity values are merged. I can't tell which is better but it's worth a look.

Perhaps the best way of combining color and intensity would be to use gouraud shading, which is used in 3D computer graphics. My own approach, which I have implemented in the past, was to put the intensity in the value layer of the HSV color of the image.

I don't think I agree with your approach of placing the intensity (tdx_n in your case) in the alpha layer of the image. This means that where the gradient is low the image will be transparent and you will see data that was plotted earlier. I think that's what's happening in your screen shot.

Furthermore I think you need to normalize your data before you pass it through the cmap, just as you normalize your intensity:

data_n=(data-data.min())/(data.max()-data.min())
img_array = plt.get_cmap('cubehelix')(data_n)

We then can use the approach of Ran Novitsky to merge the color with the intensity:

rgb = img_array[:, :, :3]

# form an rgb eqvivalent of intensity
d = tdx_n.repeat(3).reshape(rgb.shape)

# simulate illumination based on pegtop algorithm.
rgb = 2 * d * rgb + (rgb ** 2) * (1 - 2 * d)

plt.imshow(rgb[4:-3,4:-3])
plt.show()

Or you can follow my past approach and put the intensity in the value layer of the HSV triplet.

from matplotlib.colors import rgb_to_hsv, hsv_to_rgb

hsv = rgb_to_hsv(img_array[:, :, :3])
hsv[:, :, 2] = tdx_n
rgb = hsv_to_rgb(hsv)
plt.imshow(rgb[4:-3,4:-3])
plt.show()

Edit 2015-05-23:

Your question has prompted me to finish my hill shading implementation that I started a year ago. I've put it on Github here.

It uses a blending mechanism that is similar to Gouraud shading, which is used in 3D computer graphics. It's labeled RGB blending below. I think this is the best blending algorithm, HSV blending gives erroneous results when the color is close to black (note the blue color in the center of the HSV image, which is not present in the un-shaded data).

RGB blending is also the simplest algorithm, it just multiplies the intensity with the RGB triplet (it adds an extra dimension of length 1 to allow broadcasting in the multiplication).

rgb = img_array[:, :, :3]
tdx_n_exp = np.expand_dims(tdx_n, axis=2) 
result = tdx_n_exp * rgb
plt.imshow(result[4:-3,4:-3])