Pointcatcher software: analysis of glacial time-lapse photography and integration with multitemporal digital elevation models

Abstract

Terrestrial time-lapse photography offers insight into glacial processes through high spatial and temporal resolution imagery. However, oblique camera views complicate measurement in geographic coordinates, and lead to reliance on specific imaging geometries or simplifying assumptions for calculating parameters such as ice velocity. We develop a novel approach that integrates time-lapse imagery with multitemporal DEMs to derive full three-dimensional coordinates for natural features tracked throughout a monoscopic image sequence. This enables daily independent measurement of horizontal (ice flow) and vertical (ice melt) velocities. By combining two terrestrial laser scanner surveys with a 73 days sequence from Sólheimajökull, Iceland, variations in horizontal ice velocity of ~10% were identified over timescales of ~25 days. An overall decrease of ~3.0 m surface elevation showed asynchronous rate changes with the horizontal velocity variations, demonstrating a temporal disconnect between the processes of ice surface lowering and mechanisms of glacier movement. Our software, ‘Pointcatcher’, is freely available for user-friendly interactive processing of general time-lapse sequences and includes Monte Carlo error analysis and uncertainty in projection onto DEM surfaces. It is particularly suited for analysis of challenging oblique glacial imagery, and we discuss good features to track, both for correction of camera motion and for deriving ice velocities.

1. INTRODUCTION

