Skip to content

Instantly share code, notes, and snippets.

@jacobdadams
Last active December 1, 2022 16:39
Show Gist options
  • Save jacobdadams/d00215d5424a82efc07bd26c5325560c to your computer and use it in GitHub Desktop.
Save jacobdadams/d00215d5424a82efc07bd26c5325560c to your computer and use it in GitHub Desktop.
Artistic Terrain Commands and Notes

Setting up a GDAL conda environment

Create an environment

conda create -n <environment_name>

Activate the environment

activate <environment_name>

Install gdal from conda-forge

conda install gdal -c conda-forge

Terrain commands

Adding Overviews

gdaladdo -ro in.tif 2 4 8 16

  • -ro means open the tif in read only mode, which results in an external overview file in.tif.ovr, just like ArcGIS. Internal tif overviews are tricky, and I try to avoid them.
  • The 2 4 8 16 are the pyramid levels. I find these four levels to be sufficient for most of my state-wide to local maps, but this depends on your extent and the resolution of your raster.

Basic Hillshade

gdaldem hillshade in.tif out.tif

Multi-directional Hillshade

gdaldem hillshade -multidirectional in.tif out.tif

Combined Hillshade

gdaldem hillshade -combined in.tif out.tif

Slope

gdaldem slope -p in.tif out.tif

  • -p means return the value in percent slope instead of slope in degrees.

Topographic Position Index (TPI)

gdaldem TPI in.tif out.tif

Texture Shading

  1. gdal_translate -of EHdr -ot Float32 in.tif temp.flt
    • Translates our input tif to a .flt
  2. texture.exe <texture_level> temp.flt temp_texture
    • <texture_level> is a number between 0 and 2. Use 2 to get extreme edges like Tom Patterson uses for his canyon walls and alpine hatching.
  3. texture.exe <vertical_enhancment> temp_texture final_texture.tif
    • <vertical_enhancement> is a number, usually between -2 and 5, that controls the contrast. Play around for artistic value.

Skymodel using rcp.py (github.com/jacobdadams/rcp)

  1. Generate a luminance CSV from SkyLum.exe.
  2. Delete the header rows so that the first row of the csv is the first row of values.
  3. python path/to/rcp.py -s 800 -o 600 -p 6 -m skymodel -l path/to/luminance.csv in.tif out.tif
    • -s 800 uses 800 x 800 pixel chunks. Setting this higher results in exponentially slower runtimes.
    • -o 600 spececifies an overlap of 600 pixels on neighboring chunks (the skymodel algorithm shades out to a max of 600 pixels, so setting this lower than 600 will result in visible seams).
    • -p 6 processes six chunks in parallel. Set this based on the number of logical cores in your system; I recommend using n-1 or n-2 cores if you're going to be doing anything else while processing.
    • -m skymodel tells it to create a skymodel. rcp.py has several other tools, including and adjustable TPI algorithm and several smoothing algorithms.
    • -l luminance.csv is the path to the luminance csv generated and cleaned up earlier.
  4. Wait. A state-wide-plus-50-miles, 10-meter-resolution file using 250 light points (specified in SkyLum.exe) took over three days running on 11 or 12 cores of an Intel i7-8700. The roughly 9600 x 6700 pixel tif in the livestream took 50 minutes with the same luminance csv.
@YuniantoWijayanto
Copy link

YuniantoWijayanto commented Dec 25, 2020

Hello Jacob,

Nice tutorial. I tried to do the texture shading process but stopped in the script language with notifications "'texture.exe' is not recognized as an internal or external command". Please more advise about this. Thank you.

Best Regards,
Yunianto
SS

@jacobdadams
Copy link
Author

Good morning!

Texture.exe is the program written by Leland Brown to do texture shading: http://www.textureshading.com/Home.html.

It looks like the command line interpreter can't find the texture.exe program. Make sure you've downloaded it from https://app.box.com/v/textureshading/folder/1327541894 and extracted the zip file. Then, make sure you're either inside the extracted folder or calling the full path to the exe file.

