Introduction

Knowledge of the direction of ice flow in glaciers and ice sheets is useful in a variety of applications. These include defining drainage regions and delineating flowbands in mass-balance and modelling studies, reconstructing complete surface velocities from satellite interferometric synthetic aperture radar (InSAR) measurements and making along-flow interpolations.

Flow stripes are curvilinear surface features of glaciers and ice streams formed by the advection of flow-induced surface undulations (e.g. Reference Glasser and GudmundssonGlasser and Gudmundsson, 2012). These features are typically, but not always, aligned parallel to the local ice-flow direction (Reference Hambrey and DowdeswellHambrey and Dowdeswell, 1994; Reference Raup, Scambos and HaranRaup and others, 2005). Because of the laminar nature of ice flow, individual flow stripes are persistent, ranging in length from tens to hundreds of kilometres (Reference Swithinbank, Brunk and SieversSwithinbank and others, 1988; Reference Casassa and BrecherCasassa and Brecher, 1993; Reference Gudmundsson, Raymond and BindschadlerGudmundsson and others, 1998; Reference JezekJezek, 1999; Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000a), and accordingly it is possible to extract ice-flow direction information from them. Typically this is done by tracing individual lineations from the imagery (e.g. Reference JezekJezek, 1999).

Two broad settings for flow-stripe formation have been described (Reference Gudmundsson, Raymond and BindschadlerGudmundsson and others, 1998; Reference Glasser and GudmundssonGlasser and Gudmundsson, 2012). First, they may form within a glacier or ice stream, due to ice flow over basal perturbations. Modelling by Reference Gudmundsson, Raymond and BindschadlerGudmundsson and others (1998) indicates how flow stripes can result from a flow-induced surface expression of changes in bedrock topography or spatial variations in the resistance to basal slip. Alternatively, flow stripes may form in convergent flow settings. Reference Glasser and GudmundssonGlasser and Gudmundsson (2012) provide a conceptual mechanism for flow-stripe formation where flow converges around nunataks or in the confluence zone of merging glaciers or ice streams. The nature of flow stripes differs according to their mode of formation: those produced in convergent flow situations tend to be more closely spaced than those formed by flow over basal perturbations (Reference Glasser and GudmundssonGlasser and Gudmundsson, 2012).

The amplitudes of surface undulations associated with flow stripes are typically 1–2 m, but they may be up to ∼10 m high (e.g. Reference Crabtree and DoakeCrabtree and Doake, 1980; Reference Casassa and BrecherCasassa and Brecher, 1993; Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000a; Reference Raup, Scambos and HaranRaup and others, 2005). Flow stripes are visible in C-band synthetic aperture radar (SAR) images, since they produce local variations in ice-sheet properties that influence the radar echo strength (Reference Casassa and BrecherCasassa and Brecher, 1993; Reference JezekJezek, 1999). Crevassing at the shear margins of ice streams and ice shelves also produces a strong radar echo and these areas appear bright in SAR imagery (Reference JezekJezek, 1999).

Here we describe a new, partially automated technique to determine the ice-flow-direction field from the RADARSAT-1 Antarctic Mapping Project (RAMP) mosaic SAR image of Antarctica (Reference JezekJezek and RAMP Product Team, 2002) which is based on detection of variations in image intensity that are associated with small-scale flow-induced lineations. The performance of the technique is evaluated by comparing flow-direction data determined over the Lambert Glacier– Amery Ice Shelf region within East Antarctica with an existing ice velocity direction field (Reference Young and HylandYoung and Hyland, 2002) and with directions that we derive from previously traced streak lines (Reference JezekJezek, 1999).