Time-lapse imagery can provide valuable glaciological information, e.g. on glacier extents (Motyka and others, 2003), tidal interactions (Maas and others, 2006; Dietrich and others, 2007) and ice velocities (e.g. Flotron, 1973; Harrison and others, 1986, 1992; Evans, 2000; Ahn and Box, 2010; Schubert and others, 2013). Imagery can be acquired from the ground or above, with satellites capable of providing regular overpasses, typical spatial resolutions of order ~10 m, and datasets covering decadal time spans (e.g. Rignot, 1998; Heid and Kääb, 2012; Shepherd and others, 2012). For detailed analyses (e.g. of calving events: O'Neel and others, 2003; Amundson and others, 2008; Rosenau and others, 2013) repeat imagery may be required on timescales of minutes or hours. Such work is now being enabled due to the step change in spatio-temporal data resolutions resulting from increasing deployment of remote digital time-lapse cameras. However, due to the oblique view from terrestrial vantage points, perspective effects complicate quantitative data processing by varying the effective scale across the images. Furthermore, with single-camera (monoscopic) installations, ice motion towards (or away from) the camera cannot be determined, and horizontal and vertical components can only be uniquely distinguished for specific camera orientations. With vertical surface change from melting forming an important factor in glacier mass balance calculations, a technique to independently extract ice velocity and elevation change from time-lapse imagery captured from general (rather than specific) camera orientations should represent a useful glaciological tool.

Here, we present an approach that enables horizontal and vertical components of glacier surface change to be quantified from a generalised oblique terrestrial time-lapse sequence, by integrating data from multitemporal DEMs. The method is based on deriving three-dimensional (3-D) geographic point coordinates for image feature-tracks within a georeferenced time-lapse sequence, by deriving individual viewing distances for each feature in each image. The view distances are constrained using two DEMs acquired at different times, and by assuming that the planimetric path of each 3-D point is linear over the duration of the sequence (i.e. if viewed from directly above, points would appear to travel in straight lines).

We demonstrate the process using data from Sólheimajökull, Iceland (Fig. 1), where time-lapse imagery was being acquired to assess the potential influence of the Katla volcano on glacier dynamics. The Sólheimajökull sequence represents a highly challenging dataset to process, encompassing all the difficulties that tend to arise in terrestrial time-lapse images: drift due to an unstable camera, variable weather, illumination and snow cover conditions and an ice surface that evolves rapidly due to melting. These characteristics mean that approaches based on automated matching of image pairs (e.g. Messerli and Grinstead, 2015) will have limited success. To address such challenges, we developed a freely available and user-friendly software (Pointcatcher; http://tinyurl/pointcatcher) that implements georeferencing, image registration and Monte Carlo error analysis procedures in a feature-tracking application previously used for the laboratory and volcanic image sequences (Delcamp and others, 2008; Applegarth and others, 2010; James and Robson, 2014). The versatility of the resulting image processing software makes it applicable to terrestrial time-lapse sequences collected from glaciers around the world.

Fig. 1. (a) The Mýrdalsjökull area of Iceland (location arrowed in inset). The box outlines the snout and proglacial regions of the outflow glacier Sólheimajökull shown in the shaded relief map (b) derived from 2010 airborne lidar data (Staines and others, 2015). In (b), the box indicates the surveyed region shown in Figure 8a, with the black circle giving the location of the TLS and time-lapse camera. Coordinates are given in Icelandic National Grid (m).

2. CURRENT TECHNIQUES AND CHALLENGES FOR ANALYSIS OF GLACIAL TIME-LAPSE IMAGES

Significant recent progress has been made in measuring ice velocities from image sequences. Much of the progress has been driven by the increasing availability and resolution of optical and radar satellite imagery, with some of the automated image matching algorithms developed now also being explored for use with terrestrial sequences (e.g. Vernier and others, 2012; Messerli and Grinstead, 2015). Although the practical challenges associated by processing satellite data and terrestrial sequences acquired with consumer imagery can differ substantially, the same basic procedures underpin both workflows: registering sequential images or image pairs together, identifying and tracking image features that represent the glacial surface and converting results, which are initially in image coordinates (i.e. pixels) to geographic coordinates.

2.1. Image registration

Image registration (or co-registration) defines the relationships between different images that enable measurements made in each image to be described within a single reference image coordinate system. Registration transformations account for any changes in the camera or sensor parameters (such as position or pointing direction) between different images, and are usually determined by identifying and comparing image features in areas of static topography. With most modern ground-based time-lapse data acquired from either temporary or semi-permanent remote installations, camera position can generally be assumed to be constant. Thus, the registration process has to account only for small camera rotations, which can result from thermal effects, wind vibration or settling of the installation. Without effective registration, such rotations add noise (e.g. from wind vibration) or systematic displacements (e.g. from settling) to measured feature positions relative to the reference image. In some cases, where the 3-D geographic coordinates of control points are known, image registration can be combined with georeferencing, with each image registered directly to the geographic coordinate system, rather than initially to a reference image.

The quality of image registration (and the error magnitude in subsequent analyses) is a function of the precision with which static points in the landscape can be located through the image sequence and their distribution across images. Static points usually represent easier features to track than those on glacier surfaces (because they do not evolve through time) and are often amenable to accurate and fully automatic tracking. Nevertheless, registration problems arise when features change their appearance (e.g. due to varying snow cover) or are obscured entirely by cloud during periods of poor weather. Understanding the uncertainties involved in image registration is critical to overall error analyses, and a Monte Carlo approach can be used to indicate the sensitivity of the registration to the characteristics of static point measurements in any image.

2.2. Glacier feature tracking

Image-based measurement of glacier surface motion is now commonly carried out using a variety of automated algorithms that match image texture between image pairs (Scambos and others, 1992; Kääb and Vollmer, 2000; Leprince and others, 2007). These algorithms can deliver large numbers of points with matching accuracies of potentially 0.02 pixels under idealised conditions (Maas and others, 2010). However, systematic changes in the appearance of natural features (e.g. due to variation in illumination or melting) can lead to bias in the results and, as features evolve further, matching will eventually fail. Thus, for rapidly changing surfaces, fully automated matching may only be possible over relatively short durations. To follow characteristic features over longer periods, manual individual feature identification and interactive tracking can be carried out (e.g. Eiken and Sund, 2012), although it is relatively time-consuming (restricting the number of features processed) and is unlikely to deliver sub-pixel accuracies. This is the case with the Sólheimajökull data, where rapid image texture changes due to surface melting will preclude effective use of fully automated procedures for glacier surface tracking over any substantial duration. Decisions on the use of interactive or automated tracking can be guided by considering the relative magnitude of feature displacements due to ice motion and the expected measurement and registration errors.

2.3. Georeferencing and 3-D coordinates

To derive 3-D geographic point coordinates from images requires a camera model (which provides a generalised description of how the camera represents any external 3-D scene in its 2-D image), the viewing distance to each observed point, image georeferencing parameters that describe the camera position and how it is oriented within the geographic coordinate system. The application of photogrammetric techniques and stereo time-lapse installations (Eiken and Sund, 2012; Whitehead and others, 2013; James and Robson, 2014; James and others, 2014) is currently being explored to enable viewing distances to be calculated directly from multiple simultaneous images, but accurate results are difficult to attain under practical field conditions. For single-camera systems, georeferencing can be achieved by aligning the image to a contemporaneous DEM through specific control points and, by intersecting virtual rays representing the image observations with the DEM surface (e.g. Messerli and Grinstead, 2015), deriving viewing distances (hence 3-D geographic point coordinates).

In order to process sequential feature observations to derive velocities within an image sequence, an updated DEM should be used for each image, unless the positional change of the surface can be assumed to be negligible or a very specific viewing geometry, which is normal to the directions of surface change and ice movement, is employed. However, commonly, only a single DEM is available, with the implication that any surface changes in the direction of the camera cannot be distinguished. Thus, only in camera views that are perpendicular to ice motion and surface change can horizontal and vertical components of motion be independently determined (e.g. Maas and others, 2006). For more general scenarios where a component of ice motion or surface melt is towards (or away from) the camera, updated view distances are a requirement to separate horizontal and vertical motions.

Our contribution here is to address these challenges by (1) providing software that complements existing image-pair matching methods by enabling a flexible, interactive individual feature-tracking approach for the analysis of difficult time-lapse image sequences, including Monte Carlo error analysis and (2) developing a technique to extract independent horizontal and vertical velocity components from image feature tracks by integrating results with a DEM pair. (3) Finally, we demonstrate use of these in a case study at Sólheimajökull to identify variations in ice velocity and melt rate.

3. METHODS: DATA COLLECTION AND PROCESSING

To derive a high frequency record of ice movement and melting we developed a workflow to combine time-lapse imagery with multitemporal DEMs (Fig. 2). Our time-lapse image processing software, Pointcatcher v2.0, allows automatic, semi-automatic and fully interactive feature-tracking for image registration and for motion detection. Automated tracking is based on normalised cross-correlation of image patches, with selectable patch and search area sizes. The search area for a feature in a subsequent image can be centred on the initial feature coordinates, or can follow feature trajectories, with changes in camera orientation taken into account. For semi-automatic tracking, results computed (using correlation techniques) can be interactively updated by the user along the sequence. This approach exploits the advantages of visual recognition while retaining some speed and accuracy advantages from automated processing and, with difficult imagery, delivers much more sustained feature tracks than could be achieved with a fully automated approach.

Fig. 2. Workflow outline for data acquisition and processing, with greyed boxes indicating processes carried out within the Pointcatcher software.

Our methodology is based on individual feature tracking over prolonged periods, enabling measurement of cumulative feature displacements and calculation of mean velocities and within-sequence velocity variations. Image registration is carried out using observations of static features to derive corrective camera rotations. The image sequence is finally georeferenced by matching with a DEM, and 3-D point coordinates throughout the sequence are determined by integrating a second DEM. Error is assessed by determining the precision of the camera orientation based on the static points, and using a Monte Carlo approach to reproject that uncertainty onto the DEM. For the Sólheimajökull case study, DEMs were derived from terrestrial laser scanner (TLS) surveys carried out at the start and the end of the image acquisition period.

3.1. Study site

Sólheimajökull is an 8 km long, non-surging temperate glacier which drains from the Mýrdalsjökull Ice Cap (Fig. 1a). It supports a maximum ice thickness of 433 m (Mackintosh and others, 2002; Kruger and others, 2010; Russell and others, 2010; Sigurdsson, 2010) and a total area of ~78 km² which extends into the active volcanic caldera of Katla – historically the most destructive subglacial volcanic system, and responsible for routing jökulhlaups down the Sólheimajökull outlet glacier (Lawler and others, 1996; Roberts and others, 2000; Le Heron and Etienne, 2005). The subglacial drainage system seasonally establishes connectivity with the Katla geothermal zone and transports dissolved volatile gases from beneath the glacier in the summer meltwater streams (Wynn and others, 2015). The dynamic advance and retreat cycles exhibited by the glacier (Schomacker and others, 2012) are frequently asynchronous with other glaciers along the south coast of Iceland, and are likely explained through the migration of ice divides (Dugmore and Sugden, 1991). The recent establishment of an ice-contact proglacial lake may also now be having an influence on glacier motion (e.g. Carrivick and Tweed, 2013). However, high temporal resolution assessment of Sólheimajökull's motion has never been captured before.

3.2. Time-lapse image acquisition

Time-lapse images were collected for 73 days between 30 April and 11 July 2013 (days of year: 119–191) using a light weight setup originally designed for monitoring active lava flows (James and others, 2012) where speed of installation and portability of equipment are important factors, and any semi-permanent infrastructure is impractical. The images were acquired from a dSLR Canon 550 camera with a 28 mm fixed focal length lens, and triggered by an intervalometer. The camera was protected in a small, weather-proof box positioned on the ground, aligned appropriately and secured by partially burying with rocks (Fig. 3, inset). Power was supplied from an adjacent 12 V battery, recharged by a 500 mA solar panel (also secured by rocks). Image acquisition was set for an hourly interval as a deliberate oversampling to increase the probability of obtaining good quality images during periods of variable weather conditions.

Fig. 3. Panoramic view (30 April 2013) looking approximately north-east over Sólheimajökull, showing the TLS during data collection. Inset shows the time-lapse camera (left arrow) and solar panel (right arrow), when viewed from the direction of the glacier.

3.3. TLS data acquisition and processing

TLS data were acquired on 30 April and 11 July 2013 using a very-long-range laser scanner (Riegl LPM-321) previously shown to be capable of providing useful data over measurement distances of 3.5–4 km (Schwalbe and others, 2008; James and others, 2009). The scanner was used from elevated ground ~33 m from the time-lapse camera location (Fig. 3), with the camera included in the scans to enable its position to be determined to ~10 cm within the scanner’s coordinate system.

In both TLS surveys, data of Sólheimajökull and the stationary cliffs on either side of the glacier were acquired. The TLS was levelled and oriented such that its coordinate system was approximately aligned with respect to north. Although not critical to the ensuing analyses, the georeferencing of the initial survey (30 April) was then refined by registering the stable cliff area to a pre-existing airborne lidar dataset from 2010 (Staines and others, 2015). The second survey was then similarly adjusted to minimise the differences in the stable cliff regions between the two TLS surveys (using an iterative approach implemented in the Riegl processing software, RiProfile version 1.5.0).

3.4. Image registration

To account for small camera rotations during the sequence, image registration was carried out by tracking image features representing stable topography. Foreground elements of the scene could not be assumed to be stable due to ground heave and disruption by tourists, so reference points were selected covering the surrounding ridges and mountain top areas. Unfortunately, varying illumination and snow cover conditions over most of this ground made for highly variable image texture that was unsuitable to track. Thus, the most reliable features represented irregular sections of horizon, where dark topography was highlighted against light sky or cloud backgrounds. Although this restricted features to a limited region, the horizon crossed the width of the image and provided a feature distribution that was reasonable for constraining camera rotation angles. A number of additional static features were identified to use as ‘check’ points. These were not used in the calculation of the registration parameters, but can be used to give an independent assessment of registration accuracy.

Image registration comprised a two-stage process; a reference image was selected in which the greatest number of features was observed, and initial registrations were derived for all other images, using a fast, robust approach capable of identifying and rejecting outlier feature observations. This procedure is based on an image-based transform which does not account for lens distortion; so during the second registration stage remaining (inlier) features are used to define a registration for each image in terms of a physical camera model – i.e. rotations around the camera's horizontal (omega), vertical (phi) and optic axes (kappa). The camera model includes lens distortion parameters and, based on previous calibration work with the same lens (James and Robson, 2014); only one radial parameter was used for the Sólheimajökull sequence. For images in which the static features could not be observed (e.g. due to complete occlusion of the horizon by cloud), successfully derived camera orientations from a sequential image were used.

3.5. Glacier feature tracking

The glacier surface was monitored by tracking an additional ~50 natural features that were identified as recognisable throughout the sequence (even if their appearance changed significantly). Normalised cross-correlation could be used to facilitate tracking over timescales of days; however, over longer timescales, the evolving surface usually required updating of the correlation template and interactive adjustment and visual assessment of the feature positions was required.

3.6. Georeferencing and data integration

The image sequence was georeferenced to the TLS coordinate system by defining the camera position and orientation. The camera position coordinates were given by the time-lapse camera position identified in the TLS data. Camera orientation was estimated by projecting the TLS data onto the image and adjusting the camera angle until the TLS data appeared best aligned with the image scene. Due to the computer-based matching between image features and topographic data being extremely challenging, the alignment process was carried out manually in Pointcatcher, so it is not associated with formal error estimates. Nevertheless, with the camera used (with a 28 mm lens and 5 µm pixel size) an estimated miss-registration of up to 2 pixels would represent <0.005° of misalignment.

Following georeferencing, the image-based measurements can be transformed into 3-D geographic coordinates through consideration of both DEMs (Fig. 4). For the images taken simultaneously with the DEM acquisitions (at the start and end of the sequence), 3-D coordinates for observed features can be derived by reprojecting the feature locations onto the DEM surfaces. Reprojection is implemented in Pointcatcher by forming a triangular network from the DEM points then, for each feature, using a graphical approach to determine which triangle the image feature ray intersects, and calculating the point coordinates of the intersection. In order to obtain 3-D point coordinates for all other images, each point was then assumed to move within the vertical plane that contained the start and end positions, i.e. in a straight line if viewed from directly above. Thus, for each feature, the 3-D start and end points defined a vertical plane, and unique 3-D point coordinates could be calculated for all other images in the sequence by intersecting the reprojected rays of the image feature with this plane (Fig. 4e).

Fig. 4. Derivation of 3-D point coordinates for a feature observed within an image sequence. (a) Image registration throughout the sequence allows the set of feature observations to be represented in one reference camera orientation, C. (b) For the observation made closest in time to the first TLS survey, 3-D coordinates can be calculated by projecting the observation through the perspective centre of the camera, p, onto the DEM surface defined by the TLS data, DEM 1. (c) The same procedure is carried out with the last point observation and the second DEM, DEM 2. (d) The two 3-D points are then used to define a vertical plane in which the point is assumed to lie at all other times. (e) 3-D coordinates for all other image observations of that point are then calculated by intersecting their observation rays with the plane.

With interactive (rather than automated) feature tracking, noise levels may be undesirably high in individual tracks, and velocities may be determined better by averaging several points. Thus, for each point, cumulative displacement values were normalised to the total path length travelled, and then modelled with a best-fit straight line (i.e. a constant velocity model). For each image, the differences between the modelled and measured normalised cumulative displacement values were averaged for all points, to give a mean variation from the constant-velocity position, described as a fraction of full path length.

3.7. Error estimates

With the underlying objective being to assess variations in melt and flow processes, absolute georeferencing of the data was not critical for the analyses, so we focus on relative error along the sequence. Consideration of error within camera models or DEMs is also outside the scope of this work, although could be added to the analysis in the future.

The relative image registration quality represents how closely the position of a theoretical, perfectly measured static feature can be reproduced in different images. When using a camera model with a fixed location, registration quality can be characterised by the uncertainty in estimated camera rotation angles, which is a function of the number and distribution of tracked static features, and any errors in their coordinates. Thus, to assess the quality of image registration, a Monte Carlo approach can be used in which orientation values are estimated repeatedly, with different randomised errors (i.e. perturbations) added to the static points’ image position for each estimation. The perturbations are taken from a pseudo-random normal distribution with a standard deviation (SD) that reflects the precision of the image measurements. This precision can be estimated directly by considering the distances (in pixels) between the static features in a registered image and the equivalent features in the reference image; the RMS of these residual distances indicates the precision of the static feature measurements made during tracking.

Uncertainty in the image registration is thus defined by distributions of likely camera orientation angles. The overall precision of feature measurement for any particular feature within a registered image sequence, σ _mr, is then a combination of the effect of image registration error within the region of the feature, σ _r, and the image measurement precision, σ _m, representing how well the coordinates of that feature can be located in different images:

(1) $${\sigma _{{\rm mr}}} = {\left( {\sigma _{\rm r}^2 + \sigma _{\rm m}^2} \right)^{1/2}}.$$

Finally, the relative error for a point in 3-D can be derived by reprojecting this uncertainty onto a DEM surface. The reprojection process means that error magnitude in geographic coordinates increases for both increasing viewing distance and decreasing angle of incidence on the surface, with the resulting planimetric uncertainty for any one point unlikely to present a normal distribution. In Pointcatcher, the Monte Carlo implementation represents the uncertainty as a distribution of likely point positions on the DEM. Furthermore, the reprojected uncertainties from the two images can be combined and their influence on the velocity estimates assessed.

4. RESULTS

Over the 73 days between TLS surveys, 1768 time-lapse images were acquired. Due to Sólheimajökull surface velocities being ≤0.2 m d⁻¹, the sequence was down-sampled to 145 images to represent the best available at a ~12 h interval.

4.1. Image registration

Due to the relatively constant appearance of the static features (Fig. 5a), their locations could be dominantly tracked throughout the sequence using a fully automatic approach. Manual intervention was required when there was a sudden, large change in camera orientation (resulting from camera movement during maintenance and data retrieval on Days of Year 119 and 189). In this case, using a few manually tracked features to make an initial estimate of a new camera orientation then enabled successful automated tracking of the remaining features.

Fig. 5. Example features tracked, as shown by patches of 31 × 31 pixels extracted at five different times and centred on the features’ locations. Feature number is given in top left of first image for each feature, and relates to the labels in Figure 6a. The consistency of the horizon-based static features (a) used for registration contrasts with the significant evolution of the glacier surface features (b).

Fig. 6. (a) Image (30 April 2013) showing the position of tracked features adjusted for camera rotation; static points are shown by crosses located on the horizon, red for those used to derive camera orientation and white for those used as check points. (b) The number of static points visible and determined as inliers for orientation calculations, and the number of visible check points in each image. (c) The relative camera rotation angles derived for each image, with the sharp step in Phi due to camera disturbance during data retrieval. The shaded bands represent the uncertainty in the angle estimations, magnified by a factor of 10 for visibility. The largest peaks in uncertainty (particularly in kappa, rotation around the optic axis) are due to cloud obscuring static points on one side of image. (d) The quality of the image registrations are indicated by the RMS residuals on the transformed static (orientation) and check points, and are dominantly <1 pixel.

Typically, >20 static features were successfully observed in an image (Figs 6a and b), resulting in camera orientations being directly calculated for 132 images of the sequence. In the remaining 13 images, low cloud completely obscured the distant and elevated topography used for the static points and orientation information had to be propagated from preceding or following epochs. As well as the manual interventions, the calculated camera orientation angles (Fig. 6c) show initial variability early in the sequence due to ground heave, then a gradually declining rotation, presumably related to propagation of a thawing front deeper into the ground.

For each image in which the static points could be observed, the Monte Carlo analysis (2000 simulations per image) resulted in camera angle distributions that all passed Chi-squared tests for normality, indicating that a SD statistic would reasonably reflect the relative orientation precision. Thus, image registration delivered three camera angles per image, along with associated precision estimates that varied depending on the number, distribution and quality of the static point observations in each image (Fig. 6c).

The RMS of the residual magnitudes on the static points used for image registration show a mean of 0.51 pixels from along the sequence (Fig. 6d). The RMS residual for the check points is similar, but rises notably after Day of Year 183 to ~1.2 pixels. This, along with the strongly decreasing number of observed points due to cloud cover (Fig. 6b), indicates increased uncertainty in image registration for the last ~10 days of the sequence. Note that after Day of Year 190, cloud prevented observation of the static points in all but one half of one image thus, for this period estimates of camera orientation are weak, with evaluation of precision not possible for most images.

4.2. Glacier feature tracking

On the glacier surface, varying snow cover and rapid melting resulted in image texture that evolved too quickly to be tracked over sustained periods along the sequence. Consequently, most of the glacier features tracked represented the base of dirt cones (Fig. 5b) because, although these also evolved significantly, with some even disappearing to leave only faint surface traces, their junctions with the glacier surface generally provided strong image contrast. This evolution meant that semi-automatic tracking was required with frequent updates of the correlation template to maximise the length of tracks. The ‘noisiest’ tracks were those from the most rapidly changing cones, which needed operator interaction the most during tracking and in which identification of the same representative point in subsequent images, even visually, could be challenging.

For the Sólheimajökull data, we estimate that the manual measurement precision, σ _m, of the evolving natural features can be approximated as ~1 pixel. Using Eqn (1) to give a sequence-wide precision estimate by assuming that the mean RMS error on the static points (0.51 pixels) represents an indicative estimate for registration accuracy across the images and along the sequence, gives σ _mr ≈ 1.1 pixels. Inspection of the feature tracks supports this value by suggesting a general noise magnitude of ~1–2 pixels (Fig. 7b). In limited cases (e.g. the start of track 2; Fig. 5b), error magnitudes of up to ~5 pixels are shown, indicative of difficulties in feature identification.

Fig. 7. Image feature tracking for three points, using either correlation only (grey) or manually assisted correlation (black) tracking. (a) Track continuity through time is shown with the bars representing the periods in which the features were identified, and the arrows indicating when the reference template used in the automatic correlation-only tracking was either set or reset. For correlation tracking, a threshold of 0.6 was used to determine successful matches to the template. (b) Changes in feature image pixel coordinates are shown (after correction for camera orientation changes, and with tracks moved adjacent to one another for clarity).

4.3. 3-D geographic coordinates

The differences between the DEMs show a reasonably consistent −3.0 m vertical change across the scene, apart from where the Katla debris band lies (Fig. 8). Transforming the feature tracks into 3-D point coordinates indicated that they represent horizontal paths of 6–15 m in length, denoting mean horizontal velocities of 0.06–0.20 m d⁻¹, with the slower velocities located towards the edge of the glacier (Fig. 8b). Velocity uncertainties resulting from image measurement and relative registration error (Fig. 8) are calculated by determining all combinations of path length between Monte Carlo simulations for the first and last images (500 perturbed positions for each point in each image, giving 2.5 × 10⁵ simulated path lengths per point). The uncertainties are due to those of the reprojected point positions, so are independent of the glacier movement and reflect the relative orientations of the camera view and the surface onto which it is reprojected. Thus, uncertainty ellipses are aligned towards the camera, and increase in magnitude with viewing distance and with decreasing incidence to the surface.

Fig. 8. (a) Surface change in the boxed region of Figure 1b, between 30 April and 11 July 2013. Elevation change, derived by differencing 2 m resolution DEMs (generated from the TLS data using Surfer (version 9.11) software), is given by the shading. The position of points analysed in the time-lapse sequence are given by black dots, with the associated vectors showing their total horizontal displacement over the period (note the vector scale). The ellipses illustrate the uncertainty in the displacements and are magnified by a factor of 10 for visibility. Inset shows the full distribution of Monte Carlo displacements calculated for Point 2 (labelled). Plotting horizontal velocity magnitudes against Northing (b) illustrates the increase in velocity towards the glacier centreline. Velocities calculated directly from the horizontal displacements between the start and end point positions (i.e. as in (a)) are given in grey, with associated error bars representing ± 1 SD, and the vertical line giving the mean value of the points it overlaps. For each point, the black symbol represents the velocity value obtained from straight line fits to all its displacement data.

Towards the centre of the glacier, the measurements suggest a region of constant horizontal velocity, for which the 19 points furthest from the glacier margin (Fig. 8b) give an overall mean of 0.170 m d⁻¹, with a SD of 0.022 m d⁻¹. However, variability between these measurements is not fully explained by the magnitude of their error bars. Although this could be interpreted to represent fine-grained horizontal variation in surface velocity, it is important to recognise that these velocities are derived from the bases of dirt cones, and cone evolution could add variability that is not captured within the Monte Carlo-based error bars.

Separating point tracks into cumulative horizontal and vertical displacement components (i.e. displacements from their initial point positions) suggests that many points express small but systematic deviations from constant velocity (Fig. 9a). Normalising the cumulative displacements by path length then averaging all paths, highlights these temporal variations (Fig. 9b). This indicates that the variations are systematic between different points, with periods of slower-than-average ice movement illustrated by negative gradients (Fig. 9b), where point positions are gradually falling behind those given by the constant (mean) velocity models. Positive gradients indicate periods when points are ‘catching up’ or overtaking the positions derived from the mean velocity models, thus indicating periods of faster-than-average ice velocity.

Fig. 9. (a) Examples of point displacement measurements along the entire sequence (point numbers, corresponding to labels in Figure 6a, are given in the square brackets). The grey dashed lines show linear fit models to the displacement data, from which mean point velocities can be derived. For clarity, the vertical components have been offset by 1 m to separate the data. (b) Mean deviations in point displacement from the constant velocity models, with the grey bands illustrating the standard error of the mean at each epoch. Negative gradients indicate periods of slower-than-mean velocity and positive gradients indicate periods of faster-than-mean velocity. The dashed lines in the upper panel give linear fits to different periods, and are labelled with the relative change in velocity with respect to the overall mean.

Over the course of the sequence, both horizontal and vertical velocity components demonstrate periods of high and low velocity with, for the horizontal component, magnitudes of up to ~10% velocity variation over durations of ~25 days. The brief period at the end of the monitoring period, where horizontal velocity appears to almost double (Fig. 9b), occurred when the static points were generally obscured by cloud (Fig. 6b). Consequently, during this period, image registration has to be assumed to be weak, with estimates of orientation precision only available for the last image. Thus, we consider the apparent rapid velocity change as likely to be an artefact of image registration error, and it is not interpreted as representative of a glacial process.

5. DISCUSSION

The results demonstrate our semi-automated, interactive individual feature tracking methodology for image-based ice velocity measurement, as a complementary approach to the fully automated image-pair matching commonly used with satellite data (e.g. Heid and Kääb, 2012). Whereas automated techniques work well when image changes are dominated by ice flow, interactive approaches should be considered when image variations are complex (e.g. from surface melting).

5.1. Good features to track

With automated tracking techniques, the use of natural features for image registration and velocity estimation can give problems due to variations in their image texture through time. In the Sólheimajökull dataset, static features on the surrounding cliffs were too affected by illumination changes and varying snow conditions to enable reliable tracking for image registration. However, features on the horizon proved suitable, even though the saturated or near-constant brightness values in the sky give no useful areas of image texture for cross-correlation. Thus, the correlation signal is reliant on the sharp contrast at the horizon and, to give good localisation in both x- and y-directions; in any regions used the horizon should not be highly linear (Fig. 5a).

Much of the recent work on glacier feature tracking has relied on heavily crevassed surfaces to provide the image texture required. The surface of Sólheimajökull is much more varied, with image texture also presented by features such as volcanic tephra layers, cones and thrust planes. Care must be taken to select appropriate features to track, i.e. features that are not only persistent through the image sequence, but whose movement is also representative (or as representative as possible) of surface displacement. Features such as those representing thrust planes must be avoided when surface melt rate is not negligible with respect to ice movement, because calculated velocities will otherwise reflect a combination of glacial motion and melt-back along the reverse inclined plane.

5.2. Interactive individual point tracking or automated image-pair matching?

For guiding decisions on which approach to use for a given image sequence, the durations over which successful tracking can be carried out can be compared with the anticipated duration required for a given signal-to-noise ratio (SNR, i.e. a ratio of pixel displacement representing ice motion, to apparent pixel displacement due to other factors). To estimate the number of image intervals, i, to achieve a specificSNR we can consider the mean expected image displacement per image interval, d, and the overall image measurement error (Eqn (1)), which comprises both error in the image feature measurement, σ _m, and the image registration, σ _r:

(2) $$i = \displaystyle{{SNR \times \sqrt {{\sigma _{\rm m}}^2 + {\sigma _{\rm r}}^2}} \over d}.$$

For the Sólheimajökull sequence, d ≃ 1 pixel (~150-pixel-long displacement tracks over the 144-image sequence; Figs 6 and 7) and, as previously discussed, σ _m ≃ 1 pixel and σ _r ≃ 0.5 pixel. Thus, for a SNR ratio to exceed 10, measurements should be taken over intervals of >~11 images. Successfully automated image matching might reduce σ _m to typical values of ~0.1 pixel (or possibly smaller), but the overall error term will remain high, constrained by the image registration component, σ _r. In this case, 5 image intervals would still be required for SNR > 10. Note that although σ _m can be effectively reduced by averaging multiple features, this will not similarly reduce σ _r, because registration error is systematic across any one image. Thus, averaging over ~30 features also results in i ≈ 5. Such analyses can help determine processing strategy, but equally, for longer sequences, gives an indication of the expected duration over which systematic change (i.e. ice motion) can be reliably detected. Thus, at Sólheimajökull, we expect to require a duration of ~2.5 days to get reasonable ice velocity measurements, and therefore, ~25 days to determine variations of 10% with the same confidence (as reflected in Fig. 9b).

5.3. Reprojected uncertainties

Although such a broad error analysis is useful for consideration of measurement strategy, it does not specifically consider error variability within sequences and across images. Such variability is captured by the Monte Carlo approach, where camera orientation uncertainty is determined for each image, and directly applied to each specific feature measurement. Examination of the uncertainty in camera orientation angles shows that it is generally well represented by Gaussian probability distributions, and can thus be appropriately defined by statistics such as SD. However, only under restricted circumstances and with a flat DEM surface, will the use of these distributions to reproject image features onto a DEM result in similarly near-normally distributed geographic coordinates. Consequently, in most scenarios, the uncertainties associated with planimetric positions are not normally distributed so, although there are limited practical alternatives, the association of point co-ordinates and velocities with SDs should be treated with caution.

In our analyses, we focus on down-sequence variability, and thus neglect the constant errors in overall georeferencing, the DEM and camera model. For a rigorous treatment of absolute error, the use of static ground control points is recommended as one of the best ways to incorporate overall georeferencing precision into the error propagation. In many cases, repeat, high resolution TLS data may not be available, and only a single low-resolution DEM (e.g. tens of metres) may be available. In this case, Pointcatcher can still be used for tracking, registration, georeferencing and reprojection, but separation of horizontal and vertical velocity components will only be possible for specific camera orientations. The errors associated by using a low resolution DEM could be assessed by applying offsets to the DEM and considering the variation in reprojected point positions. For relative measurements over long viewing distances at favourable angles, errors may be deemed acceptable.

5.4. Ice velocity variation at Sólheimajökull

The synchronous velocity variations detected for the different points measured on Sólheimajökull indicate systematic velocity changes during the image sequence (Fig. 9). Horizontal velocities varied by ~5–10%, over timescales of ~25 days, and had changes that were asynchronous with those in the vertical component, suggesting process independence. Such independence may be due to time delays incurred from surface meltwater transit to the glacier bed where water pressurisation can promote localised reduction in basal drag, or a longitudinal coupling of ice dynamics in which local movement is driven by upstream conditions independent of local surface melt dynamics. A detailed understanding, to include any basal melting contribution from Katla will require additional measurements, which we leave for future work.

6. CONCLUSIONS

We present new software to facilitate quantitative measurement from oblique time-lapse imagery of glaciers. As well as providing a straightforward feature tracking application, the software enables integration of time-lapse data with DEMs, and includes error analysis based on the projection of Monte Carlo-based uncertainties onto DEM surfaces. Using this, we implement a novel approach for deriving independent horizontal and vertical ice velocity components using only two DEMs acquired at different times within the image sequence. Illustrating the process on a 145 image sequence from Sólheimajökull, indicates a mean ice velocity of 0.170 m d⁻¹ (19 measurements, SD, 0.022 m d⁻¹) at distances >~200 m from the glacier edge during May–July 2013. Normalising cumulative point displacements by their overall path lengths enables averaging to be used to mitigate measurement noise and to reveal systematic variations in horizontal ice velocity of ~5–10%, which were asynchronous with vertical velocity changes.

Pointcatcher, our software for carrying out time-lapse feature tracking, down-sequence image registration, Monte Carlo error analysis, georeferencing and DEM integration is freely available over the web: http://tinyurl/pointcatcher. Future work will build on this framework to include direct registration of cameras using control point resection, and implement image-only 3-D measurements through stereo image sequences (James and Robson, 2014).