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!
Ok, this doesn’t have a direct value for whatsoever, but it is just an example of how to have fun with the visualization tools in MATLAB. Note that you’ll need to download the function plot3d from github to run this example.
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM,'preprocess','carve'); S = STREAMobj(FD,'minarea',1000); S = klargestconncomps(S); RGB = imageschs(DEM,DEM,'colormap','landcolor'); [x,y] = getcoordinates(DEM); surface('XData',x,'YData',y,... 'ZData',repmat(min(DEM),numel(y),numel(x)),... 'CData',double(RGB)/255,... 'FaceColor','texturemap','EdgeColor','none'); hold on plot3d(S,DEM) axis image exaggerate(gca,2) axis off hold off