What you describe is essentially how "standard" hill shading works:
- we have a height map (our DEM data) where each pixel encodes an elevation,
- we can turn that into a normal map by assigning each pixel the 3D vector
{a-b, c-d, 2}, where a and b are the elevations at (x-1,y) and (x+1,y) respectively, and c and d are the elevations at (x,y-1) and (x,y+1) respectively,
- the screen can be considered a camera that's pointing straight down onto that height map, and
- we have a light source that we model as uniform vector field, meaning every pixel gets "lit" by the same vector.
We then perform "pixel shading", where we reflect the illumination vector over each pixel's normal, and then look at the resulting z value to see how much of the light got sent straight up, into the screen.
(width, height, pixels) = extracted from GeoTIFF
hillShaded = uint32[pixels.length]
// Standard NW light source
light = normalize({ x: -1, y: -1, z: 0.4 })
// helper function to get sensible elevation values
get_elevation = (x, y) -> {
x = constrain(x, 0, width-1)
y = constrain(y, 0, height-1)
value = pixels[x + y*width]
// assuming elevation in meters, we know the nonsense values for Earth.
if (value < -499 || value > 9000) return 0
return value + 500;
}
// "pixel shade" each elevation point in our DEM data
for (let x = 0; x < width; x++) {
for (let y = 0; y < height; y++) {
// construct the normal at (x,y):
x1 = get_elevation(x-1, y)
x2 = get_elevation(x+1, y)
y1 = get_elevation(x, y-1)
y2 = get_elevation(x, y+1)
n = normalize({ x: x1-x2, y: y1-y2, z: 2 })
// compute the amount of reflected light from our light source:
reflection = reflect(light, n);
// and then turn that into a 32 bit rgba value
i = constrain_map(reflection.z, 0, 1, 0, 255)
a = i == 0 ? 0 : 255
hillShade[x + y*width] = (i<<24) + (i<<16) + (i<<8) + a;
}
}
Where constrain_map maps a value but caps it to the specified interval, and reflect is the standard vector reflection function taking v1 and v2, and yielding 2(v1·v2)/(v2·v2) * v2.
If we apply the above code to some DEM data such as:
we get:
Of course, this basic approach has rather a big problem with noise, and using a single light source tends to yield rather "stark" looking terrain the lower you set its elevation:
So to get better looking shading you generally want to smooth out small normals, and get "diffuse" lighting by using several light sources and then use a weighted sum of the resulting reflections.
And of course we can take it a few steps further and use the original height map to create an isoband stack (with bands every 100' for example) using marching squares, and then false-color those as an underlay for our hillshading, so that we end up with a composite that actually looks good:
