Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lightsource enhancements #3291

Merged
merged 19 commits into from Sep 29, 2014

Conversation

joferkington
Copy link
Contributor

Currently, colors.LightSource preforms "hillshading" through methods that appear visually unrealistic in some cases (particularly topographic surfaces). This pull request:

  1. Refactors LightSource to allow the hillshade method to be called independently so that illumination maps can be produced.
  2. Adds support for common parameters such as vertical exaggeration and x/y spacing so that topographically accurate illumination can be calculated. (Or so that the effect of topography can be exaggerated or subdued in the final result.)
  3. Adds support for multiple blend modes to combine the illumination map with the rgb data. blend_mode='overlay' and blend_mode='soft' are both more visually appealing for complex surfaces than the default mode. Effectively, the other blend modes make the surface appear less "shiny".

All of the changes made should be fully backwards-compatible and none of the defaults have been changed. This also adds a 170Kb example data file of some SRTM elevation data (public domain) to the sample data.

As a quick example of the effect of some of the changes:

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

dem = np.load(get_sample_data('jacksboro_fault_dem.npz'))
z = dem['elevation']

# Shade from the northwest, with the sun 45 degrees from horizontal
ls = LightSource(azdeg=315, altdeg=45)
cmap = plt.cm.gist_earth

fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(8, 9))
plt.setp(axes.flat, xticks=[], yticks=[])

# Vary vertical exaggeration and blend mode and plot all combinations
for col, ve in zip(axes.T, [0.01, 0.1, 10]):
    # Show the hillshade intensity image for each vert. exag.
    col[0].imshow(ls.hillshade(z, vert_exag=ve), cmap='gray')
    for ax, mode in zip(col[1:], ['hsv', 'overlay', 'soft']):
        rgb = ls.shade(z, cmap=cmap, vert_exag=ve, blend_mode=mode)
        ax.imshow(rgb)

# Label rows and columns
for ax, ve in zip(axes[0], [0.01, 0.1, 10]):
    ax.set_title('{}'.format(ve), size=18)
for ax, mode in zip(axes[:,0], ['Hillshade', 'hsv', 'overlay', 'soft']):
    ax.set_ylabel(mode, size=18)

# Group labels...
axes[0,1].annotate('Vertical Exaggeration', (0.5, 1), xytext=(0, 30), 
                   textcoords='offset points', xycoords='axes fraction', 
                   ha='center', va='bottom', size=20)
axes[2,0].annotate('Blend Mode', (0, 0.5), xytext=(-30, 0), 
                   textcoords='offset points', xycoords='axes fraction', 
                   ha='right', va='center', size=20, rotation=90)
fig.subplots_adjust(bottom=0.05, right=0.95)

plt.show()

figure_3

@jenshnielsen
Copy link
Member

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.

@joferkington
Copy link
Contributor Author

@jenshnielsen - Thanks, good points!

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

@WeatherGod
Copy link
Member

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)
.

@tacaswell tacaswell added this to the v1.5.x milestone Jul 23, 2014
@joferkington
Copy link
Contributor Author

@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.

ls = mcolors.LightSource(315, 45)
elev = dem['elevation']
cmap = cm.gist_earth

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test raises a unclosed file warning on python3. Could you add a dem.close() statement here when you are done with the file. We are trying to reduce the number of warnings in the test suite.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it wouldn't be dem.close(), because it is a dictionary of numpy arrays.
It is the cbook.get_sample_data() that returns a file object. So
numpy.load() needs to do the closing, I think.

On Sat, Jul 26, 2014 at 5:21 AM, Jens H Nielsen notifications@github.com
wrote:

In lib/matplotlib/tests/test_colors.py:

+@image_comparison(baseline_images=['light_source_shading_topo'],

  •              extensions=['png'])
    
    +def test_light_source_topo_surface():
  • """Shades a DEM using different v.e.'s and blend modes."""
  • dem = np.load(cbook.get_sample_data('jacksboro_fault_dem.npz'))
  • Get the true cellsize in meters for accurate vertical exaggeration

  • Convert from decimal degrees to meters

  • dx, dy = dem['dx'], dem['dy']
  • dx = 111320.0 * dx * np.cos(dem['ymin'])
  • dy = 111320.0 * dy
  • ls = mcolors.LightSource(315, 45)
  • elev = dem['elevation']
  • cmap = cm.gist_earth

The test raises a unclosed file warning on python3. Could you add a
dem.close() statement here when you are done with the file. We are trying
to reduce the number of warnings in the test suite.


Reply to this email directly or view it on GitHub
https://github.com/matplotlib/matplotlib/pull/3291/files#r15432988.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jenshnielsen - Thanks for the catch!

@WeatherGod - Either one works, actually. dem is an NpzFile object. Closing it closes the file. I thought the same thing until I dug a bit deeper, though.

Either way, I went ahead and changed the test and example to use a context manager. (I didn't realize until today that np.load supported that, but it does for .npz files, and apparently has for the past few versions.) See: 896bd37

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, looks like that actually doesn't close the underlying file (or, more precisely, the file returned by get_sample_data is still open even though the NpzFile object is closed). Still getting the warnings. I pushed things a bit too soon. See: bdbc7df for the correct fix.

@joferkington
Copy link
Contributor Author

@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.

@Patrick-Cole
Copy link

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_1
Figure showing magdata using standard hillshading.

col_clip_pc
Figure showing colors clipped.

  1. 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?

  2. 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?

rad
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.

@jenshnielsen
Copy link
Member

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

@joferkington
Copy link
Contributor Author

@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()

figure_1

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()

figure_1

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()

figure_1

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.

@jenshnielsen
Copy link
Member

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

@Patrick-Cole
Copy link

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

@jenshnielsen
Copy link
Member

@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.

@joferkington
Copy link
Contributor Author

@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.

@jenshnielsen
Copy link
Member

@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.

@jenshnielsen
Copy link
Member

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

@joferkington
Copy link
Contributor Author

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.

@jenshnielsen
Copy link
Member

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?

@tacaswell
Copy link
Member

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

@jenshnielsen
Copy link
Member

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.

@joferkington
Copy link
Contributor Author

@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.

@jenshnielsen
Copy link
Member

@joferkington thanks missed that due all the other commits

@joferkington
Copy link
Contributor Author

@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!)

@jenshnielsen
Copy link
Member

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

jenshnielsen added a commit that referenced this pull request Sep 29, 2014
@jenshnielsen jenshnielsen merged commit 34ed45a into matplotlib:master Sep 29, 2014
@joferkington
Copy link
Contributor Author

@jenshnielsen - Thanks!!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants