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.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.gdaldem hillshade in.tif out.tif
gdaldem hillshade -multidirectional in.tif out.tif
gdaldem hillshade -combined in.tif out.tif
gdaldem slope -p in.tif out.tif
-p means return the value in percent slope instead of slope in degrees.gdaldem TPI in.tif out.tif
gdal_translate -of EHdr -ot Float32 in.tif temp.flt
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.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.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.
Hello Jacob,
Nice explanation and it's working. Thank you a lot, Jacob.
Happy year-end holiday :)
Best Regards,
Yunianto
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!
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!
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.
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.
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
qto view the in-program help, which describes all the different commands. My usual process is this:+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."pto generate all the points.oto 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 skymodelas 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, liketurbid_315_45_250.csvfor 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.