Wednesday, April 21, 2010

Using hillshade image as intensity (improved matplotlib shade)

Matplotlib module enables hillshade method (shade) using a LightSource class (v 0.99).
The problem is it uses the data itself as intensity and data. It is very useful for viewing a DEM but sometimes you would like the DEM as intensity underlying some other data. Another problem is that the shade method is producing a very light colored image sometimes even white where intensity is high.
I used as an example a DEM derived from SRTM v4 data acquired at the International  Centre for Tropical  Agriculture (CIAT - http://srtm.csi.cgiar.org) the hillshade production was made using LightSource class with azimuth of 165 deg. and altitude of 45 deg.)
DEM - gist_earth color scheme Hill-shade (azdeg-165,altdeg-45)
matplotlib shade method My shade method
The difference in the shading colors derived from the method used to produce it. While the matplotlib method uses "hard light" method I use a "soft light" method. the matplotlib is converting the RGB colors to HSV and then calculate the new saturation and value according to the intensity. I use a formula based on the description of ImageMagick's pegtop_light.which is much faster as it is a single formula. Another advantage is the option to use a separate layer as the intensity and another as the data used for colors.
The modified functions are hillshade and set_shade as follows:
#!/bin/env python
from pylab import *
def set_shade(a,intensity=None,cmap=cm.jet,scale=10.0,azdeg=165.0,altdeg=45.0):
''' sets shading for data array based on intensity layer
  or the data's value itself.
inputs:
  a - a 2-d array or masked array
  intensity - a 2-d array of same size as a (no chack on that)
                    representing the intensity layer. if none is given
                    the data itself is used after getting the hillshade values
                    see hillshade for more details.
  cmap - a colormap (e.g matplotlib.colors.LinearSegmentedColormap
              instance)
  scale,azdeg,altdeg - parameters for hilshade function see there for
              more details
output:
  rgb - an rgb set of the Pegtop soft light composition of the data and 
           intensity can be used as input for imshow()
based on ImageMagick's Pegtop_light:
http://www.imagemagick.org/Usage/compose/#pegtoplight'''
  if intensity is None:
# hilshading the data
    intensity = hillshade(a,scale=10.0,azdeg=165.0,altdeg=45.0)
  else:
# or normalize the intensity
    intensity = (intensity - intensity.min())/(intensity.max() - intensity.min())
# get rgb of normalized data based on cmap
  rgb = cmap((a-a.min())/float(a.max()-a.min()))[:,:,:3]
# form an rgb eqvivalent of intensity
  d = intensity.repeat(3).reshape(rgb.shape)
# simulate illumination based on pegtop algorithm.
  rgb = 2*d*rgb+(rgb**2)*(1-2*d)
  return rgb

def hillshade(data,scale=10.0,azdeg=165.0,altdeg=45.0):
  ''' convert data to hillshade based on matplotlib.colors.LightSource class.
    input:
         data - a 2-d array of data
         scale - scaling value of the data. higher number = lower gradient
         azdeg - where the light comes from: 0 south ; 90 east ; 180 north ;
                      270 west
         altdeg - where the light comes from: 0 horison ; 90 zenith
    output: a 2-d array of normalized hilshade
'''
  # convert alt, az to radians
  az = azdeg*pi/180.0
  alt = altdeg*pi/180.0
  # gradient in x and y directions
  dx, dy = gradient(data/float(scale))
  slope = 0.5*pi - arctan(hypot(dx, dy))
  aspect = arctan2(dx, dy)
  intensity = sin(alt)*sin(slope) + cos(alt)*cos(slope)*cos(-az - aspect - 0.5*pi)
  intensity = (intensity - intensity.min())/(intensity.max() - intensity.min())
  return intensity




Example of use:
One can save the code to a file named say: shading.py
now say we have a 4 byte float DEM data in a 560 lines 420 samples binary file. in a python code:

from pylab import *
from shading import set_shade
from shading import hillshade
dem = fromfile('DEM.dem',dtype=float32).reshape(560,420)
rgb = set_shade(dem,cmap=cm.gist_earth)
imshow(rgb)

will produce the "my shade method" image as above.

say we have a data to be plot using the DEM data as intensity:

replace the line before last with:
rgb = set_shade(data,intensity=hillshade(dem),cmap=cm.gist_earth)

5 comments:

  1. Thanks for a great post. It helped me a lot. :)

    ReplyDelete
  2. Thanks for sharing. Really smart output. Have you considered submitting it to mpl?

    ReplyDelete
    Replies
    1. I don't know how to do it. I wrote about it to one of the developers but I'm not sure he have done anything with it.

      Delete
  3. Thanks for sharing Ran. There is a question. The normal plot of python sets "nan" as white color. However, when I try your code to add intensity, "nan" is set as the same color of the lowest value. How do you deal with this?

    ReplyDelete
    Replies
    1. Hi Wenliang Zhao,
      You can try to use masked arrays on the output using numpy.ma.masked_where function (see: http://matplotlib.org/examples/pylab_examples/image_masked.html)
      or use matplotlib cmap.set_bad function for the colormap (see: http://stackoverflow.com/questions/2578752/how-can-i-plot-nan-values-as-a-special-color-with-imshow-in-matplotlib)

      Delete

Please Comment this Post or send me an Email