Extracting topographic swath profiles across curved geomorphic features

Abstract. We present a new method to extend the widely used geomorphic technique of swath profiles towards curved geomorphic structures such as river valleys. In contrast to the established method that hinges on stacking parallel cross sections, our approach does not refer to any individual profile lines, but uses the signed distance from a given baseline (for example, a valley floor) as the profile coordinate. The method can be implemented easily for arbitrary polygonal baselines and for rastered digital elevation models as well as for irregular point clouds such as laser scanner data. Furthermore it does not require any smoothness of the baseline and avoids over- and undersampling due to the curvature of the baseline. The versatility of the new method is illustrated by its application to topographic profiles across valleys, a large subduction zone, and the rim of an impact crater. Similarly to the ordinary swath profile method, the new method is not restricted to analyzing surface elevations themselves, but can aid the quantitative description of topography by analyzing other geomorphic features such as slope or local relief. It is even not constrained to geomorphic data, but can be applied to any two-dimensional data set such as temperature, precipitation or ages of rocks.


Introduction
Traditional cross-section analysis of topography always involves some randomness in the choice of the profile line.Swath profiles are a widespread tool to reduce this randomness and to focus on the morphologic features of interest.This technique has been used in numerous geomorphic studies.Some of them use swath profiles mainly for illustrating the topography and focus on profiles of surface elevation itself (e.g., Rehak et al., 2008;Robl et al., 2008), while several others also analyze further geomorphic properties such as slope and local relief (e.g., Fielding et al., 1994;Montgomery, 2001;Mitchell and Montgomery, 2006b).In some studies swath profiles were even used to examine the relationship between topography and other data that are potentially related to topography, for example precipitation (Fielding et al., 1994), ages of rocks (Reiners et al., 2003), and exhumation rates (Mitchell and Montgomery, 2006a).Several further examples of applying swath profiles were reviewed by Telbisz et al. (2013).As mentioned by Grohmann (2004), the basic idea even dates back to the 1920s.The implementation in geographic information systems, which was the main focus of the paper by Grohmann (2004), has also been addressed in several online tutorials and discussions.
In principle, creating a swath profile is nothing else but stacking several parallel profiles.As illustrated in Fig. 1a, computing a swath profile involves extending the original profile line (green) to a rectangle of a given width (the swath; black frame), and the elevation data are stacked along red lines normal to the profile line.In this context, stacking comprises computing at least a mean value of the elevations along the red lines, but in many cases, minimum and maximum values or the standard deviation are determined and analyzed, too.
Considering wider swaths reduces the effect of variations in topography normal to the profile line.This leads to smoother profiles and probably helps to recognize the characteristics of the along-profile topography.But in return, swath profiles are subject to a systematic bias as soon as the topography contains curved structures such as river valleys.Profiles across valleys are among the most widely used types of topographic profiles.Even if the original profile line is normal to the the valley floor, stacking along the red lines in Fig. 1a loses track of a curved valley floor depending on its curvature.As a consequence, even a perfectly V-shaped valley will look more like a U-shaped valley where the width of the artificially flattened valley floor increases with both the curvature of the valley and the width of the swath.This phenomenon is exemplified in Fig. 2 using a part of the Taugl Gorge in the Austrian Alps.This example will be discussed in more detail when demonstrating the capabilities of our new method in Sect. 4. We will find there (Fig. 4d) that the considered part of the gorge is very narrow, and that the walls become particularly steep towards the river.In Fig. 2 a rather wide (1.5 km) swath normal to a 1 km long profile across the valley is considered.Due to the nonstraight course of the river, the valley floor which is in fact very narrow appears to be about 100 m wide, and the steepness of the walls is underestimated.
In a recent paper, Telbisz et al. (2013) considered theoretical aspects of the swath profile method in detail.In this pa-per it was also shown that swath profiles are in principle not restricted to rectangular swaths, but also allow the consideration of curved structures.As an example, the application of "central swath profiles" to the morphology of a volcano was presented where the profile coordinate is the distance from a given point.Beyond this, even an extension towards curvilinear profile lines was suggested, and two different ways of stacking the data of the surface points not situated on the profile line were discussed.As an example, the application to a longitudinal profile over a mountain range was given.
In this paper we present an alternative idea of generalizing swath profiles towards curved geomorphic structures which differs from the idea of Telbisz et al. (2013) not only technically, but also in its spirit.While their method is designed to take a profile along (and in some sense parallel to) a given curved profile line, our approach generates swath profiles across a given (curved) baseline (for example, a river).The examples presented in Sect. 4 will illustrate that this concept has even a wider field of application in geomorphology than taking profiles along curved lines.