In accordance with standard fluid dynamics usage, such flow-induced lineations are streak lines, showing the cumulative history of the flow that has advected each segment of the stripe from its creation at the originating point, and only coincide with streamlines in the case of a stationary velocity field (Reference BatchelorBatchelor, 1967). The Steershead region of the Ross Ice Shelf (Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000a) and the StancombWills Ice Tongue of the Brunt Ice Shelf system (Reference Wuite and JezekWuite and Jezek, 2009) are two well-known areas where the streak lines and current velocities are not aligned. Our technique is useful for mapping in such situations, where one wants to compare the direction field of lineations with present-day velocities or modelled flow histories (Reference Hulbe and FahnestockHulbe and Fahnestock, 2004, Reference Hulbe and Fahnestock2007). Individual streak lines can be easily reconstructed using the field of orientation vectors.

Method

The orientation of image features can be calculated using a Radon transform (Reference RadonRadon, 1986; Reference CastlemanCastleman, 1996). The Radon transform encodes a spatial distribution (in our case a pre-processed version) of the image intensity, I(x, y), in terms of its values when integrated along each member of families of parallel lines through the image for all possible orientations, R(p±, θ), where θ defines the direction of the family of lines and p the perpendicular offset of each line from the origin (Fig. 1). Originally developed by Radon as an abstract mathematical quantity, it has gained prominence in recent decades as the foundation of X-ray tomography, where only the transformed quantity (i.e. variations in X-ray intensity detected using sets of parallel beams at a range of angles) is available and the computational task is to reconstruct the spatial pattern of density in the region scanned. We have deviated slightly from the usual Radon conventions in labelling the orientations by the angles of the lines rather than the directions of their normals (simply a rotation of 90°).

Fig. 1.

(a) An example image. (b, c) Slices through the Radon transform at (b) 30° and (c) 70° highlight the linear features orientated at 70° in the example image. (d) The full Radon transform shows four peaks at 70° at various offsets, corresponding to the four linear features in (a), while the uniform-strength constant-width bands correspond to the circles.

As each point, {p, B}, in the transform coordinate space represents a particular straight line in the image, peaks in the Radon transform (Fig. 1d) can be used to detect individual linear features in an image. Remote-sensing applications include the study of oceanic internal waves (Reference Ramos, Lund and GraberRamos and others, 2009), ancient Roman land boundaries (Reference BescobyBescoby, 2006) and tractor wheel marks (Reference Corbane, Baghdadi and ClairotteCorbane and others, 2011). We also start with an image, but our aim is to interrogate the Radon transforms of small sub-images to determine the dominant directions of lineations within them, viewing small enough regions for the flow features to approximate a collection of straight lines.

The steps involved in this process include: pre-processing the image, extraction of sub-images, application of the Radon transform, forming an aggregate measure of orientation, interpolation to find the dominant orientation of lineations, quality-based culling and, finally, post-processing. Parameters for these steps have been selected for use on the 125 m pixel RAMP Antarctic Mapping Mission-1 (AMM-1) SAR Image Mosaic of Antarctica Reference JezekJezek and RAMP Product Team, 2002), although the method is generally applicable to a wide variety of image types, possibly requiring some alteration to parameters.

Other landscapes, such as desert dune fields, formerly glaciated regions and geologically faulted landforms, also contain underlying directionality and could be amenable to the same technique, although usually individual lineations are less persistent than in glaciers or ice sheets, and selection of a suitable sub-image size might be challenging.

Image pre-processing

The image is despeckled by applying a median filter, with each point being replaced by the median of itself and its eight nearest neighbours. This reduces random noise, with minimal blurring of edges (Reference CastlemanCastleman, 1996). While the 125 m RAMP mosaic is derived from 25 m multi-look data and then resampled to 125 m resolution, which reduces the speckle noise from the underlying images, sufficient speckle noise remains to warrant additional despeckling (which otherwise impacts on the ability to discriminate between robust and less-reliable results, in particular the use of Eqn (4) to automatically cull results). The image is then isotropically edge-enhanced by the application of a Laplacian filter adhne, 2004). Specifically, the image is convolved with the matrix filter,

(1)

Finally, to improve the angular resolution of the Radon transform, the image is interpolated to double the number of pixels in each direction using the Lanczos2 rescaler, a windowed version of the ideal (and infinite) sinc rescaler, which offers a good compromise between aliasing, sharpness and ringing (Reference Turkowski and GlassnerTurkowski, 1990). These manipulations were carried out using bespoke Fortran routines, but are easily implemented within most image-processing software.

