**International Journal of Geosciences**

Vol.07 No.02(2016), Article ID:63892,9 pages

10.4236/ijg.2016.72011

A Set of GRASS GIS-Based Shell Scripts for the Calculation and Graphical Display of the Main Morphometric Parameters of a River Channel

Aldo Clerici, Susanna Perego

Dipartimento di Fisica e Scienze della Terra “M. Melloni”, University of Parma, Parma, Italy

Copyright © 2016 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

Received 25 December 2015; accepted 23 February 2016; published 26 February 2016

ABSTRACT

For the analysis of river evolution, the use of quantitative parameters can be quite useful in order to assess changes in the channel planform. Among the several parameters proposed by different authors in a number of papers, channel length and width, braiding and sinuosity indexes, and channel lateral shifting are proved to be the most effective ones for a quantitative analysis of river changes. However, the calculation of these parameters is time-consuming, tedious and error- prone, even where made in a GIS environment. This work describes four shell scripts that perform fast and automatic calculation of the morphometric parameters and draw curves showing the variation of the calculated parameters along the entire channel development. The scripts are based on commands of the GRASS GIS free and open source software and, as input, they require a simple vector map containing the essential features of a river channel, i.e. bankfull channel limits and longitudinal and lateral bars.

**Keywords:**

GRASS GIS, Shell Script, River Morphometry

1. Introduction

In a recent paper [1] , the morphological evolution of the unconfined reach of the Taro River (Italian Northern Apennines) in the last two centuries (1828-2011) has been evaluated. The river features as of nine different dates (1828, 1881, 1958, 1976, 1999, 2003, 2006, 2008 and 2011) were surveyed based on old historical maps and recent orthophotos. The analysis of the morphological changes was conducted on a quantitative basis by comparing the values of the morphometric parameters commonly used to define the planform characteristics of a river unconfined reach, i.e. channel length and width, braiding and sinuosity indexes, and channel lateral shifting. As the calculation of these parameters requires a long and repetitive sequence of operations, even where carried out in a GIS environment, a set of shell scripts was created, mainly based on the commands of GRASS (Geographical Analysis Support System), a free and open source GIS [2] [3] . A first script defines, for the channel as of a specific date, the channel centreline, computes the channel length and width, and the braiding and sinuosity indexes. A second script evaluates the lateral shifting of the channel centreline between two dates.

Other two scripts construct graphics, in Postscript format, showing the curves of the values of the computed parameters along the entire channel centreline. More precisely, the first one draws the curves of width, braiding and sinuosity as of a specific date and shifting over a specific time interval. The second one draws the curves of a single parameter as of several dates in a single graph, allowing punctual evaluation of the parameter changes over time.

The use of these shell scripts allows fast and automatic calculation of the main morphological parameters and the drawing of graphs which detail the continuous parameters variation along the entire channel development. In this paper, the characteristics of the four implemented shell scripts are described.

2. Channel Centreline Definition and Parameters Calculation

2.1. Input Data

The input data for the first script consist in a vector map containing the bankfull channel limits and the limits of longitudinal and lateral bars. The area between the bankfull limits, with the exclusion of the areas occupied by bars, defines the active subchannels (Figure 1).

In order to have fixed reference points, in space and time, which the channel characteristics and changes from date to date refer to, a set of ground points along both borders of the channel must also be collected. These points, which can be digitized on a reliable map or acquired by GPS on the field, are projected perpendicularly to the channel centreline and reported along the x-axis of the graphics with their distances from one another and their cumulative distances from the origin.

2.2. Centreline Extraction

In most works, the channel centreline, i.e. the line whose points are equidistant from the two limits of the channel, is traced manually. In a few works, specific algorithms have been proposed to automatically derive this feature by thinning a channel raster map pixel by pixel [4] [5] or to extract it from the bank lines in vector form

Figure 1. Portion of the vector map showing the limits of the bankfull channel and of longitudinal and lateral bars, which is the input map for the calculation of morphometric parameters.

using Delaunay triangulation principle [6] . In [7] the centreline has been defined by joining a set of evenly- spaced points along the river channel, whose equidistance from the two banks has been determined by successive approximations. In GRASS, a specific module (v. centerline) is proposed to define a line that represents an approximation of the central tendency of a series of input lines, all of which with similar trajectories and, therefore, that can be used also to define the centreline of a channel represented by its two sides. However, the module gives good results when the lines are approximately parallel to each other. Moreover, it may give unreliable results when two reaches of the channel are very close to each other, as can happen at the neck of a meander.

