### Chimaps in a few lines of code (4)

Posted on Updated on

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.

```flowpathapp(FD,DEM,S)
```

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.

`chiplot(S,DEM,A,'mnplot',true,'plot',false)`

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.

## 9 thoughts on “Chimaps in a few lines of code (4)”

Rishav Mallick said:
September 25, 2017 at 6:58 am

Hi Wolfgang,
While using chiplot to find the optimum m/n ratio for small streams, i found that the resulting structure C (created as C = chiplot(S,DEM,A);) has a predicted value as C.pred. When i plotted this value (C.distance vs C.pred) versus C.distance vs C.elevbl i found that C.pred has a a very angular shape (an obtuse angle) , especially when theta values are > 2, as opposed to the gentle logarithmic shape i expected. Is that expected? Please can you comment on that?
I was also thinking, given how noisy SRTM elevation data is (with its stair-step features), chiplot completely removes the need to compute spatial gradients. And because of that, there is no need to do spatial smoothing or resampling of the data at fixed contour lines. I think this is a very neat method 🙂 But i wonder if results of channel steepness and concavities of individual streams calculated using an integral method will be vastly different when compared to slope-area regression?

cheers,
RIshav

wschwanghart responded:
September 25, 2017 at 7:59 am

Dear Rishav,
I have tried following code to be able to reproduce your findings:

```DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM);
S  = STREAMobj(FD,'minarea',1e6,'unit','map');
S = klargestconncomps(S);
A = flowacc(FD);
C = chiplot(S,DEM,A,'mn',0.45);
plot(C.distance,C.pred)

C = chiplot(S,DEM,A,'mn',2);
figure
plot(C.distance,C.pred)
```

I do not get the angular pattern that you describe, but the predicted elevations seem at odds with true elevations for a mn-ratio of 2. The mn-ratio usually ranges between 0.2 to 0.8 (often it is set to 0.45), so this is something that I’d expect as such a high mn-ratio is unable to linearize the relation between chi and elevation.
You are completely right with the robustness of chi analysis to errors in DEMs. This is something that Royden and Perron describe in their paper as a major motivation for developing this approach. Unless you are dealing with a lot of scatter, smoothing is usually not necessary for chi analysis. Whether results from chi analysis are vastly different depends on the noise in the slope area regression. Given the robustness of the chi approach, I’d always rely on mn ratios derived from chi analysis rather than those obtained from slope area plots.
Cheers, Wolfgang

Simone said:
December 7, 2017 at 8:40 pm

Hi Wolfgang,
Is there a way to have the same plot replacing the chi index in horizontal coordinates with the ksn index, in order to obtain a ksn/z plot for different theta values?

wschwanghart responded:
December 8, 2017 at 9:37 am

Hi Simone, sure, this is possible though not out of the box. But here is the code:

```DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','c');
S = STREAMobj(FD,'minarea',1000);
S = trunk(klargestconncomps(S));
A = flowacc(FD);

theta = [0.3:0.05:0.6];
for r = 1:numel(theta)
k{r} = ksn(S,DEM,A,theta(r));
end
[~,~,z,k{:}] = STREAMobj2XY(S,DEM,k{:});
for r = 1:numel(theta); plot(k{r},z); hold on; end
```

Cheers, Wolfgang

simr87 said:
December 8, 2017 at 3:54 pm

Thanks!
The last question: how can I modify the script in order to obtain a plot with x for horizontal and ksn for vertical coordinates? Sorry for the questions, I need both ksn/z and x/ksn plots for different theta values but my matlab language knowledge is very bad

wschwanghart responded:
December 8, 2017 at 5:49 pm

Hi Simone, just modify the last two lines:

```[~,~,z,d,k{:}] = STREAMobj2XY(S,DEM,S.distance,k{:});
subplot(1,2,1)
for r = 1:numel(theta); plot(k{r},z); hold on; end
subplot(1,2,2)
for r = 1:numel(theta); plot(d,k{r}); hold on; end
```
simr87 said:
January 11, 2018 at 10:02 am

Hi Wolfgang,
I’m trying to write a script to plot in the same page long-river profiles and, for different theta values, chi-profiles, ksn/z and d/ksn profiles.
Here my script:

DEM = GRIDobj(‘srtm_bigtujunga30m_utm11.tif’);
FD = FLOWobj(DEM,’preprocess’,’c’);
S = STREAMobj(FD,’minarea’,1000);
S = trunk(klargestconncomps(S));
A = flowacc(FD);

theta = [0.4:0.1:0.6];
for r = 1:numel(theta)
k{r} = ksn(S,DEM,A,theta(r));
k{r} = smooth(S,k{r},’K’,200);
c{r} = chitransform(S,A,’mn’,theta(r));
end

[~,~,z,d,k{:}] = STREAMobj2XY(S,DEM,S.distance,k{:});
subplot(2,2,1)
plotdz(S,DEM);
subplot(2,2,2);
for r = 1:numel(theta); plot(c{r},z); hold on; end
subplot(2,2,3)
for r = 1:numel(theta); plot(k{r},z); hold on; end
subplot(2,2,4)
for r = 1:numel(theta); plot(d,k{r}); hold on; end

The problem is with chi profile that doesn’t work.What’s wrong?

Add. question: How can I smooth the long – river profile?

wschwanghart responded:
January 11, 2018 at 5:20 pm

Hi,
c{:} is a cell array of node-attribute lists. There are two options:

One:

```[~,~,z,d,k{:},c{:}] = STREAMobj2XY(S,DEM,S.distance,k{:},c{:});
```

Then you can use your code that follows.

Two:
Use the node-attribute list directly in the function plotdz

```plotdz(S,z,'distance',c{1});
```

Setting the ‘distance’ parameter in plotdz allows you to use custom distances, and not flow distances from the outlet.

You can use the function crs to smooth river profiles. See our recent paper here: https://www.earth-surf-dynam.net/5/821/2017/ , and this blog entry: https://topotoolbox.wordpress.com/2017/12/07/paper-out-now-bumps-in-river-profiles/

Hope this helps.
Wolfgang

simr87 said:
January 11, 2018 at 5:37 pm

Hi Wolfgang, here how I have solved my problem:

DEM = GRIDobj(‘srtm_bigtujunga30m_utm11.tif’);
FD = FLOWobj(DEM,’preprocess’,’c’);
S = STREAMobj(FD,’minarea’,1000);
S = trunk(klargestconncomps(S));
A = flowacc(FD);
z90 = crs(S,DEM,’K’,10,’tau’,0.9);
theta = [0.4:0.1:0.6];
theta2 = [0.3:0.1:0.6];
for r = 1:numel(theta)
k{r} = ksn(S,DEM,A,theta(r));
k{r} = smooth(S,k{r},’K’,200);
end
for s = 1:numel(theta2)
c{s} = chitransform(S,A,’mn’,theta2(s));
c{s} = smooth(S,c{s},’K’,200);
end
[~,~,z,d,k{:}] = STREAMobj2XY(S,DEM,S.distance,k{:});
[~,~,z,d,c{:}] = STREAMobj2XY(S,DEM,S.distance,c{:});
subplot(2,2,1)
plotdz(S,z90);
subplot(2,2,2);
for s = 1:numel(theta2); plot(c{s},z); hold on; end
subplot(2,2,3)
for r = 1:numel(theta); plot(k{r},z); hold on; end
subplot(2,2,4)
for r = 1:numel(theta); plot(d,k{r}); hold on; end