Radon transform

The Radon transform, R(p, θ), of a given sub-image consists of the collection of projections of the image intensity, I(x, y), along the lines with tangent vectors oriented at θ to the x-axis and offset by a perpendicular distance, p, from the origin (Reference CastlemanCastleman, 1996). It can be represented as an integration over the whole image by

(2)

where the Dirac delta function, δ(− x sin θ + y cos θ − p), restricts the two-dimensional integration to the appropriate straight line in the x-y plane. The range of the transform coordinates is 0 ≤ 0 < π and −∞ < p < ∞. In practice, p and the spatial integrals only range over the domain of the sub-image.

We are only interested in finding the orientation of the strongest features in a sub-image, rather than identifying a single line in I(x, y) or inverting recorded variations with p of the Radon transform for a set of θ directions, as in a typical X-ray tomography situation. Accordingly, as a measure of orientation signal averaged across the sub-image, we calculate the variance, σ ², in R(p, θ) for each 0 over the range of p sampled, i.e.

(3)

where R _D(p_i,θ) is the discrete sum corresponding to Eqn (2) and the extra factor of 11/n is for normalization. The variance, σ ²(θ), proved to be a good measure of the intensity contrast of multiple co-aligned linear features and is further enhanced by the Laplacian filter, which produces strong parallel positive and negative features for each edge. The value of θ corresponding to the maximum in σ ²(θ) is the relevant overall orientation angle for the sub-image.

High-resolution orientation

The angular resolution of the discrete Radon transform is limited by quantization of the sub-image into pixels. In particular, for too small an increment in Δθ, the Radon projection lines (at angles of θ and θ + Δθ and offset of p) will intersect the same pixels of the sub-image and therefore yield the same value for the discrete Radon transform, effectively discretizing the useful set of orientations.

To improve angular resolution, a sub-quantum angular resolution scheme (similar to the sub-pixel schemes used for image processing) is applied. The interpolated optimal orientation is obtained by fitting a parabola to σ ²(θ). The value of θ corresponding to the maximum value of this parabola, σ ² _max, is the interpolated orientation of the strongest linear features in the sub-image. Typical performance of this sub-quantum angular resolution scheme is shown in Table 1, and suggests that for a 46 pixel diameter Radon transform windowthe angular resolution scheme is unbiased and has a standard deviation of ∼0.2 times the selected angular resolution (Δθ ₂ in Table 1).

Table 1.

Performance of the sub-quantum angular resolution scheme applied to a 1024 × 1024 pixel auto-culled image of the Lambert Glacier region (Fig. 2). We compare the differences between measured orientations obtained for pairs of different angular resolutions, Δθ, using a 46 pixel diameter Radon transform window

Auto-culling

Provided the Radon transform has some variance with angle (i.e. σ ²(θ) ≠ constant), an estimated orientation of linear image features will be determined. As the output frequently includes a proportion of false matches, automatic culling is applied, based on the absolute magnitude of the maximum variance of the Radon transform and the relative magnitude of the maximum variance compared to the standard deviation of the variance, SD. In particular, the following three criteria for retaining points are applied,

(4)

(5)

(6)

Post-processing

As the Radon transform is only capable of discriminating angles within the range 0 ≤ θ < π, a π radian ambiguity exists in the flow direction determined from the orientation of linear flow features. To resolve this, the flow direction is compared with the slope direction of a smoothed digital elevation model (Reference Bamber, Gomez-Dans and GriggsBamber and others, 2009), and any points with a discrepancy exceeding π/2 radians have their orientation adjusted accordingly. This is a purely local treatment; more sophisticated methods to resolve the ambiguity for situations such as local reverse slope flow require examination of the regional context.

