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!
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.
A couple of weeks ago, I covered a series of blog posts about chi maps. This post links to these posts for future reference.
- Part 1: Overview and approach to creating chimaps
- Part 2: Preparing the stream network
- Part 3: The mathematical basis
- Part 4: Finding the right m/n ratio
- Part 5: Plotting the chi-map
In addition, I like to draw your attention to the tool ChiProfiler that has been developed by Sean Gallen from the ETH Zürich. ChiProfiler is based on the TopoToolbox functions and provides a visual interface to stream profile analysis and includes ksn and chi analysis, amongst other. Give it a try!
Having prepared a stream network and equipped with a reasonable value of the m/n ratio, we are now ready to plot a chimap that visualizes the planform patterns of chi. The main interest in these maps lies in chi values near catchment divides as large differences between adjacent catchment would indicate a transient behavior of drainage basin reorganization.
Some of you might have already experimented with TopoToolbox to create chimaps. Perhaps you became exasperated with the function chiplot that is restricted to calculations with only one drainage basin and has a bewildering structure array as output. The reason for the confusing output of chiplot is that it is fairly old. At this time, I hadn’t implemented node-attribute lists that are now more common with STREAMobj methods.
Realizing this shortcoming of chiplot, I wrote the function chitransform. chitransform is what I’d refer to as a low-level function that solves the chi-equation using upstream cumulative trapezoidal integration (see the function cumtrapz). chitransform requires a STREAMobj and a flow accumulation grid as input and optionally a mn-ratio (default is 0.45) and a reference area (default is 1 sqkm). It returns a node-attribute list, i.e., a vector with chi values for each node in the STREAMobj. Node-attribute lists are intrinsically tied to the STREAMobj from which they were derived. Yet, they can be used together with several other TopoToolbox functions to produce output.
Ok, let’s derive chi values for our stream network:
A = flowacc(FD); % calculate flow accumulation c = chitransform(S,A,'mn',0.4776);
In the next step, we will plot a color stream network on a grayscale hillshade:
imageschs(DEM,,'colormap',[1 1 1],'colorbar',false,'ticklabel','nice'); hold on plotc(S,c) colormap(jet) colorbar hold off
Interestingly, there seem to be some locations with quite some differences in chi values on either side of the divide. “Victims” seem to be rather elongated catchments draining northwest. Let’s zoom into one of these locations.
Are these significant differences? Well, it seems by just looking at the range of values. However, to my knowledge no approach exists that provides a more objective way of evaluating the significance of contrasting chi values and their implications about rates of divide migration. Still, we now have a nice map that can aid our geomorphic assessment together with the tectonic and geodynamic interpretation of the Mendocino Triple Junction.
Unfortunately, I must leave the discussion to you since I am quite unfamiliar with the region. If anyone wants to share his or her interpretation, I’d be more than happy to provide space here. So far, I hope that these few posts on chimaps were useful to you and sufficiently informative to enable you to compute chimaps by yourself. In my next post, I will give a short summary and show with another example that eventually chimaps can be derived really in a few lines of code.
My last post described the math behind chi analysis. By using the integral approach to the stream power model of incision, we derived an equation that allows us to model the longitudinal river profile as a function of chi. At the end of the last post I stated that this model is a straight line if we choose the right m/n ratio, that is the concavity index. Today, I’ll show how we can obtain this m/n ratio by means of nonlinear optimization using the function chiplot.
Ideally, we find the m/n ratio using a perfectly graded stream profile that is in steady state. I scanned through a number of streams using the GUI flowpathapp and found a nice one that might correspond to Cooskie Creek that Perron and Royden (2013) used in their analysis.
The menu Export > Export to workspace allows saving the extracted STREAMobj St to the workspace. Then, I use the function chiplot to see how different values of the m/n ratio affect the transformation of the river profile. Setting the ‘mnplot’ option makes chiplot return a figure with chi-transformed profiles for m/n ratios ranging from 0.1 to 0.9 in steps of 0.1 width.
Clearly, the chi-transformed profile varies from concave upward for low m/n values to convex for high m/n values. A value of 0.4-0.5 seems to be most appropriate but choosing a value would be somewhat hand-wavy. Thus, let’s call the function again, but this time without any additional arguments.
C = chiplot(S,DEM,A);
Calling the function this way will prompt it to find an optimal value of m/n. In addition, it will return a figure with the linearized profile. By default, the function uses the optimization function fminsearch that will vary m/n until it finds a value that maximizes the linear correlation coefficient between the chi-transformed profile and a straight line. So what is that value? Let’s look at the output variable C:
C = mn: 0.4776 mnci:  beta: 0.1565 betase: 2.2813e-04 a0: 1000000 ks: 114.8828 R2: 0.9988 x_nal: [581x1 double] y_nal: [581x1 double] chi_nal: [581x1 double] d_nal: [581x1 double] x: [582x1 double] y: [582x1 double] chi: [582x1 double] elev: [582x1 double] elevbl: [582x1 double] distance: [582x1 double] pred: [582x1 double] area: [582x1 double]
C is a structure array with a large number of fields. C.mn is the value we are interested in. It is 0.4776 and corresponds nicely to previously reported m/n ratios although it differs from the value 0.36 found by Perron and Royden (2013). Of course, this value might differ from river to river and you should repeat the analysis for other reaches to obtain an idea about the variability of m/n.
Ok, we have the m/n ratio now. In my next blog we will eventually produce the chimap that we initially set out for.
P.S.: We could have also used slope-area plots to derive the m/n ratio (see the function slopearea). However, along-river gradients are usually much more noisy than the profiles and power-law fitting a delicate matter. I’d definitely prefer chi-analysis over slope-area analysis.
In my previous posts (here and here), I introduced chimaps and how they are calculated in TopoToolbox. I showed that deriving a suitable stream network requires some work. Today, I will cover the dry part: the math behind.
Chi-analysis is predicated on the well-known stream power law (SPL). I will not detail the SPL but refer you to the paper of Dimitri Lague (Lague 2014). In its simplest form, the SPL shows that the energy expenditure required to incise into the stream bed is a function of stream gradient (S) and upslope area (A) and some parameters k, m and n. Combining the SPL with a simple vertical uplift rate U allows us to derive a very simple landscape evolution model for the fluvial domain:
Any change in elevation z with time t is caused by an imbalance of uplift and incision. If both processes are perfectly balanced, there is no change in elevation. We are in a state of a dynamic equilibrium or steady state:
We can rearrange this equation and solve for S=dz/dx where x is the stream distance measured from the outlet. At this time I denote that upslope area is a function of x. U and k are assumed to be spatially uniform.
The trick is then to take the integral of (dz/dx) with respect to x from a base level z(x0) (the elevation at the outlet). For convenience of units, we also introduce a reference area A0.
Now that looks complicated. Yet, by replacing the left-hand integral with z and the right-hand integral with
we can simply rewrite this equation as a function of a straight line:
A straight line? Yes… but only if we choose the right m/n ratio. I will demonstrate in the next post how to do this.
P.S.: This was really a short derivation of chi. I can recommend the paper of Perron and Royden (2013) for a more comprehensive account.
Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Process. Landforms, 39(1), 38–61, doi:10.1002/esp.3462, 2014.
Perron, J. T. and Royden, L.: An integral approach to bedrock river profile analysis, Earth Surf. Process. Landforms, 38(6), 570–576, doi:10.1002/esp.3302, 2013.