I just need to finish the year with a quick blog entry… And this one is – once again – about chi maps.
I thought about different ways to display chi maps. chi maps are usually plotted using stream networks. However, to this end, chi maps are not so much about stream networks but rather about drainage divides and how they potentially migrate in the future. In addition, the distance of potential divide migration depends on the relative difference of chi values on either side of the divide, not so much on the absolute values. How could we create a map that reflects this information
Here is an idea: First, we create a chi transform of the river network (see function chitransform). Then we map these chi-values to the entire grid. This can be accomplished with the new function mapfromnal which lets you map values from a stream network to its nearest upstream contributing pixels. Then, we will derive some-pixel wide fringes of chi values around the divides and calculate their spatial average. Dividing the chi values by these spatial averages then quantifies the degree to which some part of the divide is either a victim or an aggressor. Since victims have higher values, I call this metric the chi-inferiority.
Now here is the code:
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM); S = STREAMobj(FD,'minarea',500); S = removeedgeeffects(S,FD); [~,zb] = zerobaselevel(S,DEM); A = flowacc(FD); c = zb + chitransform(S,A); C = mapfromnal(FD,S,c); D = drainagebasins(FD,S); for r = 1:max(D) I = D==r; I = erode(I,ones(21)) & I; C.Z(I.Z) = nan; end C1 = dilate(C,ones(21)); C2 = erode(C,ones(21)); C = C/((C1 + C2)/2); imageschs(DEM,C,'colormap',ttscm('vik'),... 'colorbarylabel','\chi-inferiority',... 'ticklabels','nice'); hold on plot(S,'k') hold off
Happy New Year everone!
In previous posts (here, here, here), I have covered some of the tools and applications of smoothing in the analysis of river networks. Thereby, I have focused on the analysis of river profiles (see also our paper here). However, these tools can also be used to smooth and analyse the planform patterns of river networks.
Here is how:
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM,'preprocess','carve'); S = STREAMobj(FD,'minarea',1000); S = klargestconncomps(S); y = smooth(S,S.y,'k',500,'nstribs',true); x = smooth(S,S.x,'k',500,'nstribs',true); [~,~,xs,ys] = STREAMobj2XY(S,x,y); subplot(2,1,1) plot(S) subplot(2,1,2) plot(xs,ys)
Now what is this good for? First, simplification by smoothing might help in visualization. Second, the smoothed river coordinates can be used to derive metrics such as planform curvature of river networks. These metrics are strongly affected by the circumstance that river networks derived from digital elevation models are forced to follow the diagonal and cardinal directions of the grid. This results in local flow directions in steps of multiples of 45°. Changes in flow direction thus occur in jumps, and accordingly planform curvature of river networks are highly erratic. These erratic values may hide some of the larger-scale curvature patterns that we might be interested in.
The STREAMobj/curvature function features an example that shows how to accentuate larger-scale curvature patterns using smoothing. Following up on the code above, deriving curvature from smoothed planform coordinates can be done as follows:
y = smooth(S,S.y,'k',100,'nstribs',true); x = smooth(S,S.x,'k',100,'nstribs',true); c = curvature(S,x,y); plotc(S,c) box on % center 0-value in colormap caxis([-1 1]*max(abs(caxis))) colormap(ttscm('vik')) h = colorbar; h.Label.String = 'Curvature';
Clearly, the amount of smoothing determines the spatial scale of curvature. Here, the chosen smoothing parameter (K=100) accentuates the bends at a scale of a few hundred meters.
There are a number of applications for this kind of analysis. Lateral erosion of rivers may undercut slopes and it may be interesting to investigate if planform curvature metrics help to detect and predict locations of landsliding. Just an idea…
TopoToolbox has a new tool. Ok, it’s a small tool with limited functionality, but it might be helpful if your aim is to map points along river networks. “Mapping points can be easily done in any GIS!”, you might think. True, but as I am currently working on knickpoints in river profiles (see also the automated knickpointfinder), I wrote this tool to quickly map knickpoints both in planform and profile view of the river network.
So, here is mappingapp. Give it a try.
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM); S = STREAMobj(FD,'minarea',1000); mappingapp(DEM,S)
Currently, the GUI is limited to a few basic tools. You can map a single point which automatically snaps to the river network S. The zoom tools allow you to navigate the DEM. You can add a new point using the + button. This will make the previous point permanent and add a new row to the table. Finally, the table can be exported to the workspace. The table contains the coordinates and z-values as well as a column IXgrid. IXgrid contains the linear indices into the DEM.
I recently mentioned that MATLAB now lets you easily use different colormaps in one figure. The trick is to provide the axis handle as first input argument to the colormap function call.
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); ax(1) = subplot(2,1,1); imagesc(DEM); colorbar; ax(2) = subplot(2,1,2); imagesc(DEM); % Use the axis handle as first input argument % to colormap colormap(ax(2),landcolor); colorbar
However, how do you use multiple colormaps in one axes? That’s a little more tricky. Below, you find commented code that shows how to plot a colored and hillshaded DEM together with a stream network colored with ksn-values.
% Some data (DEM and ksn of river network) DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM); A = flowacc(FD); S = STREAMobj(FD,'minarea',1000); k = smooth(S,ksn(S,DEM,A),'K',1000); % Create a new figure hFig = figure; % Call imageschs which has some parameter/value pairs % to control colormap and colorbar appearance imageschs(DEM,,'colormap','landcolor',... 'colorbar',true,... 'colorbarylabel','Elevation'); % Get the current axes h_ax = gca; % Create a new axes in the same position as the % first one, overlaid on top. h_ax_stream = axes('position', get(h_ax, 'position'),... 'Color','none'); % Plot the stream network plotc(S,k,'LineWidth',2); % Make axis invisible set(h_ax_stream,'Color','none'); % Create a colorbar hc = colorbar(h_ax_stream,'Location','southoutside'); % ... and label it hc.Label.String = 'ksn'; % Adjust color range caxis(h_ax_stream,[0 200]) % Perfectly align both axes set(h_ax_stream,'position', get(h_ax, 'position'),... 'XLim',xlim(h_ax),... 'YLim',ylim(h_ax),... 'DataAspectRatio',[1 1 1]... ); % Link both axes so that you can zoom in and out linkaxes([h_ax,h_ax_stream],'xy'); % Resizing the figure may screw it all. Using the figure % property ResizeFcn makes sure both axes remain perfectly % aligned hFig.ResizeFcn = @(varargin) ... set(h_ax_stream,'Position',get(h_ax, 'position'));
The trick is to use additional axes. The difficulty lies in perfectly aligning them and to make alignment robust against resizing the figure.
Did you know that TopoToolbox has a function to identify knickpoints in river profiles? Well, if you don’t know, now you know.
The function is called knickpointfinder. It uses an algorithm that adjusts a strictly concave upward profile to the actual profile. Offsets between the actual and the concave upward profile occur where the actual profile has convexities. Relaxing the concavity constraint where offsets attain a maximum will adjust the concave profile to the actual profile. knickpointfinder adjusts the profile iteratively until offsets fall below a specified tolerance value. Look at the animation below which probably explains more than a thousand words.
Using knickpointfinder is easy. Just see the function help to run the example whose results are shown above.
Let us know how well knickpointfinder suits your needs. Note that this algorithm is not yet published, so please give credits to our TopoToolbox paper if you are using this algorithm in your work.
River profiles are concave upward if they are in a dynamic equilibrium between uplift and incision, and if our simplified assumptions of steady uplift and the stream power incision law (SPL) hold. The concavity derives from the SPL which states that along-river gradients S are proportional to upslope area A exponentiated by the negative mn-ratio.
I have mentioned the mn-ratio several times in this blog. Usually, we calculate it using slope-area plots or chi analysis both of which are included in TopoToolbox. However, these methods usually lack consistent ways to express the uncertainties of the mn-ratio. The lack of consistency is due to fitting autocorrelated data which elude a straightforward statistical analysis.
Today, I want to present a new function that uses Bayesian Optimization with cross-validation to find a suitable mn-ratio. While Bayesian Optimization is designed to find optimal values of objective functions involving numerous variables, solving an optimization problem with mn as the only variable nicely illustrates the approach.
Bayesian Optimization finds a minimum of a scalar-valued function in a bounded domain. In a classification problem, this value could be the classification loss, i.e., the price paid for misclassifications. In a regression problem, this value might refer to the sum of squared residuals. The value might also be derived using cross-validation, a common approach to assess the predictive performance of a model. Such cross-validation approaches might take into account only random subsets of data, which entails that the value to be optimized might not be the same for the same set of input parameters. Bayesian Optimization can handle stochastic functions.
Now how can we apply Bayesian Optimization for finding the right mn-ratio? The new function mnoptim uses chi analysis to linearize long-river profiles. If there are several river catchments (or drainage network trees), the function will pick a random subset of these trees to fit a mn-ratio and then tests it with another set of drainage basins. This allows us to assess how well an mn-ratio derived in one catchment can actually be applied to another catchment. The goal is to derive a mn-ratio that applies best to other catchments.
Now let’s try this using the new function mnoptim. Here is the code that I’ll use for entire SRTM-3 DEM of Taiwan. I’ll clip the stream network to the area upstream of the 300 m contour to avoid an influence of the alluvial low-lying reaches.
DEM = GRIDobj('taiwan.tif'); FD = FLOWobj(DEM); S = STREAMobj(FD,'minarea',1e6,'unit','map'); C = griddedcontour(DEM,[300 300],true); S = modify(S,'upstreamto',C); A = flowacc(FD); [mn,results] = mnoptim(S,DEM,A,'optvar','mn','crossval',true); % we can refine the results if we need results = resume(mn); % and get an optimal value of mn: bestPoint(results) ans = table mn _______ 0.41482
Now this nicely derives the optimal mn-value of 0.415 which is close to the often reported value of 0.45. Moreover, based on the plot, we gain an impression of the uncertainty of this value. In a transient landscape with frequent knickpoints, the uncertainty about the mn-ratio will be probably larger.
Note that mnoptim requires the latest version of MATLAB: 2017b, as well as the Statistics and Machine Learning Toolbox. It also runs with 2017a, but won’t be able to use parallel computing then.
Analyses that use upslope area usually demand that catchments are completely covered by the DEM. Values of upslope areas may be too low if catchments are cut along DEM edges, and so are estimates of discharge. How can we avoid including these catchments into our analyses?
In one of my previous posts on chi analysis, I showed a rather long code to include those catchments that have 20% or less of their outlines on the DEM edges. Here, I’ll be more strict. I’ll remove those parts of the stream network that have pixels on DEM edges in their upslope area. By doing this, we make sure that drainage basins are complete which is vital for chi analysis or estimating discharge.
Here is the code:
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM); S = STREAMobj(FD,'minarea',1e6,'unit','map'); I = GRIDobj(DEM,'logical'); I.Z(:,:) = true; I.Z(2:end-1,2:end-1) = false; I = influencemap(FD,I);
The influencemap function takes all edges pixels and derives those pixels that they would influence downstream. Again, this is going to be a logical GRIDobj where true values refer to pixels that might potentially be affected by edge effects. Let’s remove those pixels from the stream network S using the modify function.
S2 = modify(S,'upstreamto',~I); D = drainagebasins(FD,S2); imageschs(DEM,,'colormap',[1 1 1],'colorbar',false) [~,x,y] = GRIDobj2polygon(D); hold on plot(x,y); plot(S,'k') plot(S2,'b','LineWidth',1.5)
Ok, now you have a drainage network that won’t be affected by edge effects for further analysis.