'Calculates the illumination intensity for a surface using the defined
azimuth and elevation for the light source.
This computes the normal vectors for the surface, and then passes them
on to `shade_normals`
Parameters
elevation : array-like
A 2d array (or equivalent) of the height values used to generate an
illumination map
vert_exag : number, optional
The amount to exaggerate the elevation values by when calculating
illumination. This can be used either to correct for differences in
units between the x-y coordinate system and the elevation
coordinate system (e.g. decimal degrees vs meters) or to exaggerate
or de-emphasize topographic effects.
dx : number, optional
The x-spacing (columns) of the input *elevation* grid.
dy : number, optional
The y-spacing (rows) of the input *elevation* grid.
fraction : number, optional
Increases or decreases the contrast of the hillshade.  Values
greater than one will cause intermediate values to move closer to
full illumination or shadow (and clipping any values that move
beyond 0 or 1). Note that this is not visually or mathematically
the same as vertical exaggeration.
Returns
intensity : ndarray
A 2d array of illumination values between 0-1, where 0 is
completely in shadow and 1 is completely illuminated.'

def hillshade(self, elevation, vert_exag=1, dx=1, dy=1, fraction=1.0):

    dy = (- dy)
    (e_dy, e_dx) = np.gradient((vert_exag * elevation), dy, dx)
    normal = np.empty((elevation.shape + (3,))).view(type(elevation))
    normal[..., 0] = (- e_dx)
    normal[..., 1] = (- e_dy)
    normal[..., 2] = 1
    normal /= _vector_magnitude(normal)
    return self.shade_normals(normal, fraction)
