2

I'm curious how GIS systems actually generate a hill-shade from input data like a raster DEM layer. I'd like to implement my own in software (mostly for fun/learning) but I want to start with a plausible algorithm so I don't waste an incredible amount of time.

Naively, I'd think you could scan the elevation data, calculate the slope from point to point in a given direction (say left to right), then color all of either the positive or negative slopes. If the min and max slopes in the elevation tile were known, it seems like it might also be nice to scale the color based on that.

Is that a plausible way to go about the problem?

3
  • It's a bit more complicated than that. There are also a number of different implementations, even in one organization. In order to avoid having this closed as too broad or opinion-based you may need to focus on a single implementation, without looking for validation of your own algorithm.
    – Vince
    Mar 3, 2016 at 21:57
  • 1
    If you want to see some actual code have a look at gdaldem, the code can be downloaded, inspected (and even modified) here.
    – Ian Turton
    Mar 4, 2016 at 7:33
  • 1
    There is some nice .js scripting for doing Hillshades in OpenLayers, I don't understand the detail much myself but if you like .js this should give you a good intro to one method.
    – Sethinacan
    Mar 4, 2016 at 8:04

2 Answers 2

5

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

2
  • He does have an illumination source - if illuminating from left to right, then the source is directly to the left (west) at an infinite distance, based on what he described.
    – Dan
    Nov 21, 2022 at 13:47
  • Unfortunately, it looks like the "edndoc" part of the ESRI website no longer exists, and the original article does not appear to exist on the current website. Jun 20, 2023 at 15:32
0

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:

Some DEM data for a part of Vancouver Island

we get:

The same data, but hill-shaded

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:

The same data, but hill-shaded with a very low "sun"

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:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.