Finally, for convenience, we formatted the output of the orientation vectors to be compatible with the VELVIEW visualization tool that is distributed with the IMCORR software package (Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others, 1992), with the VELVIEW variables correlation strength, x error and y error being replaced by the left-hand sides of Eqns (4–6). This allows us to investigate the auto-culling parameter choices, and VELVIEW also permits easy manual culling of orientation vectors that are either associated with non-ice-flow features or are misaligned with the flow. A typical auto-culled orientation field that has been subsequently manually culled is shown in Figure 2.

Fig. 2.

(a) The Amery Ice Shelf and Lambert Glacier MODIS (Moderate Resolution Imaging Spectroradiometer) image, courtesy NASA Goddard Space Flight Center, Rapid Response. Also shown is the ASAID (Antarctic Surface Accumulation and Ice Discharge) grounding line (Reference BindschadlerBindschadler and others, 2011). The square box denotes the area used for testing the sub-quantum angular resolution scheme (Table 1), and the inset shows the location of the Amery Ice Shelf and Lambert Glacier system within Antarctica. (b) Flow-direction field from the RAMP mosaic (shown as background) for the 128 km × 128 km region bounded by the square box in (a). The auto-culled field has been further manually culled to remove correlations to non-flow direction features such as nunataks.

Evaluation

To evaluate the overall performance of the scheme we compared our flow directions with those from velocity vectors over the Amery Ice Shelf by Reference Young and HylandYoung and Hyland (2002) (Fig. 3). Over this region, we have 2670 coincident flow orientations, and the mean difference between the Radon transform orientations and those of Reference Young and HylandYoung and Hyland (2002) is −0.15°, with a relatively small overall standard deviation of 6.31°. However, this overall standard deviation does not capture the spatial pattern of deviations, as discussed below.

Fig. 3.

Difference between the flow-orientation fields of Reference Young and HylandYoung and Hyland (2002) and the Radon transform over the Amery Ice Shelf. Each square represents an individual Radon transform determination compared with co-located values from Reference Young and HylandYoung and Hyland (2002). The box indicates the location of Figure 4.

As the mean difference between the orientation fields is close to zero we assume it to be unbiased; however, the existence of regions, spanning several gridcells, with either a positive or negative difference from the Reference Young and HylandYoung and Hyland (2002) orientation field (Fig. 3) indicates there is some structure to the difference field and suggests a systematic component to the difference in orientations. This may be partially due to the difference in window sizes used when calculating the orientation fields: in determining ice velocities, Reference Young and HylandYoung and Hyland (2002) use a 0.34 km along-track by 1.14 km across-track window, but smooth with an 8 km window, whereas we have used a 3.875 km × 3.875 km window. The results presented in Figure 4, a sub-image extracted from the Figure 3 domain, demonstrate that the Radon orientations are better aligned with the streak lines than the directions from the Reference Young and HylandYoung and Hyland (2002) velocity vectors. The greater deviation from local flow direction of the flow-direction vectors of Reference Young and HylandYoung and Hyland (2002) in areas of strong streak-line curvature is likely to be due to the high degree of smoothing discussed above.

Fig. 4.

Local flow orientation in a region of high flow curvature on the Amery Ice Shelf (see Fig. 3 for location). The orientation vectors (black arrows) determined using the Radon transform conform closely to the visible streak lines. In regions of high curvature the velocity vectors (yellow lines) of Reference Young and HylandYoung and Hyland (2002) do not match the streak-line orientations as closely.

Returning to Figure 3, south of 72.5° S agreement between the flow-direction fields is better (the mean and standard deviation of this sub-region are −0.075° and 0.62°, respectively) where local changes in the flow direction are gentle. To the north of 72.5° S, at the eastern and western ice-shelf margins, greater differences between the flow-direction estimates are apparent, and these are often associated with strong local curvature in the surface features (e.g. Fig. 4 from the region of 71.5° S, 70.5° E). In more central regions of the ice shelf, differences between the orientation fields are due to both strong local curvature of the surface lineations and the presence of surface features that are neither indicative of nor related to the local flow direction.