In this work, a different approach is adopted and the centreline is extracted from the vector map containing the bankfull channel limits, manly by using the GRASS module that traces lines parallel to existing features. Starting from the bankfull channel margins, pairs of lines parallel to the channel banks are repeatedly traced out toward the channel centre. Since each of the two parallels traced out at each iteration is equidistant from the pertaining bank line, the intersection of the two parallels defines points that are equidistant from the banks (Figure 2(a1)-(a8)). The line joining these intersection points defines the channel centreline (Figure 2(b)). The length of the centreline is assumed as the length of the channel.

Figure 2. Procedure for channel centreline definition and channel transects tracing. (a1)-(a8) Iterative tracing toward the center of the channel of the lines parallel to the limits of the bankfull channel. (b) Tracing of the channel centreline by joining the intersections of the parallel lines. (c) Location of equally spaced points along the channel centreline. (d) Tracing of transects orthogonal to the centreline at the equally-spaced points.

2.3. Width Calculation

The channel width is commonly defined as the length of the line from bank to bank orthogonal to the channel centreline. In the above mentioned algorithms [4] [5] , channel width is automatically computed along the lines orthogonal to the pixels forming the computed raster centreline. In [8] , a Python script is described for the calculation of channel width at regularly-spaced transects orthogonal to the channel centreline. The same procedure is adopted in the script described in this work (Figure 2(c), Figure 2(d)).

2.4. Braiding Calculation

Many different indexes have been proposed by various authors to express the degree of braiding of a river [9] . Among them, the so-called Channel count index [10] is the most commonly used one, since it is simple to calculate and is the least sensitive to river-stage effects [11] . It is simply expressed as the number of active subchannels along a channel transect.

Similarly to what is reported in the above cited paper [8] , in the script described herein this index is computed by simply counting the number of active channels along the same transects where the channel width is defined.

2.5. Sinuosity Calculation

In its classical formulation, the sinuosity index is defined as the ratio of channel length to valley length [12] . It requires prior subdivision of the river valley into rectilinear segments. Then, the length of the channel centreline encompassed in the valley segment is divided by the length of the valley segment. The valley segmentation introduces a certain degree of subjectivity, especially where the valley shows a curvilinear path. To reduce the degree of subjectivity, in [8] a fixed length value is used for the valley segmentation. In this work, to further reduce subjectivity and obtain a more detailed assessment of sinuosity along the entire channel course, the river valley is ignored and only the channel centreline is used. More precisely, sinuosity is calculated by considering a portion of the channel centreline with fixed length, progressively shifted downstream by a constant distance. Sinuosity is obtained by dividing the fixed length of the centreline tract by the straight-line between its endpoints. The computed sinuosity value is assigned to the midpoint of the centreline tract (Figure 3).

2.6. Output

The output of the first script consists of:

1) A screenshot of three simple graphs which allow a first rough check of the results (Figures 4(a)-(c)).

2) A vector map containing the channel centreline.

Figure 3. Procedure for sinuosity calculation by dividing the fixed length of a centreline tract by the distance between its endpoints. The centreline tract is progressively shifted downstream by a fixed length along the entire centreline and the computed sinuosity value is assigned to the midpoint of the tract.

Figure 4. Example of screenshot of simple graphs of the computed width (a), braiding (b), sinuosity (c) and lateral shifting (d).

3) A report containing the main computed statistics (Table 1).

4) Three text files containing the values of width, braiding and sinuosity. In each file, the first column reports the progressive number of each measure, the second reports the distance from the origin and the third reports the parameter value (Table 2).

2.7. Lateral Shifting Calculation

The second script computes the channel centreline lateral shifting between two dates. As it requires the pre- emptive construction of the channel centrelines as of the two dates by the previously described shell script, a specific shell script has been constructed. The script, starting from a vector map containing the channel centrelines as of the two dates, computes the lateral shifting using a procedure similar to the previously described one adopted to define the centreline between the two channel banks. Indeed, the parallels to the two centrelines are traced out and their intersection joined by a line that represents the line equidistant from the two centrelines. The distances between the two centrelines are then computed along lines orthogonal to the computed axis at regular intervals and are assumed as the channel centreline shifting between the two dates (Figure 5).

