Link to a pdf copy of the workflow:

Geomorphic indices are powerful tools with which to gauge spatial change across a landscape. Most of these indices are readily derived from topographic maps, whether in paper or digital form. Hypsometry, in particular, quantifies the distribution of relief across a given drainage basin. There are two primary metrics associated with hypsometry: the hypsometric curve (HC) and the hypsometric integral (HI).The hypsometric integral is a simple, scalar value representing the general relief of a basin. The hypsometric curve, in comparison, takes elevations distributed across a drainage basin, plots them against upslope contributing area, and produces a meaningful, distinct curve that can be compared with those of other basins.

If you have a paper topo map, fine. Mark out the drainage basin's boundaries. Highlight major elevation contours. For each contour, tabulate the elevation (z) and measure the area upslope (a). Measure the highest basin elevation (Z) and the total area (A). Plot z/Z on the y-axis and a/A on the x-axis. The shape of the curve can illustrate the present stage of landscape evolution, youthful, rugged slopes to mature flatlands.

If you're using ArcGIS, the CalHypso extension can do this automatically, using information from a DEM.

But suppose you don't have ArcGIS- you're using GRASS. How can you approximate the hypsometric curve?

1) Load a DEM

Obtaining elevation data is a preliminary stage in any GIS-based topographic analysis. Several good sources include the SRTM database (30 and 90-meter SAR data), OpenTopography (lidar), and the National Map (all ranges of data resolution and type, across the USA). Once the DEM is opened in a GIS, and in the correct projection, you’re ready to start.

2) Fill the DEM

DEM construction is a careful and involved procedure. Nonetheless, even the best DEM has artificial error. ‘Pits’, a particular elevation artefact in which a given pixel is lower than all adjacent pixels, interrupt most hydrological routines (which are necessary for gathering hypsometric data). Since water flows down slopes, any pit artefacts will mislead the algorithm and interrupt the true flow path. The solution to this problem is to “fill” the pits, allowing modeled flow to proceed in its real downhill course. This function is built into most hydrological toolkits. In GRASS, it’s called “r.fill.dir”. Use this tool to make a depression-less map.

3) Run the hydrological flow algorithm

Like pit-filling algorithms, flow modeling routines are common to most GIS softwares. In GRASS, this particular tool is called r.watershed. Its functionality far exceeds what we’ll use it for in this workflow. As one of the more powerful algorithms in GRASS, it can map out flow accumulation, direction of downslope flow, estimated locations of stream channels, watershed basins, etc. The only output we care about for this particular exercise, however, is the raster map of watershed boundaries.
In order to work across landscape scales, r.watershed incorporates a user-defined area threshold for calculated basins: all basins will be larger than a set number of pixels. Here, we can use that threshold to eliminate all the smaller catchments that we do not wish to study, and trim the map to primary drainages.

4) Use r.to.vect to make vector polygons (areas) of the watersheds.

5) You’ll probably still find that the number of output watersheds exceeds that which you wish to study. Right click on the watershed vector you’ve just created, in the table of contents, and open the attribute table. There, select all catchments you wish to ignore, and delete them.

6) Using univariate raster statistics, find the maximum and minimum elevation values in your DEM (which you’ve hopefully clipped to contain only your study area). These elevations will respectively have no upslope area and maximum upslope area, and represent the two extremes of your hypsometry plots. Subdivide these elevations equally, with sufficient sampling density to produce a curve that exhibits changes (10 or so subdivisions is a good place to start).

Here’s where raster algebra comes into play. If you have a DEM named ‘mapname’, for example, the following syntax:

‘mapname’ > 1000

will yield a new, binary map of all elevations greater than 1000 (whether it’s in meters or feet depends on your data). All pixels higher than 1000 receive a value of 1, and those lower receive a value of zero. This effectively creates a mask of your study area, whose shape is bounded by a minimum elevation you’ve chosen. Make binary maps for each of your determined elevation subdivisions. Name each raster map after the elevation it’s greater than.

7) Use the v.rast.stats tool to survey each binary map within the boundaries of your catchment vector. For each run, name the column after the corresponding raster name. For example, if you surveyed the raster map where all pixels = 1 above 1500’, name it f1500. Don’t begin the column name with a number; GRASS can’t understand that.

Now, you should have a shapefile (vectorized from the r.watershed output), with one row per catchment. Each column should contain statistics for the binary “elevation above” maps. Each column title should somehow contain, but not begin with, the elevation value for the given binary map.

The columns you’ll wind up keeping contain the sums of all pixel values within a given catchment. We don’t care about minimum, mean, or maximum raster values. Nor do we care about the variance or the number of pixels surveyed (pixels of zero value still count, so n is the same for each raster surveyed).

The sum, however, of all pixels within a given catchment captures exactly what we’re looking for. Add all the zero-value pixels (those below the elevation subdivision), and you still get zero. Add all the pixels with a value of one (those above the subdivision), and you effectively count the number of upslope pixels (which approximates drainage area). Now, the data is ready for Excel.

8) Export the catchment vector as a csv file, for use in Excel.

9) Open Excel. Rename the column headers, in which you’ve concealed the elevation subdivision value, to just the numeric elevation value. Delete all columns except the ones that sum the pixels within the catchment boundaries. Remember, since elevations below the benchmark are zero, and all elevations above are set to a value of one, summing these values can tell you exactly how many pixels are upslope of each elevation benchmark (a proxy for upslope drainage area).

10) Use the min() function to find the minimum value of each row (one row of elevation values, many rows of upslope pixel values- one for each catchment). Put the minimum values in a separate column.

11) Subtract the minimum value for each row from all elements those rows, using the $ sign to keep the cell containing minimum pixel count constant. This will yield a new set of rows, with the elevations and pixel sums shifted to a minimum of zero.

12) Use the max() function to find the new maximum of each row. Do the same as before, dividing each row by its maximum value.

13) Now, you should have one row of normalized (0-1 scale) elevations, and as many rows of normalized drainage area (approximated by number of upslope pixels) as there are catchments.

14) Plot normalized drainage area (x-axis) by normalized elevation (y-axis) for each basin. The resulting graph is a comparison of the hypsometric curve for each basin (Fig. 3).