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.
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 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.hrzandsample.hrz. Mt St Helens is a good example of mountains terrain, whilesample.hrzgives 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!