@YuniantoWijayanto
Copy link

Dear Jacob,

Thank you for your response and explanation. Introduce my self, my name is Yunianto. I'm from Indonesia. I have an interest in artistic terrain processing after I saw your presentation on channel YouTube. I usually process terrain data using the processing toolbox in Arc GIS. This is my first time trying to use GDAL tools. Sorry if there are some questions due to not knowing about the GDAL process.

I have successfully carried out the steps above according to your direction. Thanks again Jacob (please find attached file result)
SS

I tried to do the Skymodel process but still confused about doing luminance extraction with SkyLum.exe. Could you tell me the steps in the Skymodel process? I saw the results of this process from your YouTube channel really charming.

Best Regards,
Yunianto

@jacobdadams
Copy link
Author

Hello Yunianto, good to meet you!

You're right, I hadn't realized there's not a lot of documentation on skylum.exe.

When you launch SkyLum.exe, you're presented with an interface for modeling the effect of different lighting conditions.

By default, the window is usually sized taller than your monitor screen and cuts off the important status line at the bottom. Grab the top of the window and resize it so it fits your monitor.

You can type q to view the in-program help, which describes all the different commands. My usual process is this:

  • Right click and select my lighting type.
  • Use the up/down (altitude) and left/right (azimuth) arrow keys to adjust the position of the sun.
  • Use + and - to adjust the number of points (I usually use the default of 250). The more points you use, the longer it takes to process, but the more realistic the lighting output. There's going to be a point of diminishing returns where it looks "good enough."
  • Once your light is set up, hit p to generate all the points.
  • With the points generated, hit o to choose an output name for the csv. Hit enter to generate the csv file.

This is the file you'll use with python rcp.py -m skymodel as described in the main instructions. Don't forget to delete the header lines in the csv so that the first line is the first set of values. Because this removes your reference for how you created the luminance csv, I usually name my csv with my settings, like turbid_315_45_250.csv for a turbid sky, azimuth 315, altitude 45, 250 points.

Because skymodeling takes so long to process, I usually start with a smaller (~2000x2000 pixel) subset of my raster to test different lighting setups before processing the whole raster.

The beauty of skymodelling through either rcp or Kennelly & Stewart's original program is that the final product is a georeferenced raster that can be loaded directly into a GIS program. An alternative method that may be faster/easier is to use the 3d-modelling program Blender to create your scene. Daniel Huffman has a great tutorial on running that at https://somethingaboutmaps.wordpress.com/2017/11/16/creating-shaded-relief-in-blender/. I haven't tried this yet myself, but I suspect that the files produced by Blender aren't georeferenced, so you'd need to do that before using them in a GIS.

@YuniantoWijayanto
Copy link

Hello Jacob,

Nice explanation and it's working. Thank you a lot, Jacob.
Happy year-end holiday :)

Best Regards,
Yunianto

@VicluUiB
Copy link

VicluUiB commented Nov 29, 2022

Hi Jacob,

I know that this is already two years late by now, but I'm hoping that you're still monitoring this, because I too have a question about this workflow.

I'm reading your reply to Yunianto, and I'm starting to think that my suspicion is right as I have not seen mentioning of what I'm about to ask anywhere else. Not even in your reply to Yunianto in December 2020. My question is this:

Do you load your own raster inside Skylum when generating the illumination.csv file, or do you just produce an illumination.csv after having selected your model of choice that is free-standing from the raster that you will later use in RCP? The reason for why I'm asking this is because the readme inside Skylum states:

"1) Use the sample digital elevation model (Mt. Hood data) provided in the geodata folder (“dem_ve5”), or select a DEM to render and multiply it by five to get a DEM with a vertical exaggeration of 5x’s (recommended). Name the output “dem_ve5”."

It was my initial suspicion that one would have to generate an illumination file that took into account the spatial extent and topography of one's own elevation model (I'm still not sure why I came to suspect this though?). I eventually found an ArcGIS toolbox for using Skymodels, albeit with a limited amount of models for illumination points:

