Hypsometric Curves, HI, R1 Index

Hypsometry and Beer: The causes of and the solutions to all my problems.
- Apologies to Homer Simpson


Hypsometry is calculated by watershed. There are a few ways to do it. You can create a binned histogram (area by bin, 100 bins) with the Reclassify tool. You can create a hypsometric curve using Jerry Davis’ Hypsometry Toolbox (cumulative area by bin, n = 100). You can do it, in part, using the R script by Veit Hoefler. You can create a map of HI values using a grid-based method provided by Shantamoy Guha.

An overview of hypsometry concepts is provided below Steps #1-3.

1a.) Reclassify DEM Into 100 Elevation Bins with an R Script
- If you use program R, the following script, below, by Veit Hoefler reclassifies DEM into elevation bins. If not an R user, skip down to the next step (“Begin with a projected DEM…”).
WESSO<-function(DGM,Bin){
min<-trunc(minValue(DGM))
max<-ceiling(maxValue(DGM))
DIFF<-max-min
class<-DIFF/Bin
break_v<-matrix(seq(min+class,max,class))
low<-cbind(seq(min,max-class,class))
numb<-seq(1,Bin,1)
rclmat<-cbind(low,break_v,numb)
rc<- reclassify(DGM, rclmat,include.lowest=T)
list<-list(a=rclmat,b=rc)
raster<-list$b
rclmat<-list$a
count<-freq(raster,useNA='no')
count_2<-hist(raster,breaks=Bin)
count_2<-c(count_2$counts)
count_2[1]
DIFFA<-as.numeric(count_2[1]-count[1,2])
laenge<-as.numeric(length(count[,2]))
count_3<-count[3:laenge,2]
DIFFA_2<-as.numeric(count[1,2])
as<-rbind(DIFFA_2,DIFFA)
nn<-matrix(rbind(count_2[2:length(count_2)]))
end<-rbind(as,nn)
end<-matrix(end)
info<-cbind(rclmat,end)
return(list(a=raster,b=info))}
</code>

1b.) Reclassify DEM Into 100 Elevation Bins with ArcGIS & Spatial Analyst
– Begin with a projected DEM clipped to the boundary of your study watershed (see Watershed Delineation).

– Convert your DEM from Float –> Integer format using Spatial Analyst > Math > Int tool.

– Use the Reclassify tool to create your elevation bins (n = 100)* from new integer DEM…

Spatial Analyst > Reclass > Reclassify tool
Input raster = watershed DEM
Reclass field = VALUE
Method = Equal interval
Number of Classes = 100

* Note: If you specify 100 bins in the tool, but only 10 bins are created, then you need to increase the maximum number of classes allowed using the "ArcMap Advanced Settings Utility.exe". Its easy to use. The maximum number of classes is 200. I found the utility here: C:\Program Files (x86)\ArcGIS\Desktop10.1\Utilities. Change settings under the "Raster" tab. You'll need to restart ArcGIS after making changes. Check the attribute table after restarting and running the tool; it should contain 100 records.

– An attribute table is available if the raster is in integer format (see Int tool).

– Open the table. The VALUE field contains the Bin ID number. The COUNT field contains pixel count for each bin.

– Optional: To graph the data in Arc: Table Options > Create Graph.

– To export the data to Excel: Table Options > Export to dBase (.dbf).

– Open the exported .dbf in Excel. Make sure to open Excel first, then find the file (show All Files *.*). Make sure to SHOW ALL FILE TYPES (don’t open the .xml file).


2.) Chart Binned Histogram in Excel

Example Dataset with Instructions & Formulas:
Example Binned Hypsometry Spreadsheet_Spring2015.xlsx

– Create a Combo Chart with both Vertical Column and Scatter types. In Excel 8, Insert menu > Recommended Charts > All Charts > Combo Chart…

Series 1 = Bin Area Bars (Column-type, Primary y-axis)
Series 2 = Cumulative Area Curve (Scatter-type, Secondary y-axis)
X-axis for both series = Bin Number

newhypso2

Example dataset created from a 30m DEM for a single watershed using the Reclassify tool in ArcGIS 10.1. Data is then exported in .dbf format and opened in Excel. Pixel count is converted to area (km2). Cumulative area is also calculated using equation formulas shown.

newhypsochart

Combo Chart: The hypsometric curve (RED DOTS) from Hypsometry Toolbox is plotted over the bin area histogram (BLUE COLUMNS) from Reclassify tool. Data is from table above.


