Morphometric analysis in Geographic Information Systems: applications of free software GRASS and R.

Development and interpretation of morphometric maps are important tools in studies related to neotectonics and geomorphology; Geographic Information Systems (GIS) allows speed and precision to this process, but applied methodology will vary according to available tools and degree of knowledge of each researcher about involved software. A methodology to integrate GIS and statistics in morphometric analysis is presented for the most usual morphometric parameters - hypsometry, slope, aspect, swath profiles, lineaments and drainage density, surface roughness, isobase and hydraulic gradient. The GIS used was the Geographic Resources Analysis Support System (GRASS-GIS), an open-source project that offers an integrated environment for raster and vector analysis, image processing and maps/graphics creation. Statistical analysis of parameters can be carried out on R, a system for statistical computation and graphics, through an interface with GRASS that allows raster maps and points files to be treated as variables for analysis. The basic element for deriving morphometric maps is the digital elevation model (DEM). It can be interpolated from scattered points or contours, either in raster or vector format; it is also possible to use DEMs from NASA's Shuttle Radar Topographic Mission, with 30m of ground resolution for the USA and 90m for other countries. Proposed


Introduction
Development and interpretation of morphometric maps are important tools in studies related to neotectonics and geomorphology, where the answers of natural landscapes to planet's interior dynamics is often masked by fast action of weathering, and the presence of drainage network anomalies and relief pattern discontiniuties may be related with recent terrain movements (Zuchiewicz, 1991;Rodriguez, 1993;Salvador and Riccomini, 1995;. Geographic Information Systems (GIS) allows this process to be fast and precise , but applied methodology depends on available tools and degree of knowledge of each researcher about involved software (Grohmann, 2004).
This work aims an evaluation of necessary proceedings to development and correlation of morphometric maps with integrated utilization of GRASS-GIS (U.S. Army CERL, 1993) and R statistical language (Ihaka and Gentleman, 1996), both distributed as free software under the GNU Public License. Examples will be presented from the Serra do Caraça Range region, eastern border of Quadrilátero Ferrífero, southeastern Brazil.

Geological and geomorphological context
The Quadrilátero Ferrífero region ( Fig.1), south São Francisco Craton, is characterized by Archaean granite-gneiss domed complexes coeval to Rio das Velhas greenstone belt (Machado et al., 1992), engaged to south and east in a Paleoproterozoic metamorphic belt (Mineiro belt, Teixeira and Figueiredo, 1991) defined as an acrescionary orogen of ca. 2.1 Ga. A dome-and-kell structure (Marshak et al., 1992), extensional, is admitted as result of orogen collapse after the collision of Minas passive margin and plate/microplate of Mineiro magmatic arc.
The study area presents four main geomorphological units ( Fig.2): Serra do Caraça Range, with average altitudes of 1400-1600m and maximum at 2064m; Serra do Pinho Range, in the eastern side of the area, with N-S trend; Chapada de Canga Plateau, in the central region of the study area, leveled at ca. 900m; Minas Gerais center-south and east Highlands, characterized by fluvial dissection. The scarp that limits the Serra do Caraça from the other units has hundreds of meters of fall, and leads to believe that not only erosional processes, but also post-cretacic tectonic movements contributed in the morphological evolution.
According to Varajão (1991), remains of planation surfaces in the Quadrilátero Ferrífero have close relations with lithostructural domains. Reactivations with vertical displacement of ancient faults would be responsible for present day altimetric differences (King, 1956;Barbosa, 1980). Evidences of Cenozoic tectonics in the Quadrilátero Ferrífero are observed in several sedimentary deposits and indicate three deformational events with distinct tension fields: the first event, extensive, oriented NNE-SSW and probably related with horsts and grabens oriented ESE-WNW; a second and more expressive event, mainly compressive NW-SE; the third event is considered to be the relaxing of previous event structures (Lipski et al., 2001).
The Quaternary Chapada de Canga Formation consists of a flat plateau formed by a succession of continental itabiritic ironstone pebble conglomerates and directly overlies the Eocene Fonseca Formation and the Precambrian basement (Sant'Anna and Schorscher, 1995). Cenozoic deposits are cut by NE and NW brittle faults and joints, related to reactivation of pre-existing structures in Precambrian basement, process that strongly influenced the development of present day landscape morphology and drainage network (Sant'Anna et al., 1997).

