Inferring ice-flow directions from single ice-sheet surface images using the Radon transform

Abstract

We present a new method for extracting the direction of surface flow for ice sheets, based on the detection of flow-induced features that are visible in satellite imagery. The orientation of linear features is determined using a Radon transform and only requires a single image. The technique is demonstrated by applying it to the RADARSAT mosaic of Antarctica, over the Lambert Glacier–Amery Ice Shelf region of East Antarctica. Comparisons with both existing flow-direction fields and traced streamlines over the same area provide an evaluation of the method. We also illustrate its application to Landsat 7 imagery.

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. Glasser and Gudmundsson, 2012). These features are typically, but not always, aligned parallel to the local ice-flow direction (Hambrey and Dowdeswell, 1994; Raup 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 (Swithinbank and others, 1988; Casassa and Brecher, 1993; Gudmundsson and others, 1998; Jezek, 1999; Fahnestock 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. Jezek, 1999).

Two broad settings for flow-stripe formation have been described (Gudmundsson and others, 1998; Glasser and Gudmundsson, 2012). First, they may form within a glacier or ice stream, due to ice flow over basal perturbations. Modelling by Gudmundsson 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. Glasser 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 (Glasser 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. Crabtree and Doake, 1980; Casassa and Brecher, 1993; Fahnestock and others, 2000a; Raup 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 (Casassa and Brecher, 1993; Jezek, 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 (Jezek, 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 (Jezek 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 (Young and Hyland, 2002) and with directions that we derive from previously traced streak lines (Jezek, 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 (Batchelor, 1967). The Steershead region of the Ross Ice Shelf (Fahnestock and others, 2000a) and the StancombWills Ice Tongue of the Brunt Ice Shelf system (Wuite 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 (Hulbe and Fahnestock, 2004, 2007). 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 (Radon, 1986; Castleman, 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 (Ramos and others, 2009), ancient Roman land boundaries (Bescoby, 2006) and tractor wheel marks (Corbane 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 Jezek 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 (Castleman, 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 (Turkowski, 1990). These manipulations were carried out using bespoke Fortran routines, but are easily implemented within most image-processing software.

Sub-imaging

A compromise must be made when selecting the sub-image size for the Radon transform: larger images improve the angular resolution of orientations of lineations, but adversely affect how well local changes in orientation can be captured. To ensure a uniform number of pixels are used in the Radon transform, irrespective of the orientation of the parallel lines of integration relative to the array of pixels, we take the maximal n × n pixel sub-image inscribed in a circle of diameter w. Specifically, the Radon transform window size is n × n and , where floor(η) returns the largest integer smaller than or equal to η. The actual n × n set of pixels used will vary with the orientation under consideration.

For the RAMP AMM-1 mosaic, with 125 m pixel size, a sub-image diameter of 46 pixels, with a corresponding window size of 3875 m, results in good levels of both angular resolution and detection of localized changes in feature orientation, and provides enough information within the sub-image to ensure adequate signal strength for a robust determination of feature orientation (see below).

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 (Castleman, 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 (Bamber 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 (Scambos 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 (Bindschadler 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 Young 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 Young 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 Young 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 Young 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 Young 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, Young 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 Young and Hyland (2002) velocity vectors. The greater deviation from local flow direction of the flow-direction vectors of Young 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 Young 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 (Jezek, 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 Jezek (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 Young 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 Jezek (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 (Fahnestock 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 (Jezek, 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 (Wu and Jezek, 2004), especially for Lagrangian-type flux-balance codes (e.g. Testut and others, 2003; Roberts 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 (Rignot 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 (Young and Hyland, 2002) and with traced streak lines (Jezek, 1999) supports the view that the Amery Ice Shelf has had a steady flow regime over several decades (King 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.