Additional flow-direction data for the Amery Ice Shelf region can be obtained from hand-traced streak lines visible in RADARSAT-1 SAR images (Reference JezekJezek, 1999). We converted these data into a flow-orientation dataset by calculating the angle between adjacent points that were separated by at least 2 km on each traced streak line (flow stripe). The resultant orientation data were subsequently averaged onto a 2 km grid. Figure 5 shows the difference between 1233 common co-located orientation estimates from Reference JezekJezek (1999) and Radon transform derived datasets over the Amery Ice Shelf. The mean difference of −0.11 ° and standard deviation of 5.94° are similar to those determined from the above comparison of the Radon orientations with the velocity directions of Reference Young and HylandYoung and Hyland (2002). Once again, the greatest difference between our orientations and those derived from traced streak lines occurs in regions where local curvature of the streak lines is high (Fig. 5), indicating the Radon transform more effectively captures changes in the ice-flow direction over small spatial scales. This may be partly due to sampling involved in manual tracing, where the aim may be to keep track of an identified streak line rather than to accurately map its every turn.

Fig. 5.

Difference between the ice-flow-direction fields derived from Reference JezekJezek (1999) and from the Radon transform over the Amery Ice Shelf.

Discussion and Conclusion

An advantage of the Radon transform technique is that orientation information can be extracted from a single image. Clearly other image types, such as MODIS (Moderate Resolution Imaging Spectroradiometer) or Landsat optical imagery, can be used to obtain feature-orientation information. Figure 6 shows feature-orientation information extracted from a Landsat 7, band 8 panchromatic image from the southeastern region of the Amery Ice Shelf in the vicinity of Clemence Massif. To improve feature-orientation detection, histogram equalization was applied to increase the image contrast prior to application of the Radon scheme. During image pre-processing, the application of a Laplacian filter was found to amplify the speckle noise associated with images having a relatively small pixel size of 15 m. Better angular discrimination was obtained without applying the Laplacian filter, due to a sharper maximum in the Radon variance (Eqn (3)). Additionally, auto-culling was reduced to only applying Eqn (5), with the threshold value on the right-hand side increased to 0.5, as the other two tests (Eqns (4) and (6)) made no further improvement.

Fig. 6.

Flow-direction field for a Landsat 7, band 8 panchromatic image sub-region from the Clemence Massif region of the Amery Ice Shelf. Histogram equalization has been used to improve image contrast. (a) Radon transform window diameter 990 m or 96 pixels. (b) Radon transform window diameter 2010 m or 192 pixels.

As the Radon transform window diameter is a function of the image pixel size, the resolution of the image type chosen for feature-orientation detection is a primary factor in determining the maximum spatial and angular resolution at which feature orientations can be determined. Selection of the Radon transform window diameter requires finding a balance between angular resolution (which only depends on pixel count), the ability to detect lineations and sensitivity to localized changes in direction (which depends on the physical spatial scale of the lineations). Our approach requires a window sufficiently large to capture the persistent lineations. It must exceed the characteristic transverse dimensions sufficiently to give directional sensitivity. However, if the lineations curve appreciably within the window, they will contribute to a range of angles in the transform. At best, an average direction may be detected, but with reduced discrimination. This balance varies according to the image characteristics (e.g. for RAMP AMM-1 SAR images with a pixel size of 125 m the Radon transform window diameter of 46 pixels or 3875 m is a suitable choice). Results of processing Landsat 7, band 8 images with a pixel size of 15 m are shown in Figure 6 for window sizes of 96 (Fig. 6a) and 192 pixels (Fig. 6b), corresponding to 990 and 2010 m, respectively. Angular resolution is sufficient in both cases, but the larger window results in a sharper peak in the Radon variance (Eqn (3)) for persistent lineations, and hence better agreement with the underlying streak lines.

