Received 25 December 2015; accepted 23 February 2016; published 26 February 2016
In a recent paper  , 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   . 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   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  . In  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   , channel width is automatically computed along the lines orthogonal to the pixels forming the computed raster centreline. In  , 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  . Among them, the so-called Channel count index  is the most commonly used one, since it is simple to calculate and is the least sensitive to river-stage effects  . It is simply expressed as the number of active subchannels along a channel transect.
Similarly to what is reported in the above cited paper  , 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  . 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  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).
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.
aRatio of the channel area to the centreline length; bMean of the width values computed at the transects orthogonal to the centerline; cRatio of the centreline length to the distance between centreline endpoints; dMean 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).
aValues are computed along transects 20 m-spaced (automatically recalculated at 19.9959 to ensure equal distance between points along the centreline); bValues 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  ).
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.
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  . 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.
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  . 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.
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.