Chimaps in a few lines of code (2)

Posted on Updated on

In my last post, I started my series on chimaps. In this and the following posts I will go through the different steps of producing chimaps using TopoToolbox. These steps are:

  1. Load DEM and derive flow directions (FLOWobj).
  2. Choose a suitable threshold for stream initiation and derive a stream network (STREAMobj).
  3. Make sure that all outlets have the same elevation.
  4. Make sure that all drainagebasins are complete, i.e. that the DEM covers their entire extent. Remove incomplete drainage basins.
  5. Choose a suitable value for the mn-ratio.
  6. Calculate the transformed variable chi (Hint: Don’t use chiplot, but check the low-level function STREAMobj/chitransform!).
  7. Plot the results.

In this post I will cover step 1-4. I will use a DEM of the Mendocino region (northwestern California) provided by Amani Al Abri from the University of Boston. Her inquiry has actually prompted me to write about chimaps. The Mendocino region hosts the Mendocino Triple Junction where the Gorda plate, the North American plate, and the Pacific plate meet. The geometry of the plates and their movements have likely favored the development of a slab window and upwelling of astenospheric material, leading to differential patterns of uplift, seismic activity and extrusion of volcanic rock.

1. Load DEM and derive flow directions

DEM = GRIDobj('mendocino.tif')

The first step is to load the DEM. While this is a trivial task, let me briefly mention a few pitfalls one may encounter. While the constructor function GRIDobj tries to identify missing values (actually, it does so quite reliably with the hack of Simon (thanks!)), I usually test this using the function


which provides several information including the minimum and maximum value. If these values are far beyond expected values, they should be replaced by NaNs. If this value was -9999, following code snippet will set these pixels to NaN:


Moreover, the info function shows the coordinate system of the DEM. TopoToolbox requires projected coordinate systems (preferably WGS84 UTM) and equal horizontal and vertical units (e.g. meters). If your DEM comes in geographic coordinates, reproject the DEM with the function reproject2utm. Note also, that you should set missing values before reprojecting the DEM since the bilinear interpolation would generate some blurred regions between non-corrected missing-value regions and regions with valid values.

Missing values usually occur at the boundaries of the DEM or offshore. These are truly missing values that you may want discard. However, data gaps may also exist as individual blobs within the valid DEM domain either as single pixels or large areas of missing values. These areas should be filled because they might otherwise generate discontinuous flow paths. Unless they are too large, I usually use the function inpaintnans to fill these gaps.

DEM = inpaintnans(DEM);

Having prepared the DEM, it is now time to derive flow directions. TopoToolbox stores flow directions in an object called FLOWobj.

FD = FLOWobj(DEM,'preprocess','carve');

FLOWobjs contain all necessary information to calculate hydrological terrain attributes (see the function help for details and some more detailed description on carving and filling here, here and here).

2. Choose a suitable threshold for stream initiation and derive a stream network

Now that you have a DEM and flow directions, you can derive a stream network. TopoToolbox stores stream networks in a STREAMobj. In essence, a STREAMobj is a reduced version of the FLOWobj being limited to the channelized domain. Unless you have the locations of mapped channelheads, the standard approach to determine the channelized domain is to assume a minimum upslope area (MUA). You should neither choose a too large nor too small MUA. A large MUA may hide chimap details close to catchment divides; a small MUA may extend from the channelized domain into the debris flow or hillslope domain. Here we simply take an upslope area of 0.05 sqkm or 500 pixels:

S = STREAMobj(FD,'minarea',500);

The resulting stream network is shown in the Fig. 1:

Fig. 1: Stream network

3. Make sure that all outlets have the same elevation

Chimaps require the stream network to have outlets with the same elevation. As DEMs often fail to cover the entire altitude range of all individual streams it may be necessary to trim the stream network in a way so that all relevant streams have the same outlet elevation. The function GRIDobj/griddedcontour and STREAMobj/modify come in handy here (see also this post). The modify function is extremely versatile. Here we use its option ‘upstreamto’ that requires a GRIDobj of zeros and ones. The modify function will then modify the stream network so that it contains streams only upstream of regions with ones. We use the output of griddedcontour and modify it slightly using the function bwmorph to produce gridded contours that are four-connected. A four-connected contour line ensures that all streams crossing the contours are reliably detected.

C = griddedcontour(DEM,[50 50]);
C.Z = bwmorph(C.Z,'diag');
S = modify(S,'upstreamto',C);

The trimmed stream network is shown in Fig. 2.

Fig. 2: The stream network reduced to those streams having an outlet at an altitude of 50 m a.s.l.

4. Make sure that all drainage basins are complete

Some of the drainage basins may be located along the DEM borders and have part of their area outside the DEM domain. To avoid edge effects on the resulting chimap, these incomplete drainage basins and their associated streams should be removed from further analysis. Testing whether a drainage basin is complete, however, is not trivial and the approach used here may not be useful in all situations.

Here I choose to identify incomplete basins by testing if their boundaries touch NaN areas or DEM edges. This requires some hard-coding since there is no ready-made function in TopoToolbox so far.

% Get drainage basins
D = drainagebasins(FD,S);

% Get NaN mask and dilate it by one pixel.
I = isnan(DEM);
I = dilate(I,ones(3));

% Add border pixels to the mask
I.Z([1 end],:) = true;
I.Z(:,[1 end]) = true;

% Get outlines for each basin
OUTLINES = false(DEM.size);
for r = 1:max(D);
OUTLINES = OUTLINES | bwperim(D.Z == r);

% Calculate the fraction that each outline shares with the NaN mask
frac = accumarray(D.Z(OUTLINES),I.Z(OUTLINES),[],@mean);

% Grid the fractions
FRAC.Z(:,:) = nan;
FRAC.Z(D.Z>0) = frac(D.Z(D.Z>0));

The resulting GRIDobj FRAC is shown in Fig. 3. Several basins share a large proportion of their outlines with the DEM boundary and we don’t want these to be part of the analysis.

Fig. 3: Drainage basins colored by the fraction that their boundaries touch missing values in the DEM.

Based on the GRIDobj FRAC we now modify the STREAMobj so that it contains only streams where drainage basins share less than 1% with the NaN mask. Again, we use the function modify:

S = modify(S,'upstreamto',FRAC<=0.01);

The trimmed stream network is shown below. Unfortunately, a large portion of the network was removed.

Fig. 4: Trimmed stream network (white) that excludes incomplete drainage basins.

Wrap up until here

As you see, preparing the DEM and stream network for chi analysis takes some effort. I may have been somewhat over-optimistic regarding the title of this series of posts. We may have many lines of code. However, having now a clean network will make the following steps much easier and straightforward. I will cover these in my next posts.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s