I usually post-process a copy of the raw hillshade in some external graphics editing program that I can directly over-write (save, instead of save-as)
Photoshop and GIMP are two that come too mind.
The process in a nutshell I follow:
Import GeoTiff'd Hillshade -> Tone -> Median -> Reduce Noise -> Gaussian Blur
- Your geotiff should be accompanied with a worldfile (like .tfw) for portability between the GIS and graphics program you use)
Use the tone tool to redistribute your min/max values of the greyscale hillshade. In GIMP you use the 'Levels' tool.
I often use a very small range with a min of 5% grey and 25% grey. Rarely do I use full white to black greyscale ranges.
Median is used to generalize. This is the part where can easily misrepresent your hillshade. In GIMP the equivalent tool is 'Despeckle'.
The median/despeckle tool will reduce the sharp edges and make then rounded, like mountian tops and valleys. Thus often making it appear as though some valleys are filled-in and some mountain tops are flatter than they appear. Use this tool lightly, or not at all.
Reduce noise is used to remove any strange artifacts that median couldn't - often I don't even use reduce noise. Almost like a cheat way of de-striping your DEM if it suffers from 'stepping' - a result of manual profiling.
Gaussian blur will introduce the final smoothing of the hillshade and is often the last thing I do (this will also help remove any introduced artifacts).
The other option, is to run the above process over an HDR hillshade.
Create hillshades from various aspects like 360 (or 0 deg), 15, 30, 45.
Those are the four angles I might use for my local region because that is the general direction the sun comes from during the morning to mid-day. HDR the four images using either GIMP or photoshop then use whatever generalizing filters you want like the aforementioned ones.
Experiment and see what works best for you.
What's good about working on a GeoTiff is that none of the georeferencing is lost as the pixel size and extent are not altered in the process.
There are lots of other filters/overlays to use and play with. The ones I've mentioned are the ones I use most frequently. Also remember that the less filters you have to apply the better so that the terrain doesn't become dramatically misrepresented.
When I apply the filters, they are done so in small doses. For instance, I find myself never applying a gaussian blur over 5%.
Useful resources:
http://www.shadedrelief.com/bumping/bumping.html
http://www.shadedrelief.com/shading/Swiss.html
Vince's comments led me to look again at Shapely in Python, and it can do it. In short, the method I used was:
- get the bounding box of the page/object you want
- calculate the diagonal of this box, as this will be the minimum length your hatch lines need to be
- draw a square array of horizontal lines centred on the bounding box's centre, each spaced suitably apart
- rotate these lines to the desired angle
- clip them to the bounding box.
Somewhat crude example code:
#!/usr/bin/python
# -*- coding: utf-8 -*-
# shapely_hatch - simple hatching function demo
# produces WKT of a 45° hatched crescent to stdout
# scruss — 2014-04-13 — WTFPL (srsly)
from shapely.geometry import box, MultiLineString, Point
from shapely.affinity import rotate
from shapely import speedups
from math import sqrt
# enable Shapely speedups, if possible
if speedups.available:
speedups.enable()
def hatchbox(rect, angle, spacing):
"""
returns a Shapely geometry (MULTILINESTRING, or more rarely,
GEOMETRYCOLLECTION) for a simple hatched rectangle.
args:
rect - a Shapely geometry for the outer boundary of the hatch
Likely most useful if it really is a rectangle
angle - angle of hatch lines, conventional anticlockwise -ve
spacing - spacing between hatch lines
GEOMETRYCOLLECTION case occurs when a hatch line intersects with
the corner of the clipping rectangle, which produces a point
along with the usual lines.
"""
(llx, lly, urx, ury) = rect.bounds
centre_x = (urx + llx) / 2
centre_y = (ury + lly) / 2
diagonal_length = sqrt((urx - llx) ** 2 + (ury - lly) ** 2)
number_of_lines = 2 + int(diagonal_length / spacing)
hatch_length = spacing * (number_of_lines - 1)
# build a square (of side hatch_length) horizontal lines
# centred on centroid of the bounding box, 'spacing' units apart
coords = []
for i in range(number_of_lines):
# alternate lines l2r and r2l to keep HP-7470A plotter happy ☺
if i % 2:
coords.extend([((centre_x - hatch_length / 2, centre_y
- hatch_length / 2 + i * spacing), (centre_x
+ hatch_length / 2, centre_y - hatch_length
/ 2 + i * spacing))])
else:
coords.extend([((centre_x + hatch_length / 2, centre_y
- hatch_length / 2 + i * spacing), (centre_x
- hatch_length / 2, centre_y - hatch_length
/ 2 + i * spacing))])
# turn array into Shapely object
lines = MultiLineString(coords)
# Rotate by angle around box centre
lines = rotate(lines, angle, origin='centroid', use_radians=False)
# return clipped array
return rect.intersection(lines)
# pipe-separated output; can be read by QGIS
print 'ID| WKT'
page = box(1000, 1000, 6000, 6000)
hatching = hatchbox(page, 45, 50)
circle = Point(2500, 2500).buffer(1000)
circle1 = Point(2000, 2500).buffer(700)
crescent = circle.difference(circle1)
crescent_hatch = crescent.intersection(hatching)
print '1|', crescent_hatch
The output, when graphed, looks like this:
Best Answer
It is important to remember that when computing hillshading, you need to have an illumination source. Using the sun as an illumination source may mean that a cell is shaded at noon, when the sun is directly overhead, but not at 4:00 p.m. Without an illumination point, your example seems more like a color coded slope map.
ESRI calculates illumination of each cell relative to its neighboring cells and has better explanations and examples of their algorithm than I can offer:
http://edndoc.esri.com/arcobjects/9.2/NET/shared/geoprocessing/spatial_analyst_tools/how_hillshade_works.htm