3.) Chart a Hypsometric Curve with Hypsometry Toolbox (Jerry Davis)

– Download toolbox here: Hypsometric Tools Script by Jerry Davis
- Click yellow Download button.
- OK the license agreement.
– Download the folder toolbox/script, which comes as a folder named AS16830.zip.
– Unzip the folder to your Desktop.
– Open ArcToolbox > Right-click on the word “ArcToolbox” at top of frame > Add Toolbox > navigate to desktop > select the unzipped folder AS16830 > Hypsometry subfolder > Hypsometry.tbx. Click OK. Toolbox will be installed.
– Expand the “Hypsometry Tools: Generic Tools” toolbox > open the Hypsometric script:
Input = watershed DEM
Output = name of raster to be created
– Run tool.
– When complete (about 10 seconds), open the attribute table Options > Export table as .dbf (dBase Table).
– Do not add table to map.
– Open the new .dbf in Excel (show All Files *.*) and chart hypsometric curve series (scatter type) over the top of binned histogram (column type).
VALUE field = Bin number (X-axis).
CUMFREQ field = Cumulative Area (Y-axis).
– Format the Primary and Secondary Y-axes carefully.
– See figure for example watershed data above.

4.) Grid-based Method for Mapping Hypsometric Integral Values (Guha Shantamoy)
With this method, you begin with a DEM of your watershed and end with a spatially distributed map of HI values.
Word document: Gridded Hypsometric Integral using Arc Map 10_Guha.doxc


OVERVIEW OF HYPSOMETRY CONCEPT

Hypsometry is a measure of the relationship between elevation and area in a basin, watershed, or catchment (Langbein, 1947; Strahler, 1952). Basin hypsometry is strongly tied to flood response and the erosional maturity of a basin. The Hypsometric Index (a value) and Hypsometric Curve (a curve) are products of hypsometric analysis. Hypsometric Curves may be generated for either Basin Area (Hb) or Channel Lengths (Hr, Hn). If you do both, then you may wish to rise beyond your slothful nature and take the next step: calculate Demoulin’s R1 Index (Demoulin, 2011). The Spatially Explicit Hypsometric Integral is another option (Cohen et al., 2008), see bottom of post for references. Taken alone or together, hypsometric indices can shed light on the local effects of denudation and tectonic uplift. Best calculated with a cold beer in at least one hand.

THE HYPSOMETRIC INDEX
A hypsometric Index (HI), sometimes called the Elevation/Relief Ratio, can be quickly calculated for any basin. The HI is a number that coarsely summarizes a basin’s relief – a general index of erosional development. The HI is useful when comparing several basins to one another as it may highlight anomalous basins in a mountain range. HI values are reported to two decimals.

HIindex

Make sure you find the mean of elevation values correctly; its not the median. Click for larger image.

HI values of <0.30 describe “tectonically stable”, “denuded”, “mature” basins according to W.M. Davis (click to learn more about Davisian ideas on landscape development). HI values >0.60 indicate “unstable”, “actively uplifting”, “young” basins. Willgoose & Hancock (1998) have a slightly different take on HI. They consider values >0.5 dominated by diffusive processes (mainly hillslope processes). HI values <0.5 are considered dominated by fluvial erosion (channel processes play a larger role). More “balanced”, flattened S-shaped, or straight hypsometric curves (HI ~0.50), suggest a relatively stable, but still developing landscape. One should be cautious: it is possible for basins with very different geomorphic histories to have identical HI values. Similar-looking curves can be produced by complex interactions of climate, tectonism, sedimentation, rock resistance (Bishop et al., 2002), amongst other things, like badgers.

Find Total Relief in Simple Math
– Min and Max elevation values are quickly found in ArcMap by looking at the DEM layer in your Table of Contents. The Min and Max elevation values are shown beside the color ramp, if you’ve chosen the “List by Drawing Order” display button. If not, then look in Properties > Source tab > scroll down for Statistics.

– Subtract Min from Max to find Total Relief for the basin. Write this down.

Calculate the Hypsometric Index (HI) Value for a Basin

HI = (Emean – Emin) / (Emax – Emin)

Emean = Mean elevation value (from summary statistics for watershed raster; not median)
Emax = Maximum elevation value
Emin = Minimum elevation value (outlet)

– In regular guy language, the HI is just the mean incision of a basin (Emean – Emin) divided by the basin’s relief (Emax – Emin).

