Making shaded relief directly from DEMs projected in degrees?

In a presentation I gave at GeoIgnite 2021, I was explaining the process of using QGIS to make shaded relief from DEMs, more or less as detailed in this post. There wasn’t time to explain the process of re-projection, so I rather boldly asserted that you don’t need to re-project your DEMs before making shaded relief. Instead, I said, you can just use the magic scaling number 111120 to convert from degrees to metres.

Really? Is this true?

Well, yes, in many cases it works pretty well. In QGIS you pull in an unreprojected DEM (its native CRS using degrees as units) but set the Project CRS to some reasonable metric projection, such as UTM. Then, opening the GDAL hillshade tool from the Toolbox, you use the DEM as your Input layer, set scale to 111120, and enter the parameters of your choice for Z factor, Azimuth of the light and Altitude of the light.

And, look: shaded relief:

This seems so much simpler than having to first re-project your DEM to square cells measured in metres. But there are some caveats. This technique fools with the terrain a bit, and could potentially distort how it looks. If you understand how it distorts it, you’ll be able to evaluate if you want to use it in your map.

The magic number

The magic number 111120 comes from the GDAL documentation, specifically the documentation for the gdaldem tool:

The gdaldem generally assumes that x, y and z units are identical. If x (east-west) and y (north-south) units are identical, but z (elevation) units are different, the scale (-s) option can be used to set the ratio of vertical units to horizontal. For LatLong projections near the equator, where units of latitude and units of longitude are similar, elevation (z) units can be converted to be compatible by using scale=370400 (if elevation is in feet) or scale=111120 (if elevation is in meters). For locations not near the equator, it would be best to reproject your grid using gdalwarp before using gdaldem.
https://gdal.org/programs/gdaldem.html

Gdaldem is the tool that I recommend using for creating hillshades in QGIS, and there in the docs they say it: unless your locale is near the equator, it’s best to re-project your grid before using gdaldem.

The obvious problem with entering a scale factor of 111,120 metres per degree is that as you move away from the equator, lines of longitude become closer to one another. For example, in Canada at 54° north, one degree of latitude is close to 111,0000 metres, but a degree of longitude is only about 65,000 metres. In other words, a DEM cell is not square on the ground. So when you enter 111,120 as the scaling factor to convert degrees to metres, you will get a hillshade calculated on the basis of a square cell that has been created by stretching a rectangular cell in an east–west direction.

For example, with a 1″ DEM, each cell is 0.0002777777° (that’s 1 ÷ 3600) on a side. If we scale both axes at 111,120, the hillshading algorithm thinks it has a square cell 30 metres on a side. However the reality at 54° north is that on the ground the cell is about 30 m north–south and 18 metres east–west.

How bad is it?

You’d expect that because the landscape is stretched east–west, the slope would be under-calculated for cells that slope east–west. Basically, these cells would appear to be on a slope less steep than it really is.

And slope is important in the hillshading algorithm: it plays a role in determining how bright a hillshade pixel is. Pixels whose aspect and slope indicate that they are orthogonal to the sun direction are the brightest white, while those facing away from the sun direction are darkest black.

So if, for example, you had a 45° slope facing west, and the sun in the west at an elevation of 45 °, you’d expect pixels on the slope to have the highest possible hillshade value, 255. But if you calculated the hillshade from a degree-DEM using a scale of 111,120, that slope would be under-estimated; maybe it would appear to be 35°. As a result it would not get assigned a 255, but something less, maybe 240. There might be a nearby area, also facing west, whose true slope was 55°, but it now gets assigned the brightest value because it appears to be 45°.

In general then, for east or west-facing surfaces, we’d expect the hillshade built from the un-re-projected degree-DEM to distribute its brightest and darkest patches differently, particularly when the sun is due east or west.

That’s just what you see.

However, when the sun is not due east or west it’s hard to see this difference in slope calculation.

You can also experiment with other, lower scale factors, such as 90,000, 80,000 or 70,000. These, in effect, cause slopes to be less under-estimated on east–west aspects, but over-estimated on north–south aspects. In other words you are distributing the error around.

This suggests that in many cases, making a hillshade directly from a DEM projected in degrees is fine.

Other questions

A second distortion arising from creating a hillshade using a degree-DEM and a scale factor of 111,120 concerns the sun’s azimuth. If you are mapping in a local UTM projection, you know that UTM grid north is usually different from true north. This was visible in the figure above with the contour lines: the cell that was “square” in degrees turns out to be rectangular in a UTM projection, and also a bit tilted to the side.

With UTM projections, this difference is not more than 3°. At worst, then, you might ask for a sun azimuth of 315° and instead get one of 312°. I doubt whether you can notice this difference.

The last problem we might wonder about is whether shadows calculated by the hillshading algorithm are distorted. Are they foreshortened in an east–west direction when the hillshade is re-projected into a metric projection? The answer to this may be surprising: the GDAL hillshading algorithm does not calculate shadows. It calculates shadowed surfaces, but this is purely on the basis of the slope and aspect of the cell it is looking at (and its eight neighbours).

In other words, a piece of terrain at the bottom of a deep gorge will be rendered purely on the basis of its slope and aspect, as if the sun could shine on it right through the surrounding mountains.

GDAL hillshade: no shadows on the glacier, despite big peaks on its north side

Hence, no, we do not need to worry that shadows will be shortened.

So…

All of these things considered, I come to the surprising (for me) conclusion that yes, you actually can often just hillshade a DEM that’s still projected in degrees, using the scale factor 111,120. It seems heretical, but when I can’t tell the difference in the results, I have to admit that for practical purposes it’s a good solution.