The generalized swath profile method
Our new approach hinges on the finding that the profile line of a topographic cross section is not a distinct morphological object in most cases.This particularly applies to swath profiles and is most obvious in the example of river valleys where the river is the distinct morphological feature.For a longitudinal river profile (i.e., a profile along the river line), swath profiles are not useful as they may be strongly biased by the cross-sectional shape of the valley even if the curvature of the river line is taken into account.In return, taking a swath profile makes sense if we are interested in the shape of the valley and choose a profile line across the valley.
In contrast to Telbisz et al. (2013) we therefore do not start from a distinct (potentially curved) profile line, but define a (curved) baseline (the river line for a profile across a river valley) that fixes the origin of the profile coordinate.The two more or less obvious ways to define a swath profile across (i.e., in some sense normal to) this baseline are discussed in the following.
The first way consists in defining several nonparallel profiles where each of them is normal to the baseline at different points on the baseline.The data of these profiles are stacked in such a way that the origin of each individual profile line is located at its intersection with the baseline.The green lines in Fig. 1b illustrate this concept.Although somewhat straightforward, this method still suffers from some problems.First, the density of the individual profiles varies within the considered domain, so that points close to the center of the baseline's curvature contribute more to the resulting swath profile than regions at the convex side.This results in an over-or undersampling depending on the curvature of the baseline.It is also easily recognized in Fig. 1b that the same surface point may even contribute to several of the stacked profiles, but at different profile coordinates.As a second problem, this method requires some smoothing of the baseline (in particular if it was automatically derived from a digital elevation model) since the orientation of the individual profile lines is very sensitive to the small-scale roughness of the baseline.
The alternative concept refrains from defining individual profiles, but directly uses the signed distance of any surface point from the baseline instead.Signed distance (sometimes called oriented distance) shall mean that all points to the right-hand side of the baseline have a positive distance, while those to the left-hand side have a negative distance.This signed distance defines the position of the considered surface point on the resulting generalized swath profile.In principle, this method stacks the elevation data along the red lines of constant signed distance in Fig. 1b.
Computing these lines of constant distance would, of course, be rather demanding.However, we do not need these lines explicitly, but simply travel through the entire domain instead, i.e., along all points of the DEM (digital elevation model) in the discrete case.We therefore only need to compute the signed distance of each DEM point towards the baseline of the profile, which is discussed in the next section in more detail.
Apart from the advantage of its simplicity, this procedure guarantees that each DEM point contributes exactly once to the profile, immediately avoiding over-or undersampling.Furthermore it is robust against the roughness of the baseline.It is easily recognized that the lines of small distances follow this roughness, while those of larger distances become smooth (see, for example, the red lines in Fig. 4a).
In view of these advantages over the method of stacking individual profile lines, we only consider the approach based on the signed distance in the following.