Certain ice-sheet surface features are not aligned with the local flow direction, yet create sufficient image contrast to generate orientation vectors during the Radon transform. Consequently, care must be taken to manually cull those orientation vectors with a high degree of misorientation relative to the local ice-flow direction that remain after auto-culling. Regions with large-scale surface rumples, crevassing or dune fields may require careful interpretation. Shadows could present an additional difficulty if applying this technique to optical imagery. In addition, when manually culling, care must be taken to identify regions where aligned features are indicative of surface wind phenomena rather than the ice-flow direction. Indeed, with our standard Radon transform window size and auto-culling parameters we do detect the orientation of Antarctic megadunes (Reference Fahnestock, Scambos, Shuman, Arthern, Winebrenner and KwokFahnestock and others, 2000b) in the RAMP mosaic, and with a suitable choice of imagery, sub-image size and spatial filtering our technique could be applied to detecting the orientation of other linear features. Furthermore, even finer-scale lineations (with pixel-scale transverse dimension) superposed on some megadunes are captured. These lineations are orthogonal to the dune fields and are believed to be aligned with the prevailing wind direction (Reference JezekJezek, 2008).

Because the Radon transform technique requires only a single image, it does not provide information on ice velocity magnitude, only its orientation. Despite this, the technique has several direct applications in the study of ice-sheet dynamics. Flow-orientation information can be used to refine ice-flux estimates in flux-balance simulations (Reference Wu and JezekWu and Jezek, 2004), especially for Lagrangian-type flux-balance codes (e.g. Reference Testut, Hurd, Coleman, Rémy and LegrésyTestut and others, 2003; Reference RobertsRoberts and others, 2011). In addition, flow-direction information is frequently required in determining surface ice-sheet velocities using InSAR techniques, where the flow direction is usually assumed to correspond to the direction of steepest descent of the ice-sheet surface, derived from digital elevation models. Flow-orientation information offers an alternative source of information for InSAR studies, and, furthermore, discrepancies between the measured flow direction and the direction of steepest descent in the ice-sheet surface indicate regions where the shallow-ice approximation is less applicable in dynamic modelling studies.

We have concentrated in this paper on presenting an objective method for detecting the orientations of lineations on the surface of the Antarctic ice sheet. This provides a means of mapping ice-flow directions, assuming steady flow and that persistent lineations typically represent streak lines (flow stripes). It might be argued that recent comprehensive mapping of ice-flow velocities (Reference Rignot, Mouginot and ScheuchlRignot and others, 2011) makes this unnecessary for Antarctica, but we consider our technique has significant uses. First, the analysis can be applied to single images from any epoch, generating fields of orientations – a more versatile attribute than manual traces of individual streak lines. Second, our scheme permits systematic comparison with present or past velocities to explore long-term steadiness in ice-flow direction. For example, the agreement between our orientations of lineations on the Amery Ice Shelf with contemporary ice-shelf velocities (Reference Young and HylandYoung and Hyland, 2002) and with traced streak lines (Reference JezekJezek, 1999) supports the view that the Amery Ice Shelf has had a steady flow regime over several decades (Reference King, Coleman, Morgan and HurdKing and others, 2007). This clearly can be extended to comparing directionality of streak lines from different epochs, even where local surface features are not sufficiently persistent to permit velocity determinations using feature tracking. Third, this approach may have value in other less well-mapped parts of the cryosphere, for example in permitting systematic mapping of flow fields for smaller ice fields and glaciers, potentially using aerial as well as satellite imagery. This could be useful in mass-balance studies to delineate drainage basins where surface velocities are sparse. Lastly, as we have suggested, by focusing on appropriate spatial scales, this technique could also be applied to mapping the orientation of other linear ice-sheet surface features not connected directly with ice-flow direction, such as megadune fields.

Encouraged by the success of the Radon transform technique presented here, we have processed the entire RADARSAT Antarctic mosaic, and manual culling is in progress. Applications for that flow-direction information include refining Antarctic mass-balance calculations, particularly for fluxes in fast-flowing outlet glaciers and ice streams. The comparison of streak-line orientations with ice velocity vectors for indications of unsteadiness in flow direction should also be particularly interesting.