I have not done a full review but this looks like great stuff.

Could you have a look at adding some more tests for the code. The only present test that I know of test_light_source_shading_color_range
in test_colors.py This test is somewhat broken since it calculates the shading of a flat surface, which results in a divide by zero error when rescaling the intensity see #3248

Perhaps you could add some more real tests? Something like your example above but shading a more simplified surface, but one with a curvature perhaps just -x^2 - y^2. You can take inspiration in the existing test_light_source_shading_color_range if you like.

I think it would be good to make sure that hillshade cannot do a divide by zero by checking that
imax > imin before doing the normalisation step. If imax==imin we should either raise an exception since there is little point in hill shading a flat surface or normalize in a different way. This will of cause require updating the existing test output.

The present test failures are pep8 style fixes which should not be to difficult to fix.

@jenshnielsen - Thanks, good points!

I made the PEP8 fixes and I'll add some tests tonight. (Should have included both of those originally!)

This is good to see. There are some places in mplot3d where LightSource is
used (namely, in plot_surface()). Usually, one uses the default
construction, but can pass in your own. Maybe a new example demonstrating
that might be nice?

On Wed, Jul 23, 2014 at 8:25 AM, joferkington notifications@github.com
wrote:

@jenshnielsen https://github.com/jenshnielsen - Thanks, good points!

I made the PEP8 fixes and I'll add some tests tonight. (Should have
included both of those originally!)

—
Reply to this email directly or view it on GitHub
#3291 (comment)
.

@WeatherGod - Good point! I added an example in 31a4239. There's an undocumented lightsource kwarg to plot_surface, but for this purpose it's easier to just calculate the rgb facecolors and pass them in.

@Patrick-Cole - I just saw your notebook on hillshading, and if you have a chance, I'd appreciate your input on this.

I'm using a completely different algorithm than yours (I hadn't seen Horn (1981) before), but it's basically accomplishing the same thing. This is a more "traditional" approach (i.e. intensity = scalar_dot_product(illumination_vector, terrain_normal_vector)) than Blinn shading. The "knobs" you can turn here are the vertical exaggeration and blending modes.