https://www.esri.com/en-us/maps-we-love/gallery/gorgeous-gorge

In this toolbox, however, one would simply select a raster as input, and then also a pre-made file with the illumination values. Something with the toolbox went wrong, however, so I wasn't able to inspect the result. Now, of course that file with already pre-made illumination values cannot have been made in reference to one's own raster data, which suggests to me that I'm on the right track with my suspicion. But, you never know.

Perhaps you can tell that I've been overthinking this somewhat. But hopefully you can bring me some clarity!

Best!
/Victor

Oh and by the way, I think the work that you do (both in terms of your skill for sharing your knowledge freely, as well as the stunning visuals
that you produce) is absolutely phenomenal!

@jacobdadams
Copy link
Author

Hello Victor! Thanks for the kind words—I'm a big believer that we all benefit by sharing knowledge. Most of what I've done builds heavily off of others' work.

In regards to the dem used for generating the luminance csv, you're correct in assuming its not using a custom dem. SkyLum.exe uses a small terrain model that lets you visualize different light conditions on example terrain and then generate a set of hill/shadowshade parameters that approximate that look.

The SkyLum folder I got from Kennelly and Stewart and bundled in the rcp repo contains two different example files: MtStHelens.hrz and sample.hrz. Mt St Helens is a good example of mountains terrain, while sample.hrz gives a good test bed for canyons and other physical depressions.

I suppose you could generate your own terrain files and load them in, but I've not tried that. The key is just to see what the light does on the example terrains and then get the csv from that for calculating your entire skymodel using either my rcp.py or the Sky Luminance tool in Esri's Terrain Tools toolbox that you linked above (rcp.py is my own implementation of Kennelly and Stewart's example code, using open-source tools and parallel processing).

I hope that answers your questions!

@VicluUiB
Copy link

VicluUiB commented Dec 1, 2022

Hi Jacob,

It sure did, thank you so much for the assistance!

However, there is now another issue that I hope you can help with...
Through the aid of a colleague of mine, who knows how to use python, the script is now working to the extent that it cuts up the raster in chunks, processes them, and then stitches them back together. So far so good! However, one recurrent issue that we have come across is highly visible seams in the final raster. We have drawn the conclusion that these seams are artefacts from the chunk processing and that their visibility relates to how one specifies the values for both "-o" and "-s". We unfortunately haven't documented our experiments that thoroughly, but the results vary from either having the entire output raster divided with highly visible seams all across, to something that almost looks perfect but with just one or two seams. How should one relate to the values of these parameters? Will there be an issue of
scaling once, or if, one would decide to either increase or decrease the spatial extent of once raster? Or should one simply go about it in a "trial-and-error" kind of way?

Our current test subject is a raster of the coast of western Norway with the following dimensions:
X: 6859 Y: 6105
Pixel size: 25x25m

The current settings for the script:
-p 8
-s 1500
-o 25

I've highlighted a few places in yellow just to illustrate what it looks like.

image

@jacobdadams
Copy link
Author

Yeah, that is a result of the -o value not being 600. The overlap determines how far outside of the chunk it reads data from for algorithms that need to know what's happening beyond the chunk (like skymodelling).

I've not documented this very well on this page, but the shading algorithm for skymodelling looks 600 pixels away from each pixel to determine whether it would be shaded for a given sun altitude and azimuth. Therefore, the overlap value needs to be 600 so that the pixels near the edge get the full shading effect.

I've always used the -s 800 -o 600 settings for 10m x 10m DEMs. You can always try going into the code and changing the 600 pixel value (the max_steps variable at https://github.com/jacobdadams/rcp/blob/master/methods.py#L324) and see what it looks like with a smaller value. Essentially, this will decrease the maximum length shadows can be, but with a 25x25m DEM a smaller value could give the same effective shading as my 600 pixels at 10m.

As you've probably seen, an overlap this increases the runtime significantly. Someday I'll try to optimize this, but that's where we're at right now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment