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)