Methods
The GIS used was the Geographic Resources Analysis Support System (GRASS-GIS), running under Linux operating system. GRASS is an open-source project (U.S. Army CERL, 1993;Neteler, 1998;Neteler and Mitasova, 2002;GRASS Development Team, 2002), freely available on the Internet [1], which offers an integrated environment for raster and vector analysis, image processing and maps/graphics creation.
The database organization is related to "locations" and "mapsets"; the "location" involves all the work area, while the "mapset" is the portion active and used for analysis, which can be smaller or have the location's same size; several mapsets can be defined for the same location (Neteler, 1998).
The "region" is an important concept in GRASS and defines, inside the mapset, the area of interest and the spatial resolution for raster maps. The coordinates of the region's enveloping rectangle can be easily changed with g.region, a command used in several steps of morphometric analysis. In this paper, we worked on regions with 1m, 25m, 50m (default), 500m and 1000m of spatial resolution.
Statistical analysis of the evaluated parameters can be carried out on R, a system for statistical computation and graphics (Ihaka and Gentleman, 1996;Grunsky, 2002;R Development Core Team, 2003), through an interface with GRASS (Bivand, 2000) that allows raster maps and sites (points) files to be treated as variables for analysis.
There are several extensions to R, which adds to the base system alternative interpolation and analyze methods, such as kriging or akima splines. The R core package and extensions, as well as related documentation, can be obtained from CRAN (The Comprehensive R Archive Network [2]).
The concept of convex hull of interpolation is still of little use in the geosciences. According to Eddy (1977), the convex hull of a points (sites) dataset is the minimum area convex polygon that contains all the data, and represents the limits of spatial validity for the interpolating function. The s.hull command provides the convex hull for a sites file; this polygon can be converted to raster and used as a mask, thus limiting the interpolation.
The basic element for deriving morphometric maps is the digital elevation model (DEM). It can be interpolated from scattered points or contours, either in raster or vector format.
Sites interpolation can be done by inverse distance weighting with s.surf.idw, or by regularized splines with tension (Mitasova and Mitas, 1993;Mitasova and Rofierka, 1993), with s.surf.rst, command that allows the empirical adjust of mathematical parameters for tension and smoothness of lines. Vector contours are interpolated in the same way, with v.surf.rst. Interpolation from rasterized contours determines elevations using procedures similar to manual methods, with r.surf.contour.
Alternatively, one might use DEMs from NASA's Shuttle Radar Topographic Mission, with 30m of ground resolution for the USA and 90m for other countries [3].
Drainage network must be in vector format and two additional layers must be prepared, one with 2 nd and 3 rd -order streams for extraction of isobase lines and the second with 2 nd -order streams from its heads, for calculation of hydraulic gradient. All needed adjusts and classifications are performed on the vector edition module v.digit; v.extract command is then applied to create a new vector layer with the selected drainage orders only.

Lineament Analysis
For interpretation of morphostructural lineaments, shaded relief maps are generated with shade.rel.sh script and tracing is made in v.digit. This approach has advantages over the use of satellite imagery, since its possible to set up the position (azimuth and inclination) of scene illumination, thus emphasize the several existing lineament orientations, due the enhancement of directions perpendicular to lighting, in spite of parallels ones (Liu, 1987;Riccomini and Crósta, 1988;Grohmann andCampos Neto, 2002,2003).
The properties of orientation and length for each lineament were calculated with DXF_Lin, a software written by the author in Delphi language (Kylix environment) which reads a DXF file exported from GRASS and writes a text file that can be used as input in structural geology software.
In this work four shaded relief maps were used, with lighting 20° above the horizon at N90° N45°, N00°and N315°. Fig. (3) shows interpreted lineaments and rose diagrams for frequency and length, at 10° interval.

