Shaded Relief using Skymodels, courtesy of Raster Chunk Processing

A couple of weeks ago I watched some excellent presentations from the How To Do Map Stuff day organized by Daniel Huffman. One that I particularly enjoyed was Jake Adams’s talk on building shaded relief out of hillshades. Toward the end of his talk he brought in something called a skymodel.

In this post, I’ll explain what Skymodels are, and how to get Jake’s Raster Chunk Processing software running in a linux environment so you can use the Skymodel feature to make your own versions of this unique kind of shaded relief.

Skymodels were introduced in a paper by Pat Kennelly and James Stewart in 2014. The basic idea behind the skymodel is that, in the real world, we see terrain illuminated by light coming from all parts of the sky in various amounts, not just from the sun. The usual hillshading algorithm, on the other hand, calculates what’s illuminated and what’s in shadow based entirely on light coming from a single point.

An overcast sky is one example of a skymodel: in the overcast sky, light comes mostly from above, but not from any particular direction.  In James Stewart’s SkyLum software, this is represented by the redder, “hotter” dots in the dome of the model, while we have yellow and then cooler blue dots down toward the horizon representing less light coming in from those directions.

skymodel_overcast

Contrast that with this model for what SkyLum calls a Type 6 sky, “partly cloudy, no gradation towards zenith, slight brightening towards the sun.”

Type 6 partly cloudy, no gradation towards zenith, slight brightening towards the sun

Or this, the Type 13, “CIE standard clear sky, polluted atmosphere.”

type 13 CIE cstandard lear sky, polluted atmosphere

Each of these will produce a different kind of shaded relief.

Baltoro_assortment
The Baltoro glacier, Pakistan. Left: overcast. Middle: type 6, “partly cloudy, no gradation towards zenith, slight brightening towards the sun.” Right: type 13, “CIE standard clear sky, polluted atmosphere.”

In Blender you can position multiple light sources, but SkyLum writes a CSV file (which I will call the illumination file) that defines 250 different light sources and assigns a relative weight to each. This is new level of complexity in lighting.

180,45,0.000962426
175,37,0.000929834
193,47,0.00118523
187,38,0.000949383
167,47,0.00143157
181,54,0.00168069
169,41,0.00141631
...

Each line in the illumination file above gives the azimuth, elevation, and weight for one illumination source. The weights for all 250 points will add up to 1.

So how does one use this to generate shaded relief?

Well, Jake Adams (thank you, Jake!) has written a clever piece of code called Raster Chunk Processing (RCP, hereafter), which he presents in a three-part blog post. RCP divides up large DEMs into smaller tiles (“chunks”), each of which is processed separately. All of the results are then merged back together for a final result. The point of RCP is that it allows you to work with DEMs that ordinarily would max out your RAM and cause processing to grind to a halt.

This is similar to the way I have divided large DEMS into tiles for processing by Blender, but Jake’s RCP code allows you to use this “divide-and-conquer” strategy to do a whole host of things. One of them is to build a skymodel hillshade, given a DEM and an illumination file from SkyLum.

The RCP code is written in python, which is platform independent, so although Jake gives us instructions for installing it under Windows, we can get it running under linux with just a few changes.  In this case I did this on a system running Ubuntu 18 Bionic, and since I was using this machine for mapping I already had QGIS and GDAL installed.

To begin, head over to Jake’s Git repository, and download the RCP software using the Clone or Download button. Once unzipped, this will produce a directory called rcp-master. Within it you will see, among other things, the main program (raster_chunk_processing.py), a subdirectory called SkyLum (which contains the SkyLum program) and a couple of other python files that will be important to us, like settings.py and methods.py.

RCP is written to be run under python 3, so from a Ubuntu point of view these are the packages you will need to install:

  • python3-numba
  • python3-astropy
  • python3-gdal
  • python3-skimage
  • python3-numpy

(Note that python3-scimage replaces the module Jake calls scikit-image.)

Another thing that RCP needs is a working copy of mdenoise, the software that implements Sun’s algorithm for denoising topographic data. Or rather, RCP needs to at least think it has a copy. So you have a choice: if you want to be able use RCP to denoise DEMs, you should compile yourself a copy of the mdenoise binary; there are instructions here, and it’s not too painful. If not, just use sudo to place an empty file called mdenoise in /bin.

Then use a text editor on the file settings.py to alter its single line about where the mdenoise executable is.

MDENOISE_PATH = '/bin/mdenoise'

The last piece you have to put in place is a way to run SkyLum. SkyLum only comes as a Windows binary, SkyLum.exe, so to run this you need to have wine installed. The good news is that SkyLum runs quite well under wine.