Implementation
The implementation of our method based on a DEM and a baseline defined by a polygonal line is in principle straightforward.In a first step, a bounding box, i.e., a tight rectangle with edges parallel to the coordinate axes around the baseline, is computed.This bounding box is extended by half the length of the desired profile in each direction, so that it includes all points of the surface that may contribute to the profile.The resulting rectangle is illustrated in Fig. 1b by the grey frame.
Then, all DEM points in this rectangle are iterated, and the signed distance of each point from the polygonal baseline is computed.As long as the number of line segments of the baseline is not too high, this can be done individually for each DEM point p without using further information.Let b 1 , . . ., b n be the points defining the baseline.We then determine the distance of p towards each segment, which is assumed to be a straight line between two adjacent baseline points b i and b i+1 .This means that we compute λ i for each segment i = 1, . . ., n − 1 in such a way that the distance between the point p and the point is minimal under the side condition 0 ≤ λ i ≤ 1.The minimum of all obtained distances yields the distance of the point p to the baseline.Practically, this algorithm is sufficient in all cases where the baseline is generated manually, so that its number of segments does not exceed some hundreds.A more efficient algorithm is only necessary for automatically generated baselines with a high resolution, as when flow channels derived from the DEM are considered.If the length of the baseline segments is similar to the mesh width of the DEM, the numerical effort increases like the baseline length raised to the power of three (or alternatively, to the DEM resolution raised to the power of three for a given baseline length), so that the procedure becomes expensive for long baselines or high DEM resolutions.
In this case, we suggest an iterative algorithm using available information on the nearest point on the baseline from the neighborhood of the considered point p.We then keep track of the index of the baseline segment that contains the presumably nearest point on the baseline for each DEM point p.When searching for the baseline point nearest to p, we consider the information on the nearest baseline segment in the neighborhood of p.Let s 1 , . . ., s 9 be the indices of the segments nearest to p among p and its eight direct and diagonal neighbors.We then search for the nearest point only in the baseline interval from segment min(s 1 , . . ., s 9 ) − 1 to max(s 1 , . . ., s 9 ) + 1.As a consequence, only a small part of the baseline has to be searched in each step.This iteration is repeated until the results do not change any more.In order to achieve a high efficiency and to avoid getting stuck at local minima, e.g., for meander-like baselines, the points of the DEM should be iterated in varying directions.A straightforward scheme would be row-wise forward, row-wise backward, column-wise forward, and then column-wise backward.Practically, only the variation between forward and backward iteration is important and provided a convergence after a few iterations in all tests that were performed.
The two end points of the baseline require a specific treatment anyway.As our approach is based on the distances towards the baseline, it extends the swath by semicircles around the end points of the baselines.These semicircles may introduce a strong bias, for example, for profiles across valleys.In order to avoid this artefact we include only these points of the DEM in the analysis where the nearest point on the baseline is none of the two end points of the baseline.The red lines in Fig. 1b already take this constraint into account.
As a last step of the analysis, the obtained pairs of elevation versus profile coordinate data must be aligned in bins.The bin width (on the profile-coordinate axis) defines the spatial resolution of the resulting profile.
Although the examples provided in Sect. 4 were derived from DEMs on regular lattices (either Cartesian coordinates or in longitude and latitude), the method can easily be applied to arbitrary point clouds without any loss of efficiency.Furthermore, it is not restricted to topographic data, but applicable to any two-dimensional data set, for example, hydrological data or the distribution of elements in a mineral.
An implementation of the method is available as an online tool at http://hergarten.at/geomorphology/swathprofile.php using data from several available DEMs with almost global coverage.Furthermore, a command line tool (C++ source code without requirement of any specific libraries) is included as a supplement, and the latest version can be directly obtained from the corresponding author.

Applications
In order to illustrate the new method for the construction of swath profiles, the applications presented in the following focus on topographic features that are inherently curved in plan view.Beyond the obvious application to valleys we present examples of a large subduction zone and an impact structure.

Fluvial and glacial valleys
The morphological analysis of valleys is probably the best example to illustrate the advantages of the method presented here.The cross-sectional shape of river valleys is often used to infer the incision process of the valley.In particular, two end members of valley shapes are usually discerned: Vshaped valley cross sections are considered to be characteristic for actively incising streams and a fluvial erosional history (e.g., Bull, 1979), while U-shaped cross sections are interpreted in terms of glacial erosion (e.g., Bonney, 1874).However, identifying these shapes quantitatively is often nontrivial: locations appropriate for representative individual cross sections are often difficult to find, for example, because of tributaries.And as discussed in the introduction, river valleys are typically curved, so that ordinary swath profiles suffer from artefacts that may even make a V-shaped valley look like a U-shaped valley.
In Fig. 3 we use the generalized swath profile method to characterize two valleys that are famous for their fluvial and glacial erosional history, respectively: a segment of the Colorado River upstream of the Grand Canyon, which is a typical actively incising bedrock channel (e.g., Pederson et al., 2006;Wernicke, 2011) and the glacially carved Yosemite Valley of California (Gutenberg et al., 1956).
For both generalized swath profiles the baseline was chosen along the actual course of the rivers draining the valleys.Our analysis reveals that both cross sections are qualitatively similar with a flat valley floor, steep faces, and well-marked ledges forming distinct shoulders on average.However, the apparently flat valley floor is just the water level for the Colorado River, while the Yosemite Valley is indeed much wider than the river itself.
Our next example (Fig. 4) refers to a V-shaped valley with a very narrow river and illustrates the power of morphological analysis of different river segments.The Taugl River is a tributary of the Salzach River that drains a major region of the Eastern Alps.It is characterized by several segments that flow through an ancient uplifted landscape, over a knickpoint in the channel profile, and in a young incised landscape.In the headwater region (Fig. 4b, c), the average cross section of the Taugl River is perfectly V-shaped with a narrow valley floor and an average topographic gradient of about 0.6 for the corresponding hillslopes.The central segment (Fig. 4d,  e) starts downstream of a distinct knickpoint (waterfall) and describes a narrow gorge with about 50 m high, almost vertical faces bordering the river.Above this very steep slope  (Gesch et al., 2002).
segment the topographic gradient declines with increasing distance to the stream forming convex hillslopes.When comparing this profile with the ordinary swath profile shown in Fig. 2, the advantage of the new method is obvious.The third segment (Fig. 4f, g) is similar to the middle segment except for an apparently flat region on the orographic left-hand side of the river.However, the standard deviations amounting to some tens of meters and the high local topographic gradients reveal that it is not a simple horizontal planation surface.Such a result may arise from the roughness of the surface, but may also occur if the overall slope of this surface parallel to the profile baseline differs from the slope of the valley itself.This uncertainty is, of course, an inevitable limitation when two-dimensional information is condensed into a onedimensional cross section.
In these profiles it is also clearly visible that the variability (both standard deviation and difference between minimum and maximum elevation) does not vanish at the baseline itself due to the downstream gradient of the river.In order to avoid this, our software optionally allows us to consider the elevations of all DEM points relative to their nearest point on the baseline and to transform the obtained profile (consisting of elevations relative to the baseline) back to absolute values afterwards.