* Note: I was getting my hair cut in Laramie one time when a hobo shuffled into the shop and asked the barber if he was giving "regular guy" haircuts that day. By "regular" he meant "free". I like to think G4G is a regular guy website.


THE HYPSOMETRIC CURVE

For hypsometric curves, cumulative area and normalized relief are typically plotted on the X-axis and Y-axis, respectively. The X-axis is scaled from 0 to 1 (normalized, cumulative area). The Y-axis is normalized elevation, also scaled from 0 to 1. Several curves for different basins may be plotted together on one graph (add multiple Series to chart in Excel).

hypsometry Davisian concept

Davisian concept of landscape age expressed by channels in 3 hypsometric curves. Click for larger image.

The shape of a hypsometric curve is something you should learn to interpret at a glance. Shape is an indicator of dominant geomorphic processes at work in a watershed (diffusive or fluvial). A convex curve (“concave down” for those of you in Rio Linda, CA Logan, UT) indicates more of the watershed’s area (or volume of rock and soil) is held relatively high in the watershed. In this case, diffusive hillslope processes such as landsliding, rainsplash, interrill erosion, soil creep, etc., play a larger role. A concave curve indicates the bulk of the basin’s area (or volume of rock and soil) resides at relatively low elevation. More material has been removed from higher areas and either transported to lower areas or advected out of the basin completely. Concave curves indicate channelized/linear/fluvial/alluvial processes dominate.

ws

Example watershed and its channel network.

new 10 bins ws

DEM clipped to the watershed boundary and reclassified into equal-interval elevation bins. 10 bins makes for easy normalization. Each of the 10 bins is shown in a separate color. An area of each bin, a group of pixels, can be calculated quickly. Bin your data in ArcGIS before moving to Excel.

 

 

 

 

 

 

 

 

Example Hypsometric Spreadsheet Hb

Example Spreadsheet for hypsometric curve for a watershed (Hb). BLUE text is for AREA. RED text is for ELEVATION. Excel formulas are shown below table. Relevant values are shown above table. You can substitute Channel Network Length data to find curve Hn. Click on the figure to open a full size version. There’s an error in the formula below column F; use (X – Min)/(Max – Min) instead.

 

Spreadsheet Set-up and Calculations
– Set up a new spreadsheet in Excel with the columns and headings shown in the example spreadsheet figure.

– Clip DEM to watershed boundary.

– For the clipped DEM, open Properties > Classified > click Classification button…

– Set 10 classes using the Equal Interval method.

Example: The range of elevations in the example watershed is 1957 to 3187m. Break Values for Bins 1-10 are: 2080, 2203, 2326, 2449, 2572, 2695, 2818, 2941, 3064, 3187.

– With the Classification window open, click on each of the Break Values. The Number of Elements (pixels) in each bin shows up at the bottom of the window. Write these down.

– Enter into Excel spreadsheet under Break Values column for Bins 1-10.

– See Example Spreadsheet figure for hypsometry calculations for the area-elevation curve.

Charting Hypsometric Curve
- Chart below shows how to plot data from the example watershed spreadsheet.

Example Hypsometric Chart2

Hypsometric curve (red) for example watershed. Area by bin is plotted as a separate series on secondary axes (gray bars). The red curve is plotted on the primary axes (top, right). Click for larger image.

 

DEMOULIN’S R1 INDEX
** R1 Index information is currently being revised **

Uplift produces an erosion wave that propagates through a basin’s river network beginning at its outlet and proceeding to its headwaters. The drainage network responds (basin changes size) as the wave rolls through, more rapidly at first. This change is recorded in the morphology of the three hypsometric factors which comprise the R1 Index (Demoulin, 2011). The index does not constrain age of relative uplift (rock uplift or base level fall), Demoulin’s R1 has been shown to correlate with basin area. The index is designed to minimize influence of bedrock lithology.

Note: Please see new developments on the R/SR Index by one of today's leading researchers on morphometric analysis, Dr. Alain Demoulin at Universite de Liege (Demoulin, 2012; Demoulin et al., 2013).

Demoulin’s R1 Index Equation:
f(lnA) = Ir / Ib

R1 Index Fig Dec 2015

The 3 curves are defined in Excel by points on normalized axes. Add a best-fit curve to the points (probably a 2nd order polynomial with r-square values of >0.9). Use the lowest order polynomial you can for the trendline. Record the equation for the best-fit trendlines. The curves should not cross. Enter these equations for Hb, Hn, and Hr into WolframAlpha along with the normalized axis integral limits to find the areas Ib and Ir. Limits are 0 and 1. The ratio of Ir to Ib is the R1 Index for the watershed is in km2 units.

