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

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)

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

ReplyDeleteThanks for sharing. Really smart output. Have you considered submitting it to mpl?

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

DeleteThanks 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?

ReplyDeleteHi Wenliang Zhao,

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