Subduction zones and related oroclines
Subduction zones are a key to understanding plate tectonics (Stern, 2002).Their topographic and bathymetric representation often shows similar characteristic features over several hundred to some thousand kilometers along strike, but strong variations perpendicular to these zones.Abyssal plains, forebulge and the deep sea trench are located on the subducting lower plate while the related mountain ranges are located on the upper plate.In general, these geomorphic features are located at a certain distance to the subduction zone and profiles normal to the strike of these zones are useful to characterize the morphological features.However, small-amplitude features like the elastic forebulge on the subducting plate may easily be missed because they are swamped by irregularities along strike.The construction of swath profiles is therefore an essential aid in the topographic analysis, but may be hampered by the curvature of the subduction zone.
In order to illustrate the value of our new method, we present a topographic characterization of the curved plate boundary between the Nazca Plate and the South American Plate (e.g., Isacks, 1988) in Fig. 5.As the subduction of the Nazca Plate is clearly related to the topographic development of the Andes (e.g., Gephardt, 1994), the bathymetry and to-pography perpendicular to the plate boundary should exhibit characteristic features at similar distances to the plate boundary although the age of the Nazca Plate and its mechanical properties are spatially diverse (e.g., Capitanio et al., 2011).For our analysis we have chosen the deep sea trench between the two plates as the best feature to align the baseline of the general swath profile.
Even for the long baseline of more than 2000 km along strike used here, the standard deviations of the averaged elevations are quite small, especially for the submarine part of the swath.The analysis shows that the average surface elevation of the abyssal plains of the Nazca plate is −4.3 km with a standard deviation of less than 0.5 km.Local extrema in the curves of minimum and maximum elevation (light blue area) are predominantly seamounts (local maxima of the maximum elevation) and graben structures (local minima of the minimum elevation).The distribution of the seafloor topography is positively skewed due to the existence of some high seamounts that do not necessarily have a strong influence on the mean elevation and its standard deviation.In return, the distribution of the land surface elevations is negatively skewed because even the entire line of minimum elevation may arise from a single deep valley.
The forebulge at a distance of −50 to −100 km to the trench accounts for an average increase of about 0.2 km of surface elevation relative to the abyssal plains.The deep sea trench is on average 6.5 km deep with a standard deviation of about 0.7 km, indicating a very uniform depth distribution on the subsiding plate.On the side of the continental South American Plate, the surface elevation reaches sea level at an average distance of about 100 km from the trench.The standard deviation is small (approximately 0.7 km), indicating a uniform active continental margin over more than 2000 km in strike.The highest average surface elevation reaches about 4 km (with peaks above 6 km) and appears at a distance of 250 km from the trench.Here the standard deviation amounts to 0.8 km going along with the variable width of the Andes with a maximum west-east extent approximately at the maximum curvature of the trench.At a distance of 250 km and beyond the average surface elevation declines, while the standard deviation increases rapidly.Obviously, topographic build-up is still related to the subduction of the Nazca plate, but no longer uniformly distributed along strike.In summary, the swath profile constructed with our method shows that there are indeed characteristic topographic features perpendicular to the subduction zone and at similar distances to the trench.

Impact structures
In the analysis of impact structures, morphological features are essential to infer aspects of the impact process.The diameter and the depth of the crater are interpreted in terms of size and velocity of the impactor.The asymmetry of the crater is interpreted in terms of the impactor's shape and the obliquity Earth Surf.Dynam., 2, 97-104, 2014 www.earth-surf-dynam.net/2/97/2014/ of the impact angle (Littlefield and Dawson, 2006).Fig. 6 shows our topographic analysis of Meteor Crater in Arizona, also known as Barringer Crater.It was formed by the impact of the Canyon Diablo iron meteorite about 50,000 years ago in sedimentary rocks consisting predominantly of dolomites with minor sandstones (Barringer, 1905;Nishiizumi et al., 1991;Phillips et al., 1991).The crater shows a rounded rectangular shape and deviates from a circular or elliptical form due to preexisting joints.The baseline of our analysis follows the topographic high of the crater rim.About four-fifths of the entire crater contributes to the swath profile shown in Fig. 6b.The rest of the rim has been omitted just in order to illustrate the shape of the swath in this case.Interestingly, the surface elevation declines uniformly at increasing distance from the rim (base-