Hypsometry, slope and aspect
The hypsometric map is quickly obtained from reclassification of DEM with r.reclass; colors are adjusted with r.colors and contours are extracted in vector format with r.contour. DEM 3-D visualization (Fig. 4) can be done with d.3d or in NVIZ, where one can interactively adjust positioning of view, lighting and several others parameters. Slope (Fig. 5A) and aspect (Fig. 5B) are calculated with r.slope.aspect; slope can be expressed in degrees or percentage and aspect is counted counterclockwise from east (N90º) in Cartesian angles, by default. To adjust aspect values to compass angles, one can use a script (Appendix A) written by M. Hamish Bowman (Department of Marine Sciences, University of Otago, New Zealand).
In Fig. (5B), the three highlighted subsets represent regions used for analysis of slope as function of aspect. First, subregions are created with g.region; then, r.resample is used to create raster layers with data only inside the subregions. Density plots can then be made to show aspect variation in each subregion (Fig. 6A).
Aspect map can also be reclassified and crossed with slope. In this work, aspect was classified into four categories: NE (0-90º), SE (90-180º), SW (180-270º) and NW (270-360º), in separate layers. Multiplying these layers with slope, using r.mapcalc, provides a way to analyze slope variations in each aspect category, as in Fig. (6B, C, D).

Swath profiles
Altimetric swath profiles (or projected profiles of Baulig, 1926, cited in Tricart andCailleux, 1957) are those were intersections of contours with equally spaced profile lines are marked within a swath, or band. This kind of profile can provide a broader view of altimetric behaviour, and help to determine inclination of large topographic features in planaltic regions (Meis, 1982). The d.profile command outputs up to four ASCII files with cumulative length in meters, raster value (elevation), east and north coordinates; these files can be imported into common spreadsheet software for profile construction. In this work, three swath profiles were made, with 500m of distance between each profile line (Fig. 7).

Drainage and lineament density
Drainage density maps (Horton, 1945) express the ratio between accumulated length per area (km/km²) and was calculated for 1000x1000m cells (Fig. 8A). The same specifications were used in the lineament density map (Fig. 8B).
To determine drainage and lineaments densities, one must convert the vector maps to raster with high ground resolution (g.region), to ensure details preservation. This map is reclassified, to assign all the lines a unitary value, and then converted to sites.
The resolution is changed for that chosen to calculate density (e.g. 1000m) and s.cellstats is used with the "sum" option to count how many points are in each cell. The result is converted to raster and multiplied, with r.mapcalc, for a value that represents the ratio of the cell size in high resolution and the area of low-resolution cell, to achieve the desired relation of length/area, in km/km². For example, if the vector map was converted to raster with 25m resolution, and one searches density for 1000m cells: density = (sites sum) * (25m/1000) (1) or, simply density = (sites sum) * 0.025 (2) This map is finally converted to sites; the resolution is turned back to default value and the sites are interpolated to achieve a continuous surface. The process is summarized in Fig. (9).

Surface roughness
The surface roughness map (Hobson, 1972;Day, 1979) was calculated as the ratio between the surface (real) area and the planar (flat) area of cells with 1000x1000m (Fig. 10A, 11). This method is useful for morphological division since it considers the shape and not the altitude. Thus, tilted relives have their expression showed, while it could be masked in a hypsometric map, as consequence of altimetric variations (Grohmann andCampos Neto, 2002, 2003).
The surface area of cells is gotten from the geometric relations between the regular grid and slope (Fig. 10B).
The map for cells planar area can be achieved reclassifying any preexisting raster layer with a unique value.
Both maps are then converted to sites; grid resolution is adjusted to the desired value (e.g. 1000m) and s.cellstats is used with the "sum" option to get the surface and planar areas; the results are converted to raster and ratio is calculated in r.mapcalc. This result is converted to sites and interpolated with default resolution.
The map of isobase was made from the intersections of contours with 2 nd and 3 rd -order stream channels (drainage orders according to Strahler, 1952a,b).
First, vector contours and streams are converted to raster, to produce maps with values 1 (one) and null for the streams, and altitude and null for topography. This allows that maps, when multiplied (r.mapcalc), present result only where there is an intersection of features (Fig. 12). The resulting raster is converted to sites and interpolated (Fig. 13A). Additionally, manual interpolation was made (Fig. 13B) and maps compared with linear regression (Fig. 13C).