The output consists of a report (Table 3) and of a text file with the shifting values at regular interval along the computed centreline axis (Table 2). A rough simple graph is also displayed on the screen (Figure 4(d)).

3. Graphs Construction

For the creation of graphs showing the variations of the parameters along the river channel, two scripts were built. The first one constructs graphs of the four calculated parameters (width, braiding, sinuosity and lateral shifting) as of a specific date. The input includes the vector map containing the channel centreline, the vector map containing the reference points and the four text files containing the parameters previously computed by the first script. In addition to the curves of the original computed values, curves of simple moving average of various orders can also be drawn. Obviously, vertical and horizontal scales, size and colours of the lines and characters can be defined by the user. Examples of the curves of the four parameters are given in Figure 6.

Table 1. Example of output report containing the main computed statistics.

^{a}Ratio of the channel area to the centreline length; ^{b}Mean of the width values computed at the transects orthogonal to the centerline; ^{c}Ratio of the centreline length to the distance between centreline endpoints; ^{d}Mean of the sinuosity values computed at the fixed length tracts along the centreline.

Table 2. Example of text files with the computed values of channel parameters (only the first and last 10 values of each file are reported).

^{a}Values are computed along transects 20 m-spaced (automatically recalculated at 19.9959 to ensure equal distance between points along the centreline); ^{b}Values are computed for centreline tracts 5000 m long and shifted 50 m downstream.

Table 3. Example of an output report resulting from the centreline lateral shifting calculation between two dates (the example refers to the time interval 1958-1976 in [1] ).

Figure 5. Procedure for lateral shifting calculation. The axis of the centreline position as of the two dates is calculated and the distance along the lines orthogonal to the axis is assumed as the shifting value.

(a)(b)(c)(d)