Click Here for PDF of the Figure Above

Demoulin’s R1 has three inputs: a.) the area hypsometric curve for the basin (Hb), b.) total length of the basin’s drainage network (Hn), and c.) the length of the long profile for the trunk stream (trunk stream is always shorter than total network length, thus Hn > Hr). You will need to find the equations for each curve and the areas between the curves (Calculus!). We’ll get around the calculus by letting Excel construct best-fit curves. We’ll copy those equations into WolframAlpha and let it do the integrals for us.

* Note: If you are comparing R1 values calculated for different basins, make sure you use a similar channel initiation threshold value for all. The threshold value controls the elevation channel heads, thus the lengths of the trunk stream and total network.

To begin, we need the X,Y points in Excel that define two of the length curves (Hb and Hn). X is normalized elevation, Y is normalized channel length (downstream distance).

Total Channel Network (Hb)
The basin’s stream network is created via steps described in the Watershed Delineation. Elevation bins and segment lengths were found in ArcGIS in earlier steps. Total network length is tabulated near the end of the Watershed Delineation lesson.

selectedchannel

Selected stream segments comprising the long profile of the trunk stream.

profilepoints

Trunk stream ready for intersection with contours and elevation extraction.

Trunk Stream Long Profile (Hr)
You need one continuous channel for the long profile (a polyline). Since there are often several tributaries in the headwater region to choose from, you’ll need to find the segment with the highest-elevation channel head. The tributary with the highest-elevation channel head must be part of Hr. Here’s how…

– Make a copy of the stream network layer in your Table of Contents. Rename the copy in ArcCatalog.

– Use the Editor tools to pare down the network in the new layer to a single channel. Delete the unwanted segments. When done, stop editing and save the edits.

– You should now have a polyline representing the trunk stream. The channel may be segmented, so use the Dissolve tool to merge them into one polyline (Data Mgmt > Generalization > Dissolve).

– Create a Contour shapefile from the clipped watershed DEM (Spatial Analyst > Surface > Contour tool). Use the equal interval value found previously for the 10 elevation bins (Ex: 123m interval for the example watershed).
– Name output shapefile = contour100m

– Create a point shapefile where the contours and the channel polyline intersect
Analysis > Overlay > Intersect tool
– Inputs = contours and trunk stream shapefile (both shapefiles)
– Name output shapefile = int_pts

– Return to the Editor toolbar, this time open an editing session on int_points. Arc will likely have created more than one point at each intersection. Zoom in to each crossing and delete the extra points. When done, stop editing and save edits.

10m contours, streamline for profile, and intersection points.

10m contours, streamline for profile, and intersection points. Make sure to clean up the points (see below).

– Next, extract elevations values for the intersection points from the DEM. Use the Extract Values to Points tool, the clipped DEM, and the intersection points as the inputs. The tool creates a new point shapefile with the extracted elevations in a field called RASTERVALU.
– Spatial Analyst Tools > Extraction > Extract Values to Points…
– Input point features = int_pts
– Input raster = watershedDEM
– Output point feature = int_points_extracted_elevs

Zoomed in view of points (red) created at intersections between Contours (black) and Streamline (green).

Zoomed in view of points (red) created at intersections between Contours (black) and Streamline (green).

– When complete, check to see if the Attribute Table of the output file (int_points_extracted_elevs) contains the extracted elevation values in a field called RASTERVALU. If you wish, copy these data to a new field with properly named Elev_m.

– Export the attribute as .dbf to Excel. Open Excel, File > Open, show All Files(not just all Excel Files), open the .dbf table. Save file as .xlsx.

– Normalize the long profile’s X and Y axes (see RI Index figure above).

– Plot curves for Hr, Hn, and Hb.

multiple intersections

Arc will often create several points (red dots) at locations where contours (black lines) and the streamline (green line) intersect. Clusters at the right and left in the figure are examples. Be sure to clean up the point shapefile by deleting the extras so there is just one point at each intersection. Use the Editor toolbar.

…more to come…

 3.) MatLab Method
See Wobus et al. (2007). FIlling DEM, burning in streams, then extracting the channel profile.

