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.
@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