Figure 6. Examples of the curves of the four channel parameters: width (a), braiding (b), sinuosity (c) refer to the 1976 date, shifting (d) to the time interval 1958-1976 [1] . 57 reference points were collected. In addition to the curves of the original values, the moving average curves of 21 terms are traced (thicker red lines). In the lower x-axis the progressive distances (in meters) from the origin are reported. On the upper x-axis the positions of the reference points (numbered from R01 to R57) projected on the channel centreline are reported, together with their partial and cumulative distances. The parameter value is on the y-axis (for high resolution images, click: (a) http://www.cler.unipr.it/IJG/fig6a1.jpg; (b) http://www.cler.unipr.it/IJG/fig6b1.jpg; (c) http://www.cler.unipr.it/IJG/fig6c1.jpg; (d) http://www.cler.unipr.it/IJG/fig6d1.jpg).

The second script draws the curves of each single parameter as of several dates in a single graphic, allowing an overview of the progressive variations of each parameter over time. If many dates are considered and/or the differences among the curves are small, the graph proves easier to read and interpret if only the moving average curves are traced. An example showing moving average curves of the channel width as of several dates is reported in Figure 7.

The input and output of the four shell scripts are schematically summarized in Table 4.

4. Conclusion

The use of shell scripts based on the GRASS GIS commands allows fast and automatic calculation of the main morphological parameters and the drawing of graphs which detail the continuous parameters variation along the

Figure 7. Example of multiple dates graphic. Curves refer to the channel width as of the nine dates from 1821 to 2011 considered in [1] . For simplicity, the curves of the original width values have been omitted and only the moving average curves of 21 terms have been traced. In the lower x-axis, the progressive distances from the origin are reported. On the upper x-axes, the positions of 57 reference points (R1-57) projected onto each channel centreline are reported, together with their partial and cumulative distances. Since the channel length changes from date to date, the projections of the reference points change their positions along each x-axis. The difference in distance between each couple of reference points from one date to the previous one, which is a measure of the change in channel length, is reported on each axis. The different dates are marked with different colours, for both the curves and the axes of the upper part of the graph, proceeding from the oldest (1828), in red, to the most recent (2011), in black. The channel width value is on the y-axis. All values are in meters. For a high resolution image, click http://www.cler.unipr.it/IJG/fig7.jpg.

Table 4. Input and output for the four shell scripts used for the calculation and graphical representation of the channel characteristics.

entire channel development and a punctual evaluation of the river planimetric changes in space and time. Furthermore, the methodology here adopted reduces the possibility of errors and, most of all, makes the common and very subjective practice of prior channel segmentation unnecessary. The scripts here described simply require digitization of the bankfull channel limits, of the longitudinal and lateral bars and of a number of reference points.

Acknowledgements

The authors are grateful to D. Peis for his kind assistance in solving computer problems and to Interconsul s.r.l. for their valuable aid in translation.

Cite this paper

AldoClerici,SusannaPerego, (2016) A Set of GRASS GIS-Based Shell Scripts for the Calculation and Graphical Display of the Main Morphometric Parameters of a River Channel. *International Journal of Geosciences*,**07**,135-143. doi: 10.4236/ijg.2016.72011

References

- 1. Clerici, A., Perego, S., Chelli, S. and Tellini, C. (2015) Morphological Changes of the Floodplain Reach of the Taro River (Northern Italy) in the Last Two Centuries. Journal of Hydrology, 527, 1106-1122. http://dx.doi.org/10.1016/j.jhydrol.2015.05.063
- 2. GRASS Development Team (2015) GRASS: Geographic Resources Analysis Support System.

http://grass.osgeo.org - 3. Neteler, M. and Mitasova, H. (2008) Open Source GIS: A GRASS GIS Approach, 3rd Edition, Springer, US. http://dx.doi.org/10.1007/978-0-387-68574-8
- 4. Fisher, G.B., Bookhagen, B. and Amos, C.B. (2013) Channel Planform Geometry and Slopes from Freely Available High-Spatial Resolution Imagery and DEM Fusion: Implications for Channel Width Scalings, Erosion Proxies, and Fluvial Signatures in Tectonically Active Landscapes. Geomorphology, 194, 46-56. http://dx.doi.org/10.1016/j.geomorph.2013.04.011
- 5. Pavelsky, T.M. and Smith, L.C. (2008) RivWidth: A Software Tool for the Calculation of River Widths from Remotely Sensed Imagery. IEEE Geoscience and Remote Sensing Letters, 5, 70-73.
- 6. Haiyong, L. and Lihong, S. (2009) Double-Line River Axis Extraction Based on Delaunay Triangulation. In: Yaolin, L. and Xinming, T., Eds., Proceeding SPIE 7492, International Symposium on Spatial Analysis, Spatial-Temporal Data Modeling, and Data Mining, Wuhan, China, 13 October 2009.
- 7. Lauer, J.W. (2006) Channel Planform Statistics Toolbox. National Center for Earth-surface Dynamics. Minneapolis, MN 55414. http://www.nced.umn.edu/system/files/PlanformStatisticsTools_0.ppt
- 8. Casagrande, L., Cencetti, C., De Rosa, P., Fredduzzi, A., Martinelli, A. and Minelli, A. (2011) L’Utilizzo dei GFOSS nel Calcolo dell’Indice di Qualità Morfologica (IQM) dei Corsi d’Acqua. Geomatics Workbooks, 10, 57-78. http://geomatica.como.polimi.it/workbooks/n10/GW10-FOSS4Git_2011.pdf
- 9. Thorne, C.R. (1997) Channel Types and Morphological Classification. In: Thorne, C.R., Hey, R.D. and Newson, M.D., Eds., Applied Fluvial Geomorphology for River Engineering and Management, John Wiley & Sons Ltd., 175-222.
- 10. Ashmore, P. (1991) Channel Morphology and the Bed Load Pulses in Braided, Gravel-Bed Streams. Geografiska Annaler, 73A, 37-52. http://dx.doi.org/10.2307/521212
- 11. Egozi, R. and Ashmore, P. (2008) Defining and Measuring Braiding Intensity. Earth Surface Processes and Landforms, 33, 2121-2138. http://dx.doi.org/10.1002/esp.1658
- 12. Schumm, S.A. (1963) Sinuosity of Alluvial Rivers on the Great Plains. Geological Society of America Bulletin, 74, 1089-1100. http://dx.doi.org/10.1130/0016-7606(1963)74[1089:SOAROT]2.0.CO;2