Today, I read about the Tanaka method, an approach to depicting relief. The method relies on filled contours and coloring contour lines according to their direction and light origin. I have seen solutions in QGIS (here and here) which made me thinking whether the Tanaka method can also be implemented in MATLAB.
One difficulty with Tanaka’s method is that coloring of contour lines requires knowledge about whether higher values are found to the left or the right of a contour line. GDAL’s contouring algorithm has consistent behavior in this respect. However, MATLAB’s algorithm doesn’t. My trick is to calculate surface normals and their angle to the light source, and then interpolate the gridded directions to the vertex locations of the contour lines.
Now here is my approach which uses some undocumented features of the HG2 graphics. An excellent resource for undocumented features in MATLAB is Yair Altman’s blog.
function tanakacontour(DEM,nlevels) %TANAKACONTOUR Relief depiction using Tanaka contours % % Syntax % % tanakacontour(DEM,nlevels) % % See also: contourf, hillshade % % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de) % Date: 8. December, 2017 if nargin == 1; nlevels = 10; end % calculate solar vector azimuth = 315; azid = azimuth-90; azisource = azid/180*pi; [sx,sy,~] = sph2cart(azisource,0.1,1); % get coordinate matrices and z values [Z,X,Y] = GRIDobj2mat(DEM); % calculate surface normals [Nx,Ny] = surfnorm(Z); N = hypot(Nx,Ny); Nx = Nx./N; Ny = Ny./N; H = GRIDobj(DEM); H.Z = Nx*sx + Ny*sy; % Get coordinate matrices and Z [~,h] = contourf(X,Y,Z,nlevels); axis image drawnow % colormap cmap = gray(255)*255; cval = linspace(0,1,size(cmap,1))'; % the undocumented part for r = 1:numel(h.EdgePrims) % the EdgePrims property contains vertex data x = h.EdgePrims(r).VertexData(1,:); y = h.EdgePrims(r).VertexData(2,:); % interpolate hillshade grayscale values to contour vertices c = interp(H,double(x(:)),double(y(:))); % interpolate to colormap clr = interp1(cval,cmap,c); % adjust color data to conform with EdgePrims ColorData property clr = clr'; clr = uint8([clr;repmat(255,1,size(clr,2))]); % set all properties at once set(h.EdgePrims(r),'ColorBinding','interpolated','ColorData',clr,'LineWidth',1.2); end
Ok, now let’s try this on some data. In this code, I have used GRIDobj, the class for gridded, georeferenced data.
[X,Y,Z] = peaks(500); tanakacontour(GRIDobj(X,Y,Z),10)
And now for the Big Tujunga DEM:
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); tanakacontour(DEM,5)
Now that looks ok, but I am not totally convinced by these unrealistic terraces, to be honest. Moreover, my approach so far doesn’t include changing line widths. Contour lines should actually be thicker for slopes that face both towards and away from the light source. However, I didn’t figure out how to do this so far. Any ideas?
Some changes and additions to TopoToolbox are subtle and may not be visible at first glance. Hence, I report them here. One such addition is part of version 2.2 (which Dirk and I are about to finalize shortly) and is found in the folder IOtools (input/output tools): readopentopo. The function uses opentopography’s RESTful Web service for global raster data access for seamless download of global digital elevation models. readopentopo offers access to SRTM-3, SRTM-1, and ALOS 3D World data at your finger tips. By default, the downloaded *.tif-file is stored in your system’s temporary folder and is deleted as soon as it is read into MATLAB. If you want to permanently save the file on your hard drive make sure to define a file path and to set the ‘deletefile’-option to false. The tif-file comes in geographic coordinates. Hence, don’t forget to reproject the DEM to UTM (see reproject2utm) before you start working with it.
Many thanks to members of the opentopography-team who make the data and service available. Please use these services reasonably and don’t excessively download data!
Today’s guest post comes from Philippe Steer, assistant professor at the University of Rennes. Philippe contributed two new functions and I asked him to describe them here to show potential applications. Thanks, Philippe, for this contribution!
I recently needed to compute river network orientation. Indeed river network orientation can reveal distributed tectonics deformation at the vicinity of strike-slip faults (e.g. Castelltort et al., 2012) or can enable to identify developing vs paleo river networks. One way to obtain orientation is first to extract the river segments that are linking confluences to outlets or confluences to channel heads and then to compute the orientation of the resulting segments. The streampoi function allows to extract these feature points of the network such as channel heads, confluences, b-confluences (the points just upstream any confluence) and outlets. To connect these feature points into segment, one can use the propriety that each segment belongs to a unique watershed that is given by the function drainagebasins using as outlets the true river outlets and the b-confluences. By definition, each segment also corresponds to a unique value of strahler order.
Once the river segments are identified, segment orientation is easily obtained. What is interesting is that defining these segments also paves the way to compute several other geometric metrics averaged over the river segment such as segment direct or flowpath length, sinuosity, direct or flowpath slope, mean drainage area and even steepness index (if assuming a value of concavity).
To cleanly do the job, I have written two MATLAB functions:
Here is a piece of code showing you the use of these functions (assuming you have already loaded the DEM of interest):
%% Flow direction and drainage area DEMf = fillsinks(DEM); FD= FLOWobj(DEMf); A = flowacc(FD); %% Extract rivers % Set the critical drainage area A_river=1e7/(res^2); % m^2 W = A>A_river; % Extract the wanted network S = STREAMobj(FD,W); %% River segments % Extract segment and compute their orientations mnratio=0.5; [segment]=networksegment(S,FD,DEM,A,mnratio); % Plot results plotsegmentgeometry(S,segment)
To illustrate the results, I thought the South Island of New-Zealand was a good example. Indeed, as first suggested by Castelltort et al. (2012) the river catchments at the south of the Alpine fault are tilted, possibly suggesting that they have recorded a distributed shear deformation. This tilting is also quite obvious when looking at the river segments in the following figure. Another interesting feature is that the river segments north of the Alpine fault, facing the tilted rivers, display a higher steepness than their south neighbours. This could be due to a long list of forcings, including active tectonics, orographic precipitation gradient or the transition from glacial to fluvial erosion.
Castelltort, S et al. (2012). River drainage patterns in the New Zealand Alps primarily controlled by plate tectonic strain. Nature Geoscience, 5(10), 744-748. [DOI: 10.1038/ngeo1582]
The Eiger north face (46.34°N, 8.00°E) is one of the most spectacular climbing routes in the Alps. Loose rocks, avalanches and severe weather conditions make it particularly challenging. And, as pertinent to north faces, it is shadowed during most of the daytime.
When and where you can catch a glimpse of sunlight, you can calculate using a sun position calculator, TopoToolbox, and a digital elevation model (DEM). A sun position calculator has been made available by Vincent Roy on the file exchange. The function calculates for a given latitude, longitude, date and time the azimuth and zenith angle of the sun. You can then use this input into the function castshadow which calculates the shadowed areas of a DEM for any azimuth and altitude (90°- zenith angle).
The animation below shows the Eiger north face on September 1, 2015. I looped through the hours from 5:00 to 19:00 and saved the figures as gif movie. Of course, there are edge effects that cause that parts of the DEM close to the DEM borders lack the cast shadow of mountains nearby. But the Eiger north face should be captured quite right.
Do you have any applications for the function castshadow? Let me know and comment!
In a just published article in Geology by Jan Blöthe, Oliver Korup and me, we show that the Himalaya-Karakoram Range exhibits an intriguing vertical stacking of erosional domains driven by mass-wasting and (peri-)glacial activity. We show that large bedrock failures preferably detach below the limit of sporadic alpine permafrost adjusting hillslope inclinations to rapid fluvial incision and postglacial stress fields in oversteepened valley flanks. Intense frost cracking above the permafrost limit in turn may result in high-frequency/low-magnitude rockfalls and avalanches, thus limiting the occurrence of low-frequency/high-magnitude rock-slope failures.
As part of the analysis, we introduced ‘excess topography’ as a new measure to characterize mountain topography. Here is an excerpt taken from the supplementary material that describes excess topography and how it is calculated from a DEM:
Excess topography measures for each grid cell in a digital elevation model (DEM) the rock-column height above an arbitrarily defined threshold slope surface. The threshold slope surface is an idealized surface that we derive from a given DEM. The threshold surface only contains inclinations that are lower or equal than the arbitrarily set threshold inclination. The elevations of the threshold surface are constrained by the DEM elevations and the threshold inclination such that (1) the elevation of the threshold and DEM input surface coincide for a given grid cell if the slopes between this grid cell and all others are less or equal the threshold slope; and (2) elevations of the threshold surface fall below the DEM input surface in grid cells rising above the surrounding topography at angles steeper than the threshold inclination. In case (2), we calculate the threshold surface elevation as the minimum elevation of all surrounding locations plus the maximum allowed topographic rise (threshold inclination) along the distance separating the locations.
A picture is worth a thousand words. Here you go:
Excess topography can be calculated using the new function excesstopography in TopoToolbox. The function is quite simple. The most important algorithm that it requires – image morphology with a grayscale structuring element – is made available by the function ordfilt2 in the Image Processing Toolbox. Did you know this function? This function is extremely powerful and I can think of quite a lot of nice applications of DEM analysis. Since I said that graydist is my favorite function, ordfilt2 is my second favorite function.
Blöthe, J.H., Korup, O., Schwanghart, W. (2015): Large landslides lie low: Excess topography in the Himalaya-Karakorum ranges. Geology, 43, 523-526. [DOI:
** edited on May 23, 2015 to update reference **
** edited on June 23, 2015 to include free pdf download **
In my last posts I wrote about carving and filling. As noted, both algorithms are implemented in the constructor function FLOWobj and thus are primarily used to derive flow directions through topographic depressions. These options, thus, don’t change the DEM itself. In this post, I will show how you will carve and fill the DEM.
Filling is easy and implemented in most terrain analysis software. A number of different algorithms exist that perform flood-fill operations; in our paper, Dirk and I compared processing speeds to fill an image. We found that the speed at which one and the same DEM is filled varies widely between different software. The Matlab’s image processing toolbox’s function imfill (the function that is called by TopoToolbox function fillsinks) outperforms other software by 8-10 times faster processing speed. However, I like to stress that speed is not the only indicator of performance and that other software may have better ways to manage memory. However, I don’t want to go into detail here…
Here is how to fill a DEM (the Big Tujunga DEM is included in the distribution of TopoToolbox 2)
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); DEMf = fillsinks(DEM);
One way to judge whether filling is a good option is to calculate the “energy” required to modify the DEM. When we fill the DEM, the energy is the total volume added to the DEM.
sum(DEMf.Z(:)-DEM.Z(:))*(DEM.cellsize)^2 ans = 18801000
A volume of ~19×106 m3 is required to force all flow towards the DEM edges. So how about carving? This requires us to calculate a FLOWobj and use it subsequently in a function called imposemin. Imposemin takes a FLOWobj and a GRIDobj and calculates a downstream minimum operator, i.e., it lets minima flow downstream.
FD = FLOWobj(DEM,'preprocess','carve'); DEMc = imposemin(FD,DEM); sum(DEM.Z(:)-DEMc.Z(:))*(DEM.cellsize)^2 ans = 7848000
Carving the DEM requires removing a volume of ~8×106 m3 only. This is less than 1/2 of the energy that filling needs.
Based on these results, I would better carve the DEM to make it hydrologically correct. And, as John Lindsay pointed out in his comment on a previous post, carving is probably the better method in most situations. There may be, however, situations where filling might prove better, e.g., when the DEM has erroneously low valued, single scattered pixels. I had this problem once when working on a DEM in the Himalayas. This is why you find another function in TopoToolbox called elevateminima. This one allows you to raise single pixels to the height of their surrounding pixels.
Another note: Hybrid methods exist that partially fill and breach. I have not implemented these algorithms until now but perhaps I’ll do so in the future.