Right-click SkyLum.exe, and choose Open With… >Other Application. In the list of application choose, “winebrowser.”

winebrowser

SkyLum should open right up and show you a piece of terrain with a sky dome over it.

fresh skylum

The complete instructions for SkyLum are in the README file included with it, but I will give a summary here of what I find useful:

  • By default when you open SkyLum the sun (sun position is shown at the bottom) is at 45° elevation and azimuth 180° (due South). 0° is north, and azimuths increase clockwise, as is standard
  • Hit ? for a help screen
  • Move the sun with the arrow keys.
  • Rt-click and choose an illumination model. See Kennely and Stewart’s paper for more on these.
  • Hit p to have SkyLum.exe calculate points after you have positioned the sun and chosen a skymodel
  • Hit o to have SkyLum write a sky model file. It’s conventional to give it a .csv extension.
  • Use a text editor to delete the header lines. All you want left are the comma-separated lines with the azimuth, elevation and weight, as shown above.

Now, how to run RCP?

Let’s assume you have a DEM  called myDEM.tif, and an illumination file you made with SkyLum called illum.csv. You’ve already deleted the header lines from illum.csv. You want to divide your DEM into 1000x1000px chunks, with an overlap of 200 px. We’ll also assume you have 4 processors and you want to use 3 of them for this operation.

Drop both myDEM.tif  and illum.csv in the rcp_master directory. Open a terminal there.

The general form of the RCP command line is:

python3 raster_chunk_processing.py -m [method] {general options} {method-specific options} input_file output_file

Notice that the first element in the command line is python3, not just python.

For my skymodel in this case the command line will be…

 python3 raster_chunk_processing.py -m skymodel -s 1000 -o 200 -p 3 -l illum.csv --verbose myDEM.tif myDEM_skymodel.tif

What you’ll see in the terminal is a bunch of information about the job to be done, and then you’ll see RCP submitting the sub-jobs (the chunks) to be skymodelled. You can walk away, or watch, fascinated, as the chucks get worked on…

Preparing output file myDEM_skymodel.tif...
Output dimensions: 4046 rows by 4515 columns.
Output data type: <class 'numpy.float32'>
Output size: 69.7 MiB
Output NoData Value: -32767.0

Processing chunks...
Tile 0-0: 1 of 25 (4.000%) started at 0:00:00.701178 Indices: [0:1200, 0:1200] PID: 8777
Tile 0-1: 2 of 25 (8.000%) started at 0:00:00.701253 Indices: [0:1200, 800:2200] PID: 8778
Tile 0-2: 3 of 25 (12.000%) started at 0:00:00.701284 Indices: [0:1200, 1800:3200] PID: 8779

All of the switches and parameters are explained quite well in Jake’s post. My additional notes are:

  • the input file is generally a geotiff, but it can also be a TIF/TFW pair. It has to have a NoData value defined or you will get an error No NoData value set in input DEM.
  • RCP will exit with an error if the output file already exists
  • By default RCP applies a vertical exaggeration of 5X to the DEM, because this was what Kennely and Stewart did in their paper. However, you can change this if you would prefer a different vertical exaggeration. Open methods.py in your text editor and go to line 272. This line originally says  in_array *= 5. You can change that 5 to whatever number you would prefer. (No recompile necessary.)
  • the overlap (-o) parameter is based on how far you think shadows may stretch. The shadowing algorithm checks outward for 600 pixels to see if a given point is shadowed by anything. For this reason, there is no point in making overlap larger then 600.
  • the chunk size parameter (-s) is based on how much RAM each process requires while running. You can experiment with this and watch on the system monitor to see how close you are to spilling over into swap.
  • When shadows are being calculated in the skymodel,

The output file can be dragged into QGIS, where I often find I want to increase brightness.

adjusdting_brightness
Left: hillshade image fresh out of RCP. Right: with brightness +100

Another idea that seems to produce some nice results is to combine two or more skymodels with different transparencies and styling. You might also want to check out a gallery of the results of applying all of the different skymodels in SkyLum to the same piece of terrain.

Now that you have RCP running, of course you’ll need to try all of the different skymodels on your favourite DEM to see which you like best. But it’s also worth checking out the other kinds of processing RCP can do on DEMs. It is a very fast conventional hillshader (method hillshade). It runs a nice, quick gaussian blur (method blur_gauss), much faster than the SAGA Gaussian filter module in QGIS. And I haven’t even tried the Clahe and TPI processing yet!