Hydraulic Gradient
In an attempt of determine areas with similar hydraulic behavior, a map of hydraulic gradient was elaborated (Fig. 14), according to the proposition of Rodriguez (1993). This parameter is evaluated for each 2 nd -order stream as the ratio of the altimetric difference between head and mouth with the plan length; the value is attributed to the mid-portion of the stream.
In this work, an automation of this process was achieved with Grad_Hidro, a software written in Delphi language (Kylix environment), which relates the spatial position of drainage vector lines with altitude in DEM (Fig. 15).
Initially, the vector layer containing only 2 nd -order streams is exported in ASCII format. The program reads this file and returns a sites file with points corresponding to vertices located at the extremities and middle of each line, which must be imported into the system and converted to raster with 1m resolution. This raster is multiplied with the DEM, to obtain a layer with information only where the drainages vertices are; this is made to save disk space and processing time in the following steps. The result of the multiplication is converted to sites and the program is executed again; it assigns DEM´s altitude to vertices, calculate the hydraulic gradient and attribute the value to the mid-point of each stream. The resultant sites file is imported into GRASS and interpolated with default resolution.

Discussion
Morphostructural lineaments analysis indicate the presence of main trends oriented NE-SW and NNE-SSW through the whole area, with secondary trends oriented N-S, NNW-SSE and NW-SE (Fig. 3). In the northern subregion, NE-facing slopes prevail, and have values ranging mainly from 0-10º; the opposite direction, SW, have more slopes within 10-20º, which leads to asymmetric ridges oriented NW-SE with steepest slope in the SW face (Fig. 6B); low intensity of slopes can be seen as indicative of a very eroded area. The western subregion, which comprises essentially the Serra do Caraça Range, can be divided in north and south; to the north, slopes face mainly NE, with lower dips than those facing NW and SW; to the south, majority of slopes face SW (Fig. 6C). The southern subregion has lined ridges oriented mostly NE-SW, but the NW-SE trend is also present. Drainage (Fig. 8A) and lineament (Fig. 8B) densities agree with general orientation of these structures. In general, the south portion of the entire area has fewer lineaments than the north portion; in the center of the area, over the Chapada de Canga Plateau, there is a lack of both drainage and lineaments. In the southern portion of Serra do Caraça Range, low drainage density and high lineament density is associated with a subterranean drainage system, responsible for an extensive cave network development. These caves, such as the Centenário Cave, the deepest cave of the world in quartzite, are developed in WNW-SSE vertical faults (Dutra et al., 2002).
Swath profiles (Fig. 7) allowed identification of altimetric leveling at ca. 2050m (profile 2) and ca. 1650-1700m (profiles 1 and 2) in the Serra do Caraça Range. A step in the Chapada de Canga Plateau between 900 and 850m (profile 3) is related to a E-W normal fault with south block down; this fault can also be seen as an anomaly in orientation of isobase lines, as pointed out in previous works (Sant'Anna, 1994;Sant'Anna et al., 1997).
Surface roughness (Fig. 10A) and hydraulic gradient (Fig. 14) show a similar schema, with a relatively flat area in the central region, surrounded by areas of higher roughness and gradient. In both cases, the Serra do Caraça Range represents a very fast change of regional trend and has the highest values of all area. This trend change is also present in the isobase map (Fig. 13A), associated with others anomalies, such as alignment and approach/separation of isobase lines.
Although an important morphological feature, the Chapada de Canga Plateau is not well distinguished in hydraulic gradient and isobase maps, since these parameters are based on the relations between topography and drainage network, and the Plateau has very low drainage density.
Linear correlation models were applied to all evaluated morphometric parameters. The best results were achieved among hydraulic gradient and surface roughness (Fig. 16). Despite the observed dispersion, the graph shows a positive linear correlation between parameters.

Conclusions
Within the study of recent evolution of the eastern border of the Quadrilátero Ferrífero region, several morphometric methods were applied in a GIS environment, and proposed methodology is presented.
The integration between GIS and statistics in morphometric analysis allows agility and precision in determination of necessary parameters. The methods presented in this work can be adapted according to necessities and available tools; documentation of input and output file formats in the source code makes easy to implement custom-written software. Use of free and open-source tools guarantees access to everyone, and its increasing popularization opens new development perspectives in this research field. GIS mailing list, in special Hamish Bowman, Glyn Clements and Markus Neteler; Roger Bivand for reviewing the manuscript.