Find Best-fit Equations for Hn, Hb, Hr in Excel
1.) Start with your hypsometric curve charted in Excel.
2.) Add a Trendline curve that best fits the long profile (i.e., nth-order Polynomial or Power).
3.) In the Trendline options box, check the boxes for both ‘Display the equation on chart’ and ‘Display the R2 value on chart’. Write these down. You will need them later.

Note: It's best to use the lowest-order polynomial when selecting a best-fit curve. A stats person would be the one to ask why. My wife is one, but she won't tell me. Google 'parsimony' and 'curve-fitting‘.

4.) Repeat for the other 2 curves.
5.) Add curves to same chart as new Series.

Find Area Between 2 Curves
6.) Review the general method in the figure below.
7.) With the equations for all three best-fit curves in hand, you can move to WolframAlpha.
8.) Find roots if needed in Wolfram –> roots (f(x))-g(x))
9.) Enter in equations, beginning with “int”. See example syntax in figure below.

Finding Area Between Curves_R1 Index

General example for calculating the area between Line f(x) and a Curve g(x) between limits 0 and 6 (x-axis).
Click for larger image.


Using WolframAlpha Website
WolframAlpha is pretty cool. Exactly what the internet should be. The syntax for the example above –  finding the area between a line and a curve – is shown below. The line equation is y = 2x. The curve equation is y = x2-4x. We integrate over the range 0 to 6 on the x-axis. “int” prompts the program to calculate the integral. The solution, the area, is 36.

R1 example equation in WolframAlpha

Calculate R1 Index
With integral areas Ir and Ib as well as the integral values for Hb, Hn, Hr in hand, recall the equation:
R1 = Ir / Ib
Ir / Ib = (Hn – Hr) / (Hb – Hn)
(Hb – Hn) = Ib
(Hn – Hr) = Ir

Example:
Hb =
Hn =
Hr =
Ib =
Ir =

……….……….……….……….……….……….

Optional: Calculate the Spatially Explicit HI
Cohen et al. (2008) describe a method for calculating a pixel-by-pixel HI value. Their study examined HI patterns at the sub-watershed scale and compared them to 1-km resolution soil maps, geologic maps, and curvature values. The same form of the HI equation is kept, but the inputs change because you are calculating for each cell.

SEHI = Zlocalmean – Zpixel / Zmax – Zpixel

Zpixel = Elevation of a pixel
Convert DEM pixels to points. VALUE field in attribute table contains elevations.

Zmax = Maximum elevation of the watershed
This is a fixed value from watershed DEM (Properties > Source tab).

Zlocalmean = Mean elevation of contributing area for each pixel
Consider each pixel an outlet and the upstream area contributing flow to it its watershed. The local mean will therefor change with each pixel. Substitute the following into the main one, above.

Zlocalmean = ( Zmax / Zpixel ) / 2

……….……….……….……….……….……….


Refs:

Demoulin (2011) Geomorphology 126 <– R1 Index introduced by author here
Demoulin (2012) Geophysical Research Letters 39 <– New developments on R/SR method
Demoulin et al. (2013) Tectonophysics 608 <– New developments on R/SR method
WolframAlpha LINK
Stack Exchange – Hypsometric Curve LINK
Stack Exchange – Extract a portion of a DEM based on elevation LINK
Ecosystems Data Observatory – Notes on hypsometry LINK
Langbein (1947) USGS Water Supply Paper 968-C, p. 125-157 LINK <= first use of hypsometry?
Strahler (1952) GSA Bulletin 63, p. 1117-1142
Anderson & Anderson (2010) Geomorphology: The mechanics & chemistry of landscapes, 637 pgs.
ArcGIS Hypsometric Tools Script by Jerry Davis LINK; Here’ help for running python scripts LINK
Ascione et al. (2008)
Bishop et al. (2202)
Cohen et al. (2008) Journal of Geophysical Research 113
Gomez, L.H. (2008) MS Thesis on Gargathy Lagoon VA, Old Dominion University
Luo, W. (1998) Computers & Geosciences 24, p. 815-821
Oertel, G.F. (2001) Journal of Coastal Research 17, p. 775-783
Perez-Pena et al. (2009)
Pinter (1996)
Pike and Wilson (1971) <– Formula for HI, which they called E, or elevation-relief ratio
Weissel et al. (1994)
Davis Cycles of Erosion LINK
Willgoose and Hancock (1998) Earth Surface Processes and Landforms 23
Hancock (2005) Hydrological Processes 19
Hancock and Evans (2006) Earth Surface Processes and Landforms 31

Comments are closed.