(Side note: a vertical exaggeration of 1 is only accurate if the correct dx and dy are supplied. In the figure at the top of the page, V.E.=0.01 is actually the closest to a true V.E. of 1 because the cellsize is about 90 meters (1/90 ~ 0.01). To do this properly, you'd specify dx=90, dy=90 to ls.shade and then use realistic vertical exaggeration numbers.)

Your scipy talk was a large part of the impetus for this pull request. I'd written a lot of this years ago, but never thought to contribute it back upstream, as I figured I was the only one using it.

Hi,

I apologize for the late reply. I have been playing with your hillshading and I quite like the results. We tend to use this as a tool for interpreting data sets such as magnetic data, since it is useful at highlighting small structures in a directional way. I do have some comments (which may or may not lead to possible enhancements).

1. Magnetic data unfortunately has large outliers. The data set I tested with had a mean of 28,661, a minimum value of 21,156 and a maximum of 46,286. Unfortunately the standard deviation is only 243. This results in a colorless image. I get around this by applying some sort of clip to the colors of the shading only. If the clip is applied to the shading, it results in clipping of the peaks and troughs, which is not always desirable. Is there a way to do this here?

Figure showing magdata using standard hillshading.

Figure showing colors clipped.

2. When I tried to have a colorbar, it gave me a range of 0. to 1. I think it is useful to see the actual range applied to the color portion of the hillshading. Is this possible?

3. On some occasions, an interpreter may which the shading to be elevations, but the colors to be another dataset, say radiometrics. Is the mixing of datasets something that can be considered?

Figure showing radiometric total counts overlain on elevations. You can see the alluvium draining down the valleys in the north.

I hope this helps. If it is out of the scope of this addition to matplotlib, then no worries.

@Patrick-Cole Is it possible to use the vmin and vmax arguments to imshow to get rid of the outliers?

@jenshnielsen - Not directly, because imshow just gets passed a pre-rendered RGBA array in this case. Instead, you'd pass in nom=Normalize(vmin, vmax) to ls.shade. It's a bit clunky, but I didn't want to change the API to LightSource any more than I had to. It's a good idea to add vmin and vmax kwargs to ls.shade, though...

@Patrick-Cole - Thanks for testing things! Actually, you've revealed my ulterior motive in sending in this pull request.

What you're asking for is all entirely reasonable, but is rather clunky to do as things are currently structured. In the long run, I'd like to put together an MEP (enhancement proposal) for a shadedrelief axes method (and therefore pyplot function). To enable the colormap, etc, to work correctly, it needs to return a new subclass of AxesImage. This would also allow more efficient updating and separation of the data being shaded and the data being colored.

A new function like this would provide a simpler interface to the shading functionality. However, I wanted to get this pull request accepted before I proposed larger changes. (Also, like I mentioned, I think a new axes method arguably needs a "full-blown" MEP, rather than just a pull request.)

At any rate, to answer your questions, here's how you'd do what you're asking with the current (i.e. this pull request) state of things:

Avoid outliers

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource, Normalize

y, x = np.mgrid[-4:2:200j, -4:2:200j]
z = 10 * np.cos(x**2 + y**2)

# Add some outliers...
z[10, 5] = 2000
z[3, 15] = -9000

ls = LightSource(315, 45)
rgb = ls.shade(z, plt.cm.copper, norm=Normalize(-10, 10)) # Note the norm!!

fig, ax = plt.subplots()
ax.imshow(rgb)
plt.show()

Shade different data

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource, Normalize

y, x = np.mgrid[-4:2:200j, -4:2:200j]
z1 = np.sin(x**2) # Data to hillshade
z2 = np.cos(x**2 + y**2) # Data to color

norm=Normalize(z2.min(), z2.max())
cmap = plt.cm.jet

ls = LightSource(315, 45)
rgb = ls.shade_rgb(cmap(norm(z2)), z1)

fig, ax = plt.subplots()
ax.imshow(rgb)
plt.show()

Display a correct colorbar

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource

y, x = np.mgrid[-4:2:200j, -4:2:200j]
z = 10 * np.cos(x**2 + y**2)

cmap = plt.cm.copper
ls = LightSource(315, 45)
rgb = ls.shade(z, cmap)

fig, ax = plt.subplots()
ax.imshow(rgb)

# Use a proxy artist for the colorbar...
im = ax.imshow(z, cmap=cmap)
im.remove()
fig.colorbar(im)

plt.show()

At any rate, it's all possible, but often clunkier than it should be. A separate axes method and a separate artist for shaded plots would make this a lot simpler, assuming there's enough interest.

@joferkington I think it would be great to add your examples to the documentation as part of the pull request.

Thanks for the examples. I agree with them being added to the documentation. Its definitely looking good so far!

@joferkington Any progress on adding the examples to the documentation?

I would like to add some tests in the style of the original test that you removed (i.e.
array comparison rather than image comparison) since these are much more isolated.
I am happy to add these as a pull request to your pull request, just let me know.

@jenshnielsen -
Sorry for the delay on this. Work went from being busy to being absolutely frantic for the past month or so. The array comparison I removed didn't actually test anything (It was a flat plane), so I added a different plane comparison test in its place. I incorporated a couple of different array comparison tests in abf8b87. Let me know what you think.

@joferkington Great I will have a more careful look later. I agree on removing the original array tests since they were basically broken. So adding new array tests seem sensible to me. The main point is that the array tests will test this piece of code more directly independent from any other code in matplotlib.

BTW: I think that this is a significant improvement so you may want to add a small entry to the whats new folder too.

Gah! I managed to mess up my fork and branch somehow... Just a bit, let me see if I can fix this. In the meantime, ignore all the commits at the end.

The tests in abf8b87 fail due to pep8 issues in the included np arrays. I don't think that fixing this is a good idea since the arrays will be less readable if pep8ed. Personally I think that the best solution is to store the arrays in separate npz files rather than in the source code. @tacaswell @mdboom any thoughts on that?

I don't think we should be pep8ing the tests.

I think that I agree with that. We can remove the pep8 on that file but in principle I think that it is a good idea to separate tests from reference data. Also when the reference is not images.

@jenshnielsen - On a side note, I fixed the pep8 issues on the printed arrays in fe474e7. The result is still quite readable, though I agree I might be better not to have the reference data inline.

@joferkington thanks missed that due all the other commits

@jenshnielsen - Yeah, sorry about those. I believe I have all the duplicate commits fixed now. (On a side note, rebase can do some odd things under certain circumstances!)

This looks good to me. A great addition with good documentation and testing so I will merge this. Thanks for all your work.

@jenshnielsen - Thanks!!