Conclusions
We have extended the concept of swath profiles towards curved baselines allowing the topographic analysis of river valleys, curved subduction zones, crater rims, and perhaps many other morphologic structures.Our method directly uses the signed distance of all surface points to the baseline.It immediately avoids artefacts of alternative approaches, namely the principal bias due to the curvature occurring in ordinary swath profiles and the over-and undersampling that occurs when individual, nonparallel profiles are stacked.Furthermore, the new method is robust against the roughness of the baseline and can be implemented efficiently in a numerical code.The high efficiency even persists when surfaces consisting of unstructured point clouds instead of DEMs on a regular lattice are considered.An implementation is available as an online tool

Figure 1 .
Figure 1.(a) Illustration of the idea of ordinary swath profiles.Green: profile line.Red: lines along which the elevation data are stacked.(b) Two ways of generalizing swath profiles towards curved baselines (blue).Green: profile lines normal to the baseline which might be stacked.Red: lines of constant signed distance from the baseline (the method proposed here).The gray frame defines the region of all points that potentially contribute to the profile as discussed in Sect.3.

Figure 2 .
Figure 2. (a) Ordinary swath profile across the central part of the Taugl Gorge (Austria).The curves in (b) comprise the mean elevations (black line), extreme values (hull of the light blue area) and mean values ±1 standard deviation (hull of the dark area).The data were derived from lidar based DEM data (spatial resolution 1 m) in the MGI (Military Geographical Institute) Austria M31 coordinate system provided by the federal government of Salzburg.

Figure 3 .
Figure 3. Fluvially versus glacially eroded valleys.The maps show sections of (a) the Colorado River formed by fluvial incision and (c) the Yosemite Valley shaped by glacial scouring.Constant distance lines are red, only the baseline is dashed red and white.The pale region covers all values included in the swath profiles shown in (b) and (d).Mean elevation (black line), extreme values (hull of the blue area) and mean elevation ±1 standard deviation (hull of the dark blue area) are plotted.Bin width is 10 m.Maps and swath profiles were derived from the National Elevation Dataset in the UTM (Universal Transverse Mercator) coordinate system with a spatial resolution of 10 m(Gesch et al., 2002).

Figure 4 .
Figure 4. (a) Average geomorphic parameters of three morphologically different sections of the Taugl Gorge (Austria).The baseline (white and red dashed line) follows the course of the Taugl River.The blue transparent regions cover all values included in the profiles as annotated.Topography (b, d, f) and topographic gradient (c, e, g) are shown, comprising mean values (black line), extreme values (hull of the light blue area) and mean values ±1 standard deviation (hull of the dark area).The bin width of all profiles is 5 m.Maps and profiles were derived from lidar based DEM data (spatial resolution 1 m) in the MGI Austria M31 coordinate system provided by the federal government of Salzburg.

Figure 5 .
Figure 5. Topographic and bathymetric representation of the converging Nazca and South American plates as an example of a curved subduction zone.The baseline of the generalized swath profile follows the trench (white and red dashed line).Red contours indicate the signed distance to the baseline.The pale region covers all values included in the profile.The mean elevation is shown by the black line.Extreme values and mean values ±1 standard deviation are indicated by the light-and dark-blue areas, respectively.Bin width is 1 km.Map and profile were derived from the ETOPO1 digital elevation model (http://www.ngdc.noaa.gov/mgg/global/).

Figure 6 .
Figure 6.Generalized swath profile analysis of Meteor Crater (Arizona).The baseline (white and red dashed line) follows the crater rim, and contours indicate the signed distance to the baseline.The blue transparent region covers all values included in the profile.Mean elevation, (thick black line), extreme values (hull of the blue area), and mean elevation ±1 standard deviation (hull of the dark area) are shown.Bin width is 1 m.The map and the profile are based on the lidar point cloud processed by the National Center for Airborne Laser Mapping (http://www.ncalm.org)and interpolated onto a regular Cartesian grid (UTM coordinate system) with a spatial resolution of 2 m.
line) as indicated by a very low standard deviation ranging from 1 m at the floor to 13 m at the steepest faces.