diff --git a/examples/chimaps.ipynb b/examples/chimaps.ipynb new file mode 100644 index 00000000..68b7eb57 --- /dev/null +++ b/examples/chimaps.ipynb @@ -0,0 +1,482 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "481cc45e-82d8-47ba-b3e3-aefeb5dd9860", + "metadata": {}, + "source": [ + "# Chimaps in a few lines of code\n", + "\n", + "Based on the posts [Chimaps in a few lines of code](https://topotoolbox.wordpress.com/2017/08/18/chimaps-in-a-few-lines-of-code-final/) by Wolfgang Schwanghart (2017)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97c69d62-d347-4d9c-a1f4-53d8f26bf3c5", + "metadata": {}, + "outputs": [], + "source": [ + "import topotoolbox as topo\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "id": "30bae8dc-1834-44c6-9e9b-20154a172502", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "Since Willett’s et al. (2014) paper on the growth and decay of drainage basins, chimaps have become a popular tool to illustrate the direction of divide migration. Whether a drainage basin is a victim or aggressor, i.e., will loose or gain catchment area in the future, is indicated by the variable chi in the vicinity of the divides. chi (Χ) is a transformed horizontal coordinate measured from the outlet of a drainage basin. The transformation uses upslope area to linearize the concave upward profile for a well-chosen value of the mn-ratio. If all drainage basin outlets have the same elevation, then chi values should be roughly the same on both sides of the divides in a steady-state situation. In this case, no drainage basin reorganization should be expected. Yet, if chi values differ, then one should expect that divides migrate from the basins with lower chi values (i.e. the aggressors) towards those with higher chi values (i.e. the victims) (see Figure below).\n", + "\n", + "![](files/willett2014_fig1.png)\n", + "\n", + "Now how is it possible to calculate chimaps using TopoToolbox? It is actually pretty easy, but there is no ready-made function. And there is a reason for this. The description above shows that plotting chimaps involves several modular steps, each of which should be carefully done. These steps are:\n", + "\n", + "1. Load DEM and derive flow directions (FLOWobj).\n", + "2. Choose a suitable threshold for stream initiation and derive a stream network (STREAMobj).\n", + "3. Make sure that all outlets have the same elevation.\n", + "4. Make sure that all drainagebasins are complete, i.e. that the DEM covers their entire extent. Remove incomplete drainage basins.\n", + "5. Choose a suitable value for the mn-ratio.\n", + "6. Calculate the transformed variable chi (Hint: Don’t use chiplot, but check the low-level function STREAMobj/chitransform!).\n", + "7. Plot the results.\n", + "\n", + "\n", + "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.\n", + "\n", + "## Load DEM and derive flow directions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f78744d-60d0-413d-84d3-55050910d13c", + "metadata": {}, + "outputs": [], + "source": [ + "DEM = topo.read_tif('files/mendocino2.tif')" + ] + }, + { + "cell_type": "markdown", + "id": "08ec7d70-7831-4269-aae2-3d6b9a85fa6b", + "metadata": {}, + "source": [ + "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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cbef232c-ff3a-4a41-8922-99c6b75b798d", + "metadata": {}, + "outputs": [], + "source": [ + "# TODO: Show min/max elevation in info\n", + "DEM.info()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "749396f9-9759-45b3-9381-371ff9df8c23", + "metadata": {}, + "source": [ + "which provides several information including the minimum and maximum value. If these values are far beyond expected values, they should be replaced by NaNs.\n", + "\n", + "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.\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4087409e-e01b-4f8a-9c3c-b84af0969cd9", + "metadata": {}, + "outputs": [], + "source": [ + "# TODO: Implement inpaintnans\n", + "DEM.inpaintnans()" + ] + }, + { + "cell_type": "markdown", + "id": "cf83e596-5120-459b-94a2-8fe3aca7bc99", + "metadata": {}, + "source": [ + "Having prepared the DEM, it is now time to derive flow directions. TopoToolbox stores flow directions in an object called FLOWobj." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69b52f56-5e35-47fc-8f58-dff0ffaa86cd", + "metadata": {}, + "outputs": [], + "source": [ + "FD = topo.FlowObject(DEM)" + ] + }, + { + "cell_type": "markdown", + "id": "39d0e4b7-ec1d-46ce-9d24-8d635b5b989c", + "metadata": {}, + "source": [ + "FLOWobjs contain all necessary information to calculate hydrological terrain attributes (see the [function help](https://topotoolbox.github.io/pytopotoolbox/_autosummary/topotoolbox.FlowObject.html#topotoolbox.FlowObject) for details and some more detailed description on carving and filling [here](https://topotoolbox.wordpress.com/2014/04/26/the-difference-between-carving-and-filling-1/), [here](https://topotoolbox.wordpress.com/2014/05/04/the-difference-between-carving-and-filling-2/) and [here](https://topotoolbox.wordpress.com/2014/05/13/the-difference-between-carving-and-filling-3/)).\n", + "\n", + "## Choose a suitable threshold for stream initiation and derive a stream network\n", + "\n", + "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:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c53787e9-4816-4b42-8206-d36b856ae0ad", + "metadata": {}, + "outputs": [], + "source": [ + "S = topo.StreamObject(FD,units='km2',threshold=0.05)\n", + "# TODO: This is very hard to see. Hillshade?\n", + "S.show(overlay=DEM)" + ] + }, + { + "cell_type": "markdown", + "id": "7a1bb0ef-c1f8-4437-b4b8-208d0cead30b", + "metadata": {}, + "source": [ + "## Make sure all outlets have the same elevation\n", + "\n", + "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](https://topotoolbox.wordpress.com/2015/06/22/demo-on-stream-network-modification/)). 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](https://en.wikipedia.org/wiki/Pixel_connectivity#4-connected). A four-connected contour line ensures that all streams crossing the contours are reliably detected." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "364ce1d4-d796-4eb2-a91b-a6abb6e1524f", + "metadata": {}, + "outputs": [], + "source": [ + "# TODO: Implement these processing functions\n", + "C = DEM.griddedcontour(levels=[50],connectivity='4')\n", + "C.Z = bwmorph(C.Z,'diag')\n", + "S = S.modify(S,upstreamto=C)\n", + "C.show()\n", + "S.show(overlay=DEM)" + ] + }, + { + "cell_type": "markdown", + "id": "9aa01c9c-c0c9-48bc-bbd0-4ec2a8c69e92", + "metadata": {}, + "source": [ + "## Make sure that all drainage basins are complete\n", + "\n", + "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.\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4894225f-5248-4863-bc99-ec0fbf6608e8", + "metadata": {}, + "outputs": [], + "source": [ + "# Get drainage basins\n", + "# TODO: Implement drainage basins\n", + "D = FD.drainagebasins(S);\n", + " \n", + "# Get NaN mask and dilate it by one pixel.\n", + "I = np.isnan(DEM.z);\n", + "I = dilate(I,ones(3));\n", + " \n", + "# Add border pixels to the mask\n", + "I.z[[1,-1],:] = true;\n", + "I.z[:,[1,-1]] = true;\n", + " \n", + "# Get outlines for each basin\n", + "OUTLINES = np.zeros(DEM.shape,dtype=bool)\n", + "for r in 1..max(D):\n", + " OUTLINES = OUTLINES | bwperim(D.Z == r);\n", + " \n", + "# Calculate the fraction that each outline shares with the NaN mask\n", + "frac = accumarray(D.Z(OUTLINES),I.Z(OUTLINES),[],np.mean);\n", + " \n", + "# Grid the fractions\n", + "FRAC = GRIDobj(DEM);\n", + "FRAC.z[:,:] = nan;\n", + "FRAC.z[D.z>0] = frac[D.z[D.z>0]];\n", + "FRAC.show()" + ] + }, + { + "cell_type": "markdown", + "id": "93a48743-d43a-49b7-a998-cfac6b3033e1", + "metadata": {}, + "source": [ + "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.\n", + "\n", + "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:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4544b68f-b521-4ba2-8b37-d4ad281f6e8d", + "metadata": {}, + "outputs": [], + "source": [ + "S = S.modify(upstreamto=FRAC<=0.01);\n", + "S.show(overlay=DEM);" + ] + }, + { + "cell_type": "markdown", + "id": "ac1c7b7f-cdc8-413a-8d57-3252082bc5b9", + "metadata": {}, + "source": [ + "The trimmed stream network is shown above. Unfortunately, a large portion of the network was removed.\n", + "\n", + "## The mathematical basis\n", + "\n", + "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](http://onlinelibrary.wiley.com/doi/10.1002/esp.3462/abstract)). 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:\n", + "\n", + "$$\n", + "\\frac{\\text{d}z}{\\text{d}t} = U - kA^mS^n\n", + "$$\n", + "\n", + "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:\n", + "\n", + "$$\n", + "U - kA^mS^n = 0\n", + "$$\n", + "\n", + "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.\n", + "\n", + "$$\n", + "\\frac{\\text{d}z}{\\text{d}x} = \\left(\\frac{U}{kA(x)^m}\\right)^{1/n}\n", + "$$\n", + "\n", + "The trick is then to take the integral of ($dz/dx$) with respect to $x$ from a base level $z(x_0)$ (the elevation at the outlet). For convenience of units, we also introduce a reference area $A_0$.\n", + "\n", + "$$\n", + "\\int \\frac{\\text{d}z}{\\text{d}x}\\ \\text{d}x = z(x_0) + \\left(\\frac{U}{kA_0^m}\\right)^{1/n} \\int \\frac{A_0}{A(x)^{m/n}}\\ \\text{d}x\n", + "$$\n", + "\n", + "Now that looks complicated. Yet, by replacing the left-hand integral with z and the right-hand integral with\n", + "\n", + "$$\n", + "\\chi = \\int \\frac{A_0}{A(x)^{m/n}}\\ \\text{d}x\n", + "$$\n", + "\n", + "we can simply rewrite this equation as a function of a straight line:\n", + "\n", + "$$\n", + "z = z(x_0) + \\left(\\frac{U}{kA_0^m}\\right)^{1/n} \\chi\n", + "$$\n", + "\n", + "A straight line? Yes… but only if we choose the right m/n ratio.\n", + "\n", + "This was really a short derivation of chi. I can recommend the [paper of Perron and Royden (2013)](http://onlinelibrary.wiley.com/doi/10.1002/esp.3302/abstract) for a more comprehensive account.\n", + "\n", + "### References\n", + "\n", + "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.\n", + "\n", + "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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "a5cd86c9-bfac-4956-8e92-377398cce25f", + "metadata": {}, + "source": [ + "## Finding the right m/n ratio\n", + "\n", + "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. I stated that this model is a straight line if we choose the right m/n ratio, that is the concavity index. I’ll show how we can obtain this m/n ratio by means of nonlinear optimization using the function chiplot.\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "310df9f1-11b7-4122-835f-8eab3a690151", + "metadata": {}, + "outputs": [], + "source": [ + "# This won't work, but probably shouldn't...\n", + "flowpathapp(FD,DEM,S)" + ] + }, + { + "cell_type": "markdown", + "id": "910f501f-ef8b-4871-b6af-14dafe4bd757", + "metadata": {}, + "source": [ + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "07465c1f-3cf7-4d2d-88cc-95b2fc962ba5", + "metadata": {}, + "outputs": [], + "source": [ + "S.chiplot(DEM,A,mnplot=true,plot=false)" + ] + }, + { + "cell_type": "markdown", + "id": "8e6e3b92-fa0f-4cd5-80de-733e97691d20", + "metadata": {}, + "source": [ + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a96aafed-2138-4b8a-acee-e1db7fd42812", + "metadata": {}, + "outputs": [], + "source": [ + "C = S.chiplot(DEM,A);" + ] + }, + { + "cell_type": "markdown", + "id": "2fa5bf68-b559-4154-a1e7-8199ce5da806", + "metadata": {}, + "source": [ + "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:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3dc595a5-1055-4a4c-b5d3-a43e01f38a1d", + "metadata": {}, + "outputs": [], + "source": [ + "C" + ] + }, + { + "cell_type": "markdown", + "id": "b123ac99-81b5-4485-ae34-332f3b66472a", + "metadata": {}, + "source": [ + "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.\n", + "\n", + "Ok, we have the m/n ratio now. Next we will eventually produce the chimap that we initially set out for.\n", + "\n", + "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.\n", + "\n", + "## Plotting the chi-map\n", + "\n", + "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.\n", + "\n", + "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.\n", + "\n", + "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.\n", + "\n", + "Ok, let’s derive chi values for our stream network:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdcc71f0-cb55-497d-abf0-d66bf00e230a", + "metadata": {}, + "outputs": [], + "source": [ + "A = FD.flow_accumulation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6a828a6-7c24-4358-938b-bcdd2da21f50", + "metadata": {}, + "outputs": [], + "source": [ + "c = S.chitransform(A,mn=0.4776);" + ] + }, + { + "cell_type": "markdown", + "id": "9576cdfe-9252-4bfe-9509-abe05b144201", + "metadata": {}, + "source": [ + "In the next step, we will plot a color stream network on a grayscale hillshade:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ad7d2ac-3a01-45e5-ac1e-dc51391838ee", + "metadata": {}, + "outputs": [], + "source": [ + "S.show(overlay=DEM,color=c)" + ] + }, + { + "cell_type": "markdown", + "id": "de09ffe7-5ddb-4deb-a9dd-c1c52efb50c2", + "metadata": {}, + "source": [ + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27a77b44-5963-4071-bd24-271a95d40385", + "metadata": {}, + "outputs": [], + "source": [ + "# TODO: Zoom in to this location" + ] + }, + { + "cell_type": "markdown", + "id": "9c0b0b45-047c-44fd-a7cd-69b23f9a4b60", + "metadata": {}, + "source": [ + "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.\n", + "\n", + "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." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}