A semi-automated method for mapping glacial geomorphology tested at Breiðamerkurjökull, Iceland

Remote Sensing of Environment 163 (2015) 80–90

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

A semi-automated method for mapping glacial geomorphology tested at Breiðamerkurjökull, Iceland Ciaran Robb a,⁎, Ian Willis a, Neil Arnold a, Snævarr Guðmundsson b a b

Scott Polar Research Institute, Lensfield Road, University of Cambridge, Cambridge CB2 1ER, UK Institute of Earth Sciences, University of Iceland, Sturlugata 7, Askja, 101 Reykjavik, Iceland

a r t i c l e

i n f o

Article history: Received 7 October 2014 Received in revised form 11 March 2015 Accepted 12 March 2015 Available online 8 April 2015 Keywords: Geomorphological mapping Lidar OBIA Land surface parameters Glacial

a b s t r a c t Despite advances in data collection and analytical software, mapping of glacial geomorphology is still largely achieved through manual techniques. Manual mapping is a relatively slow, inefficient and potentially subjective process, particularly given the extensive spatial coverage of remotely sensed datasets. Thus, glacial geomorphic mapping would benefit from semi-automated techniques able to map large areas quickly and efficiently. Here, we develop novel semi-automated techniques for mapping glacial geomorphology, based on lidar derived data collected from Breiðamerkurjökull, Iceland using Object Based Image Analysis (OBIA). The landforms classified were moraines, eskers, flutes/drumlins, and glacio-fluvial outwash. Several land surface parameters (LSPs) were calculated during image segmentation and analysis, representing the characteristics of the landforms of interest. The segmentation and classification routines were developed in Matlab, whereas most attempts so far have used specialised OBIA software. Our techniques utilised exploratory data analysis, where statistical groupings within the data define different landform types. Accuracy was assessed against a reference map derived using traditional manual mapping techniques. Landforms were classified to an overall accuracy of 77%. Our semi-automated mapping techniques should, with further algorithm development, be applicable using other datasets commonly used in palaeo-glaciology such as radar derived elevation data and swath bathymetry. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Geomorphological mapping forms a key component of palaeoglaciology. Developing new and improved techniques for assessing the characteristics of previous glacial environments is helpful for improving models of glacier response to future climate change. Remotely sensed data has provided an excellent resource from which the extent and dynamics of past glaciers and ice sheets have been inferred from their geomorphological signatures (e.g., Bradwell et al., 2008, Clark, Hughes, Greenwood, Jordan, & Sejrup, 2012; Clark et al., 2004; Evans, Clark, & Mitchell, 2005; Ottesen et al., 2008; Salcher, Hinsch, & Wagreich, 2010). The data tend to be used in two forms; passively sensed multispectral data and actively sensed elevation data, the latter of which is often in the form of a digital elevation model (DEM). With advances in technology, improved data acquisition and greater data availability, glacier geomorphology has been mapped with increasing detail (Napieralski, Harbor, & Li, 2007). Light detection and ranging (lidar) is emblematic of these advances and has the potential to represent the land surface with sub-metre spatial resolution. With lidar data, there is the potential to identify and map smaller scale landforms than those

⁎ Corresponding author. E-mail address: [email protected] (C. Robb).

http://dx.doi.org/10.1016/j.rse.2015.03.007 0034-4257/© 2015 Elsevier Inc. All rights reserved.

that have been mapped hitherto, such as flutes and annual push moraines. Image processing and Geographic Information Systems (GIS) software such as Matlab, Erdas, ArcGIS, QGIS and Whitebox provide the platforms on which to compile and analyse elevation data. The advent of GIS has facilitated the compilation of maps and databases of glacial landforms such as the BRITICE (Map and GIS database of glacial landforms and features related to the last British Ice Sheet), which has since been utilised in a number of studies (Clark et al., 2004, Evans et al., 2005). Despite the advances in data acquisition and analytical software, traditional manual techniques remain the most widely used approach for mapping glacial geomorphology from remotely sensed data. However, manual techniques are time consuming and potentially subjective, particularly with larger datasets and more ambiguous landforms. A semi-automated method of terrain interpretation and landform detection may serve to streamline the process of geomorphological mapping in palaeo-glaciology, opening up larger and less accessible areas to research into past glacier and ice sheet extent and dynamics. Such methods may serve to improve consistency in geomorphological mapping by providing precise metrics to differentiate between landforms, particularly those of a more ambiguous nature. Whilst not applied in palaeo-glaciology specifically, many attempts have been made to automate landform mapping in other disciplines, including geology, geomorphology, soil science, hydrology, ecology and

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

conservation (Bishop, James, Shroder, & Walsh, 2012). These use techniques developed in image processing and remote sensing and can be divided into: • Pixel-based methods (e.g., Iwahashi & Pike, 2007; MacMillan, Martin, Earle, & McNabb, 2003; Reuter, Wendroth, & Kersebaum, 2006; Tagil & Jennes, 2008). • Object Based Image Analysis (OBIA) (e.g., Anders, Seijmonsbergen, & Bouten, 2011; Drăguţ & Blaschke, 2006; Saha, Wells, & Munro-Stasiuk, 2011; Schneevoigt, van der Linden, Thamm, & Schrott, 2008; van Asselen & Seijmonsbergen, 2006).

Pixel-based methods are based purely on the individual cell values within an image whereas OBIA involves analysis of the properties of objects within the image. For the purposes of this paper, we use the term “image” to refer to any raster-based dataset, including a DEM. We elaborate on these two techniques in turn below. Pixel-based techniques begin with the generation of land surface parameters (LSPs) using moving window calculations, which pass a window of specified size and shape systematically through an image, performing a statistical operation within the window (Hengl & Reuter, 2009; Lillesand, Kiefer, & Chipman, 2008). LSPs characterise terrain based on measurements such as slope angle, aspect, curvature, flow accumulation and relative elevation, as the variability in elevation alone is usually insufficient for quantitative geomorphological mapping (Hengl & Reuter, 2009). Specifying LSP threshold values to isolate specific areas within an image is a commonly used approach within this method (Hengl & Reuter, 2009). Thresholding of LSPs has been the mainstay of many studies aiming to characterise the land surface into landform elements, the constituent parts of landforms (e.g., Pennock, Zebarth, & De Jong, 1987; Reuter et al., 2006). Landform elements include saddles, shoulders, backslopes and footslopes, based on surface discontinuities and geometric attributes (Hengl & Reuter, 2009). Whilst landform elements are effective for characterising basic land surface characteristics, combining landform elements to map discrete landforms is more problematic. In the context of glacial geomorphology, many landforms exhibit similar characteristics and are likely to be represented by similar sets of landform elements, making their separation into definitive landforms difficult. Given data of increasingly fine spatial resolution (e.g., lidar), a landform element based method would likely oversegment a DEM, requiring merging steps to produce a representative segmentation for landform mapping. The pixel-based techniques discussed so far use an expert's judgement to guide the classification process, whereas there are other pixel-based techniques that are not dependent on user input. Iwahashi and Pike (2007) use a nested means algorithm to map landforms based on slope, surface texture and local convexity. The number of output classes can be altered to accommodate low relief landscapes by increasing the number of classes (Iwahashi & Pike, 2007). This method represents a computationally efficient way of dividing up a landscape using only three parameters. Stepinski and Bagaria (2009) build on this work to map terrain on Mars. They add a secondary process whereby the initial classification is re-clumped based on the frequency and texture of landforms, producing a terrain map of more clearly defined boundaries with homogenous interiors (Stepinski & Bagaria, 2009). Both these techniques are applied at a coarse spatial scale (30 m +) and identify large scale landforms types (e.g., mountain ranges, plains). They are also unsupervised, with the resulting classification requiring interpretation after processing. The automated techniques described so far analyse individual pixel values to classify an image and are best suited to relatively homogenous, coarse spatial resolution data (Iwahashi & Pike, 2007). Given that the overall morphometry of a landform provides one of its most important defining characteristics, the lack of analysis of the resulting classes is a notable shortcoming of pixel-based techniques. OBIA provides an

81

analytical framework akin to human interpretation, where raster segments are analysed for not only their cell-based but also their geometric properties in order to classify remotely sensed imagery (Blaschke, 2010). In this way, OBIA presents a solution to many of the shortcomings associated with pixel-based techniques. OBIA is particularly well suited to fine spatial resolution data as it produces an image consisting of homogenous segments for analysis, whereas pixel-based analysis of fine spatial resolution can appear noisy (Blaschke, 2010). OBIA is a two stage process, beginning with segmentation of input imagery, followed by analysis of the resulting segments using underlying pixel and geometric properties. Either a hierarchical decision tree-based system or machine learning algorithm is used for classification (Blaschke, Lang, & Hay, 2008). Segmentation requires input parameters, such as the degree of homogeneity (the extent to which constituent pixels are similar), to exert some control on the size, shape and number of segments (Blaschke, 2010; Saha et al., 2011). Finding optimal segmentation parameters usually calls for trial and error experimentation, although some attempts have been made to semi-automate this process (for examples, see Anders et al., 2011; Drăguţ, Eisank, & Strasser, 2011). Furthermore, the detail present in fine spatial resolution data limits the effectiveness of homogeneity-based segmentation techniques, often resulting in over- or under-segmentation of the data (Derivaux, Forestier, Wemmert, & Lefèvre, 2010). However, over-segmentation can be resolved through merging segments based on pixel, geometric or proximity characteristics. Additionally, careful choice of LSPs and moving window sizes prior to segmentation may reduce the need for subsequent merging routines. Anders et al. (2011) attempted to optimise segmentation parameters by using the constituent pixel statistics from manually delineated landforms in a training dataset. For a given landform, the parameter set that produces a segmented dataset (such that the constituent pixel characteristics of segments best matches those of the manually delineated landforms in the training dataset) is subsequently used to map the landform type across the entire dataset (Anders et al., 2011). A key limitation of this optimisation approach is that it relies on an element of manual interpretation in the first instance; the training examples of each landform type must be derived from a geomorphological map. Given the variability of scale and shape within various categories of glacial landform, obtaining representative training samples involves an element of subjectivity. Drăguţ et al. (2011) provide an alternative method of segmentation parameter optimisation based on the work of Woodcock and Strahler (1987), where the local variance within successive segmentations is used to determine the best parameter set. Both methods require multiple segmentations of different LSPs using different segmentation criteria, which is a time consuming process. Eisank, Smith, and Hillier (2014) optimise segmentation parameters for delineating synthetic drumlins based on a dataset developed by Hillier and Smith (2012). However, it is questionable whether idealised data such as synthetic DEMs and drumlins are representative of reality, and whether such a segmentation is then applicable to a real dataset. Furthermore, the optimisation techniques discussed are all reliant on specific software (Trimble eCognition) and its specific algorithms which cannot be used on other platforms with different segmentation algorithms. Post-segmentation merging techniques provide a relatively simple alternative to these optimisation approaches, where the segments are merged during the classification process to more suitable representations of the target classes. Importantly, this process is not specific to one software platform or algorithm. Thus, rather than optimising segmentation parameters based on training samples, our study will calibrate the OBIA process based on the statistics present within the data using segment merging techniques, clustering techniques and exploratory data analysis, as well as detailed knowledge of glacial geomorphology. For a semi-automated methodology to be more effective than traditional methods, it must be a more efficient process as well as an accurate one.

82

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

The objectives of this study are threefold: • To develop a robust and efficient semi-automated methodology for mapping the geomorphology of glacier forelands. • To use and assess the potential of exploratory data analysis and clustering techniques on object properties to categorise landforms based on their statistical similarities and differences. • To assess the accuracy of the resulting classification using the confusion matrix method as used in most remote sensing research.

2. Study area Glaciers and ice caps currently cover ~11% of Iceland's total land area (Jóhannesson et al., 2011). Iceland's largest ice cap, Vatnajökull is one of the largest ice masses outside Greenland & Antarctica (Björnsson & Pálsson, 2008). Vatnajökull's forelands are particularly informative in glacial geomorphic research because they contain a varied, well preserved suite of erosional and depositional landforms from both surgetype and non-surge-type glaciers. Thus, they provide an excellent test bed for developing novel techniques for mapping glacial geomorphology. Of these glacier forelands, Breiðamerkurjökull's is amongst the most extensively studied, most notably by Evans and Twigg (2002). Breiðamerkurjökull is a large outlet glacier on the south eastern side of the ice cap characterised by 4 lobes marked by large medial moraines (Fig. 1A). The foreland contains many lakes, the largest of which is Jökullsarlon, into which the eastern lobe calves. Landform assemblages are typified by areas of fluted/drumlinised till, punctuated by ice marginal depositional areas of push moraines and glaciofluvial material.

Smaller push moraines typically demarcate annual limits of the glacier front, whilst larger ones indicate stationary periods of sustained sediment deposition (Evans & Twigg, 2002). The foreland also includes numerous simple and anabranching eskers formed during recessional phases (Evans & Twigg, 2002). Evidence of glacier surging exists on the eastern lobe section of the foreland in the form of crevasse fill ridges, with historical records suggesting numerous surges since the late 18th century (Björnsson & Pálsson, 2008). The combination of extensive previous research and geomorphological mapping, showing evidence of varied and well preserved landform assemblages, make Breiðamerkurjökull's foreland an ideal area for benchmarking a semi-automated geomorphological mapping methodology. It is hoped that methods developed here will be applicable to other outlet forelands around Vatnajökull and beyond. Development is focussed on distinguishing between discrete landform types around the margins of Vatnajökull beginning with the foreland of Breiðamerkurjökull, as existing geomorphic mapping data of this area is amongst the most detailed and accurate available (Evans & Twigg, 2002). 3. Methods 3.1. Data preparation The following glacial landforms will be mapped semi-automatically as originally mapped manually by Evans and Twigg (2002) who used aerial photographs and field observations. • Moraines • Flutes and drumlins

Fig. 1. (A) Study area (light green) with glaciers (white) superimposed on a lidar hill shade. The black arrow represents the camera position for the corresponding photograph in C. The blue box is the subset area for subsequent LSPs shown in Fig. 3. (B) The geographical context of the study area (red inset box) within southern Iceland showing the major ice bodies superimposed on a Shuttle Radar Topography Mission (SRTM) DEM of the area. (C) A typical scene looking north-west across the Breiðamerkurjökull foreland.

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

83

• Eskers • Glacio-fluvial outwash

The scale of the Evans and Twigg (2002) map is 1:30,000. One of our goals is to assess the accuracy of our semi-automated method, using the manually produced map of Evans and Twigg (2002). However, there has been over 600 m of glacier retreat exposing new glacier forefield since Evans and Twigg (2002) produced their map. Furthermore, the Evans and Twigg (2002) map is highly generalised and illustrative in places (for example flutes/drumlins are shown as ridge crest lines rather than as specific geomorphic features) and a quantitative comparison with a semi-automated digital map would therefore not be possible. For these reasons, we have produced a new manually digitised map, guided by the Evans and Twigg (2002) original, but using new airborne lidar data (acquired in 2011) of Breiðamerkurjökull's proglacial area. The lidar data were collected by Topscan Gmbh for the Icelandic Meteorological Institute. Point cloud data spacing was approximately 0.3 points per m2, with a vertical accuracy of 0.5 m. The point cloud data were gridded to give a DEM of 2 m spatial resolution using a modified version of the linear prediction method as implemented in the SCOP.DTM software (Institute of Photogrammetry and Remote Sensing, 2002). Areas not relevant to the analysis, notably the glacier, surrounding hills and the proglacial lakes of Jökullsarlon and Breiðalorn, were masked to produce a DEM of the proglacial area (Fig. 1A). Hillshades of this DEM and a slope raster were used to pick out the main geomorphic features. Guided by the Evans and Twigg map and by 2012/13 fieldwork undertaken in areas of recently exposed forefield by two of the current authors (Robb and Guðmundsson), the DEM was digitised manually to produce an updated version of the Evans and Twigg map that could be compared quantitatively with a map produced using the semi-automated methods. Fig. 2 presents a summary of the workflow associated with the semiautomatic mapping procedure. It is based on the 2 m resolution DEM described above. 3.2. Land surface parameters (LSPs) The glacial landforms in our study area have distinct morphologies which may aid in their identification. Push moraines typically have asymmetrical, often discontinuous profiles, with gentle ice proximal and steep ice distal slopes (Bennett, 2001; Evans & Twigg, 2002). Flutes and drumlins are smooth ridges in profile and elongate in planform. Eskers are steep and symmetrical in profile often with anabranching ridges (Evans & Twigg, 2002). Glacio-fluvial outwash fans are low lying, relatively flat areas, characterised by many sinuous channels. Four LSPs were chosen for this study to represent the various characteristics of glacial landforms: i) elevation percentile to represent relative elevation for the delineation of ridge, plain and pit forms; ii) slope, a proxy indication of amplitude, asymmetry and discontinuity (Hengl & Reuter, 2009), to distinguish the differing slope characteristics of flutes/drumlins and moraines; iii) D8 based stream network (Tarboton, Bras, & Rodriguez-Iturbe, 1991), from which to calculate sinuosity; and iv) sinuosity (Dilts, 2009), to identify areas of fluvial from non-fluvial incision/deposition, based on the presence of many sinuous channels. Elevation percentile was calculated in Whitebox Geospatial Analysis Tools (Lindsay, 2014). The D8 valley stream network was calculated using TAU DEM software (Tarboton, 2001) and sinuosity was derived using an ArcGIS toolbox developed by Dilts (2009). Of these, elevation percentile and sinuosity had a prescribed moving window size. In the case of elevation percentile, this was chosen qualitatively to maximise internal homogeneity of the landforms of interest without blurring their boundaries with other landforms. The sinuosity radius was chosen qualitatively to be large enough to cover densely spaced valley networks without including more sparsely distributed valleys. Sinuosity is calculated by calculating the ratio of the Euclidean distance between

Fig. 2. Summary of the workflow methodology, starting at the green panel (Lidar DEM) and finishing at the red panel (Accuracy assessment).

stream extrema to the actual stream length within the moving window (Dilts, 2009). Slope and stream sinuosity were used for analysis purposes only later in the study and were not used in the segmentation process. The D8 stream network was generated as an intermediate step in calculating sinuosity and was not used in further analysis. Fig. 3 shows subsets of the LSPs used in this study. Having defined the LSPs, the next stage was to choose an LSP to produce a segmentation image appropriate for each landform. 3.3. Segmentation For OBIA to be effective, the segmentation on which it is based must be representative of the objects of interest, in this case various glacial landforms. Thus, the LSP from which the segmentation is derived must represent adequately the variance exhibited by the landforms within the elevation data. Elevation percentile fits this requirement as it represents relative height within the DEM and thus peaks, ridges, pits, valleys and flat areas are represented by relatively homogenous values. Multiple elevation percentile rasters were generated with successively increasing window sizes in order to find a segmentation which preserved the approximate borders of individual landforms without merging them with adjacent ones. An elevation percentile raster generated from a moving window radius of 31 pixels was assessed qualitatively as providing the optimal segmentation. Image segmentation was performed using the mean shift algorithm, which segments a raster into homogenous regions based on a feature space clustering technique (Comaniciu & Meer, 2002). Mean shift segmentation was carried out in Matlab using a wrapper of the original authors' Edge Detection and Image Segmentation system (EDISON) C++ code. The following equation describes the multivariate kernel

84

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

Fig. 3. Subset examples of LSPs extracted from a lidar derived DEM of the proglacial area of Breiðamerkurjökull, Iceland. (A) Elevation percentile, (B) slope, (C) D8 stream network and (D) sinuosity. See Fig. 1A for subset location.

K(x) used in the mean shift calculation in feature space, which is defined by the following inputs: • (hs) the pixel search bandwidth — all pixels within this radius are considered for the mean shift calculation; • (hr) the pixel range bandwidth — all the pixels within (hs) whose values are within the range (hr) are used for calculating the mean Khs hr ðxÞ ¼

C k h2s hpr



2 !
x8




k
hs



r
2 !
x


;
h
r

where, x8 is the spatial part and xr is the range part of a feature vector, k(x) is the common profile used in both domains and C is the corresponding normalisation constant (Comaniciu & Meer, 2002). Mean shift segmentation is a two stage process involving first, an edge preserving mean shift filtering (the basic mean shift calculation) and second, a

concatenation of similar regions into segments based on local gradients (Comaniciu & Meer, 2002). These are defined further below. 1. Specifically, the mean shift calculation works by treating feature space as a probability density function, fixing the kernel described above on each data point and calculating the mean of all the data values within the kernel. The kernel is then ‘shifted’ to the mean data value within feature space and this process is repeated until convergence is achieved (i.e., the kernels no longer move within feature space). Ultimately, the algorithm will converge on dense regions of data points (the pixel values) within feature space, with data points altered based on the trajectory of the kernels, concluding the filtering element of the segmentation process. 2. Segmentation is achieved by joining all homogenous regions within basins of attraction (based on image gradients) associated with mean shift convergence points (Comaniciu & Meer, 2002).

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

Upon completion of the algorithm, the output image consists of homogenous regions, each with a mean shift convergence value. The minimum size of segment (in pixel values) may also be specified (Comaniciu & Meer, 2002). Segmentation parameters (hs, hr) were increased successively from (hs = 3, hr = 3) by increments of 2 until a balance was reached between preserving salient landform boundaries and avoiding over-segmentation. The variable hs defines the “footprint” in which the calculation takes place, leaving hr, which dictates the range of values used in the calculation, to control indirectly the size and shape of resulting segments. Examples of segmentations using different input parameters are displayed in Fig. 4. The segmentation with parameters 10, 7, and 30 was judged qualitatively to provide the best representation of landform boundaries whilst avoiding over-segmentation (Fig. 4C). 4A & 4B are examples of oversegmentation whilst 4D is under-segmented. 3.4. Object Based Image Analysis Each segment was analysed for various pixel-based and geometric properties using algorithms developed in Matlab. The algorithms used

85

to this end are widely implemented in many programming languages for image processing applications (Gonzalez, Woods, & Eddins, 2009). The properties considered in this study are shown in Table 1. Having calculated the various segment properties, they were examined for potential groupings. Colour maps and histograms were used to visualise the segment properties, aiding in the interpretation of potential groupings at each stage of the classification. The classification used a decision tree framework, where object primitives such as ridges were extracted from the segmented data, which were in turn classified into sub-objects, and so on until the target landforms were isolated. At each stage of the decision tree, a clustering algorithm or interactive threshold was used to extract candidate objects. For clustering, the fuzzy c-means algorithm (Bezdek, 1981; Pedrycz, 1990) was used to identify segments (and therefore potential landforms) with similar properties. Unlike hard clustering algorithms, the fuzzy c-means algorithm accounts for the potential soft boundaries between clusters. The fuzzy c-means algorithm requires only two inputs prior to initialisation; the number of clusters, and a fuzzy weighting exponent. An appropriate cluster or pixel threshold value was then chosen to isolate the landform(s) of interest; for example, during classification, the cluster group that corresponded to ridges. Such a cluster group was

Fig. 4. Examples of segmentations using different parameters overlaid on elevation percentile. Parameters are in hs, hr, and minimum segment size in pixels. (A) 3, 3, and 30 (B) 5, 5, and 30, (C) 10, 7, and 30 (D) 15, 10, and 30. See Fig. 1A for subset location.

86

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

Table 1 Selection of segment properties examined in the OBIA process. Pixel intensity refers to the digital number assigned to each pixel. Grey level co-occurrence matrices (GLCMs) are a measure of segment pixel texture in space, based on how often pixel values are adjacent to each other based on a user specified offset. Example offsets could be 0, 45, 90 or 135° from the pixel of interest. The statistical measures of texture are calculated from the resulting GLCM. For further discussion of GLCM based texture measures see Gonzalez et al. (2009). Segment property Basic stats Mean pixel value

Description Mean intensity value of LSP pixels within the segment (i.e., Elevation Percentile, Slope and Sinuosity)

Segment texture measures based on grey level co-occurrence matrix Entropy Randomness of segment pixel values Contrast Intensity contrast based on spatial distribution and intensity values of segment pixels Correlation How correlated pixels are to their respective neighbours within the segment Energy The sum of squared elements within grey-level co-occurrence matrix Homogeneity The closeness of the pixel value distribution within the segment Geometric segment properties Area The number of pixels within the segment Major axis length The major axis of the smallest ellipse enclosing the segment Minor axis length The minor axis of the smallest ellipse enclosing the segment Eccentricity The eccentricity of the smallest ellipse enclosing the segment (the ratio of Minor Axis to Major Axis) Orientation The angle (in degrees) between the horizontal axis of the image and the major axis of the smallest enclosing ellipse of the segment Equivalent diameter The diameter of a circle with the same area as a segment Solidity The proportion of pixels within the convex hull that are also in the segment Extent The ratio of pixels within the segment to its bounding box Centroid position The centre of mass of the segment expressed as a set of coordinates

best defined by segment mean relative elevation values, therefore these values were clustered to produce a map where each cluster corresponded to ridges, flat areas and pits (k = 3). These clusters provided the object primitives for further classification. For thresholding, we adjusted parameters until a single threshold value or multiple values delineated only the landforms of interest. For example, to separate moraine or esker candidates from drumlins, a mean slope threshold N3° was applied. If required, contiguous segments were merged to produce larger segments during analysis, such as contiguous segments that constituted larger landforms such as eskers and moraines, within which minor discontinuities had resulted in over-segmentation. These segments were merged based on their similarity with contiguous neighbours defined as a similarity in mean slope or relative elevation values. Groups of landforms with different morphometric properties were revealed by this combined clustering and thresholding process. For example, moraines tended to have steeper, more asymmetric profiles than flutes or drumlins, and this was reflected in distinct differences in the values of various segment properties, such as the mean slope pixel value. For each landform, the final classified segments were converted to vector format in ArcGIS. The final map of all landforms was then compared with the updated (as described earlier) manually derived map. In this way, the accuracy of the semi-automated mapping technique was assessed. Subsets of the original Evans and Twigg (2002) map, the lidar hillshade, the manually digitised reference map and the final map produced using the semi-automated techniques are displayed in Fig. 5. Accuracy assessment was carried out in Matlab using the confusion matrix method as used in many remote sensing applications, in which the columns represent counts of pixels in the predicted classification

and rows the counts of pixels in the reference map (Lillesand et al., 2008). Individual class accuracies were calculated through two measures; the “producer's accuracy” and “user's accuracy” (Congalton, 1991). Producer's accuracy is the number of correctly classified pixels (columns) divided by the number of reference pixels for the same class (rows), in other words, a measure of omission error. User's accuracy is the total number of pixels correctly classified divided by the total number classified, in other words, a measure of commission error. Overall classification accuracy is calculated as the percentage of pixels correctly classified as compared to the reference map (Congalton, 1991). Additionally, the kappa coefficient was used to measure the level of agreement between the model predictions and reality, as well as whether the results in the confusion matrix are significantly better than random (Lillesand et al., 2008). 4. Results 4.1. Classification process & mapping In this section, we present the results of the semi-automated glacial geomorphological mapping. The processing steps used in the mapping procedure are summarised in Table 2. As explained in the Methods section, the processing steps were, in all cases, performed on the segmented elevation percentile LSP. Landform classes are presented in the order they were extracted from the data. The resultant map is displayed in Fig. 6, with numbered insets providing more detailed views. The landforms (moraines, eskers, flutes/drumlins and glacio-fluvial outwash) were extracted from the same segmentation to avoid overlap between classes. Moraines were classified in five stages by extracting a ridge subclass based on clustering of mean elevation percentile values, thresholding mean slope values within the ridge subclass, followed by localised clustering of orientation values. Eskers were classified in six stages. A mean-slope threshold of 5° separated eskers along with larger moraine ridges. Clustering solidity, followed by clustering of centroid position and localised clustering of orientation, isolated eskers. From the remaining ridge segments, flutes/drumlins were classified by the clustering of centroid position followed by localised clustering of orientation. The remaining low-lying, unclassified segments were used to map glacio-fluvial outwash. Glacio-fluvial outwash was classified using stream/valley network sinuosity as a proxy indicator of glacio-fluvially derived deposits. First, the stream/valley networks were mapped using hydrological techniques developed by Tarboton et al. (1991) and Peucker and Douglas (1975). A D8 flow accumulation raster was computed and a stream/valley network was produced from this using drop analysis (see Tarboton et al., 1991 for a discussion). Second, sinuosity was measured within a 60 m radius moving window, using the stream/valley network raster as the ratio of the actual stream length to the Euclidean distance of the stream, and then mean sinuosity was calculated for each segment. Segments with high sinuosity values were classified as fluvial deposits. Overall, relatively simple measures based on basic pixel statistics (e.g., mean pixel value) and segment geometry (i.e., orientation, position, eccentricity) appeared to be most useful during the classification process at distinguishing between different types of glacial landform. GLCM texture descriptors (e.g., contrast and energy) did not offer any useful distinctions during the classification process. 4.2. Accuracy assessment Table 3 shows the results of the accuracy assessment in confusion matrix format, together with the overall accuracy as a percentage and the kappa statistic. Aside from the accuracy scores, the units represent the number of pixels in the classified and reference datasets for each landform type.

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

87

Fig. 5. Subset examples of (A) Evans and Twigg (2002) map, (B) lidar hillshade, (C) manually digitised reference map, and (D) semi-automated map. See Fig. 1A for subset location. Till is included in C and D to provide visual continuity with A.

For the four landform types considered in this study, an overall classification accuracy of 76.62% was achieved with a kappa coefficient of 0.54. 5. Discussion 5.1. LSPs & segmentation parameters LSPs form an essential part of characterising terrain, as elevation data is insufficient to represent the variability present. Elevation percentile best characterised the variability within the glacial foreland as it represents relative elevation and thus ridges, pits and flat areas, but it required a moving window size to be prescribed. Thus, the first stage in the semi-automated mapping process was the choice of window size in generating the elevation percentile raster for segmentation. We relied on a qualitative approach of trialling multiple moving window scales from a radius of 3 in increments of 2 until the landforms of interest became visually blurred within the raster. However, a more quantitative approach may be desirable in other applications. A possible solution may be to combine multiple window scales of the same LSP and take the average value of each pixel to produce a single band

image with average values. A similar solution may be the application of principal component analysis to multiple window scales of an LSP, whereby the principal component(s) with the dominant Eigen value(s) provides the multi-scale image. Improvement to the mapping techniques in this and previous papers will require a robust empirical analysis of LSPs and generative window sizes to identify optimal parameters for delineating particular types and scales of landforms. During this study, an indirect control on segment size was performed qualitatively on the spatial and range input parameters of the mean shift algorithm. However, a more objective approach would be desirable in future studies. A variable band width mean shift algorithm coded in C++ has been implemented previously by Comaniciu, Ramesh, and Meer (2001) and this could be used within the context of our methodology in the future. Anders et al. (2011) and Drăguţ, Eisank and Strasser (2011) both provide possible solutions to estimating segmentation parameters for an appropriate global scale within the landscape; the former using training samples and the latter using the increase of local variance with scale. These are both based around the input parameters of the multi-resolution segmentation algorithm (Trimble eCognition) which requires a scale parameter input as a control on the average segment heterogeneity within the imagery (Benz, Hofmann,

88

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

Table 2 Summary of segmentation and OBIA processing steps for each landform type. Landform

Simplified processing steps summary

1. Moraine ridges

1. 2. 3. 4. 5. 1. 2. 3. 4. 5. 6. 1. 2. 3. 4. 5. 1. 2. 3.

2. Eskers

3. Flutes/drumlins

4. Glacio-fluvial outwash

Cluster mean elevation percentile to extract ridge class Slope N 3° Cluster centroid to create smaller localised groups Merge segments Cluster orientation for each group Cluster mean elevation percentile to extract ridge class Slope N 5° Cluster centroid to create smaller localised groups Merge segments Cluster solidity Cluster orientation Cluster mean elevation percentile to extract ridge class Slope ≤ 3° Cluster mean elevation percentile Cluster centroid to create smaller localised groups Cluster orientation for each group Exclude previously classified areas from segmented image Calculate channel sinuosity based on D8 stream network Classify segments based on high mean sinuosity value

Willhauck, Lingenfelder, & Heynen, 2004). Given that we use a different segmentation algorithm, we were unable to trial these methods. Additionally, it is worth noting that the input LSP is usually the result of a moving window calculation, and consequently has had a scale

imposed on it by the moving window size. With increasing moving window size, variance in the resulting LSP is reduced. As this will inevitably affect the resulting segmentation, it demonstrates that objective segmentation parameter estimation is not sufficient to estimate scale within the landscape, unless it is applied to the DEM itself. Thus, more focus should be placed on a post-segmentation, iterative classification procedure, where techniques such as segment merging, mathematical morphology, further segmentation and active contours may be used to refine the mapping product. In this case, initial segmentation parameters only need to be approximate. This is particularly important given the variety of software platforms and algorithms available for image segmentation as well as the variety of scales present within remotely sensed data. 5.2. Classification process Moraines were mapped in five stages using a mean slope threshold and fuzzy c-means clustering. Most moraines were classified at the global scale from a simple threshold of mean slope values. A source of error was the apparent identification of drumlins with high slope values as moraines. A localised procedure was needed to eliminate these features from the moraine classification. Furthermore, drumlins which seemingly had resulted from slight modification of moraines by overriding ice were also eliminated based on their local orientations. This suggests the degree of modification by ice could be quantifiable based on

Fig. 6. Classification results for the entire study area superimposed on a hillshade.

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

89

Table 3 Confusion matrix results, overall percentage accuracy & kappa coefficient for the classification accuracy test. Numbers of correctly classified pixels and overall accuracy (%) are in bold. Landform

Drumlins

Moraines

Eskers

Outwash

User's total (pixels)

User's accuracy (%)

Drumlins Moraines Eskers Outwash Producer's total (pixels) Producer's accuracy (%) Overall accuracy

265,146 159,080 6456 248,475 679,157 39.04

33,582 493,168 10,923 174,195 711,868 69.28

1687 57,055 38,460 28,942 126,144 30.49

41,296 206,837 7961 2,403,525 2,659,619 90.37

341,711 916,140 63,800 2,855,137

77.59 53.83 60.28 84.18

76.62

Overall accuracy: 76.62%; kappa: 0.54.

mean slope values. Moraines received the lowest accuracy values, but visual inspection suggests at least parts of most moraines were correctly classified. The lower accuracy is at least partly attributable to only the crests of larger moraine ridges being identified, whereas the manual map also includes the slopes of the moraine. This results in the greater errors of omission and the relatively low errors of commission for moraines (see Table 3). Using multiple scales of elevation percentile may be a means of improving results in future work. Eskers proved the most difficult landform to classify, with the greatest number of steps required to isolate them from other features. Most eskers had greater slope values than most of the moraines, but this global threshold still retained larger moraines such as those deposited during the Little Ice Age maximum. Further steps of classifying features of low solidity, followed by localised orientation, were required to isolate eskers. Many flutes and drumlins were classified correctly, with relatively low errors of omission. However, the flutes and drumlins class had high errors of commission, particularly with moraines and outwash. This highlights a statistical similarity between these features based on the segmentation parameters and classification techniques used. At the outset, orientation appeared to be a powerful descriptor for the separation of flutes/drumlins, but in reality it was only effective at the local scale due to the radial pattern of these landforms in the study area. Thus, to achieve the classification in this study they needed to be first clustered for position, and then clustered for orientation within each position cluster. Glacio-fluvial outwash was classified based on channel sinuosity, as outwash is typically characterised by numerous sinuous river channels. This approach yielded the highest user's and producer's classification accuracies. This indicates that modelled channel network sinuosity is an effective proxy indicator of glacio-fluvial outwash. Glacio-fluvial outwash has a distinctive texture characterised by the braided channels within it; thus a method of isolating this statistically from the other landforms would perhaps improve results further. With all of the landforms, the omission of some features was partially attributable to the segmentation parameters used, due to compromises on size and shape of the resulting segments. The errors of commission for landforms in this study despite using fuzzy techniques suggests that further investigation is required into fuzzy membership parameters for them to be used even more effectively. This ambiguity also raises the possibility that some form of interpretation may always be required to distinguish some landforms from each other, as beyond their context in a glacial landscape little else separates them based on parameters ultimately derived from elevation alone. This underlines the importance of supporting field-based study in ascertaining the identity of landforms and their constituent sedimentology. Improvements in extracting contextual information, such as proximity to different landforms and local orientation, may reduce errors of commission and omission. Beyond lidar elevation derived data, the potential to identify landforms and their constituent materials may also be enhanced by the use of corresponding multi/hyperspectral imagery (Gilvear, Tyler, & Davids, 2004; Marcus, Legleiter, Aspinall, Boardman, & Crabtree, 2003; Schneevoigt et al., 2008).

Whilst the overall classification into object primitives (e.g., ridges, steep ridges, gentle ridges) is quite powerful, it is more difficult to distinguish the specific genetic forms within each object primitive category (e.g., moraines and eskers within the ‘steep ridge’ object primitive). For this, clustering the dataset into localised groups was required, as even other characteristics such as orientation (transverse of parallel to inferred ice flow, for instance) are spatially variable at a scale greater than that of the landforms themselves. 5.3. Comparison with other studies The only study, to our knowledge, that uses semi-automated methods to map glacial landforms specifically is Saha et al. (2011), which maps a drumlin field in North America. Their method consists of mapping the entire drumlin field manually, then using the statistics from their manual map to inform their classification of the drumlins. This would be more informative if the same method was applied to an unmapped area, but is unlikely to be effective in a landscape that contains other ridge-like forms such as eskers and moraines. This study does not offer an accuracy assessment for comparison. Anders et al. (2011) and Schneevoigt et al. (2008) are the next most comparable studies but are not aimed at mapping glacial landforms per se, rather all mappable land forms within their datasets. Nevertheless, our classification methods achieve greater accuracy (+ 5%) than Anders et al. (2011), whose study is based on elevation derived data only. Schneevoigt et al. (2008) achieve greater accuracy (92%), which may be partly attributable to their use of multi-spectral data, but their accuracy results must be treated with caution as their reference dataset only covered small areas of their map. 6. Summary and conclusions We set out to explore the potential of OBIA techniques for mapping glacial geomorphology from fine spatial resolution lidar data. Our methodology consisted of: • the generation of four LSPs (elevation percentile, slope, D8 stream network, stream sinuosity) that were used to represent the characteristics of the target glacial landforms. • image segmentation of elevation percentile. • calculation of segment geometric and underlying LSP properties. • decision tree-based classification using thresholds and fuzzy-c means clustering.

We have demonstrated that various pixel descriptors and object metrics can be used to map glacial geomorphology with relatively high accuracy. We are able to appraise each stage of the classification process and modify it as required, rendering it a more intuitive process akin to human interpretation. Clustering techniques show promise in subdividing and defining landform types from candidate image segments, particularly when multiple thresholds are required or the defining boundaries between landforms are not obvious.

90

C. Robb et al. / Remote Sensing of Environment 163 (2015) 80–90

We mapped moraines, eskers, flutes/drumlins, and glacio-fluvial outwash to an overall accuracy of 77% compared with a map produced manually using traditional methods. Our study demonstrates the potential of OBIA to map large areas of glacial geomorphology quickly and accurately. With further algorithm development on different lidar datasets, we are confident that accuracy results may be improved further. A multiscale LSP for the initial segmentation may aid in finding the appropriate scale at which to segment the landscape in question. Future study should also focus on developing more effective descriptors to differentiate between the more ambiguous landforms (such as moraines and eskers). The inclusion of lidar intensity and multi-spectral data within the classification process in future work may aid in the identification of depositional units and less morphologically distinct landforms. The success of our approach so far, suggests that semi-automated mapping of glacial geomorphic features in other contemporary glacierised settings or former glaciated environments should be viable. This could include places where the initial elevation data are collected at different spatial scales using different platforms, e.g., DEMs derived from satellites or from ship-borne swath bathymetry. Acknowledgements CR acknowledges the financial assistance of the Cambridge Home & European Scholarship Scheme, BB Roberts Fund, Scandinavian Studies Fund and RGS Dudley Stamp Memorial Fund. We also thank Professor David Evans for providing the original map, and Dr. Tomas Jóhannesson, Dr. Finnur Pálsson and Professor Helgi Björnsson for providing the lidar data and for their help and encouragement. References Anders, N.S., Seijmonsbergen, A.C., & Bouten, W. (2011). Segmentation optimization and stratified object-based analysis for semi-automated geomorphological mapping. Remote Sensing of Environment, 115, 2976–2985. Bennett, M.R. (2001). The morphology, structural evolution and significance of push moraines. Earth-Science Reviews, 53, 197–236. Benz, U.C., Hofmann, P., Willhauck, G., Lingenfelder, I., & Heynen, M. (2004). Multi-resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information. ISPRS Journal of Photogrammetry and Remote Sensing, 58, 239–258. Bezdek, J.C. (1981). Pattern recognition with fuzzy objective function algorithms. Kluwer Academic Publishers. Bishop, M.P., James, L.A., Shroder, J.F., Jr., & Walsh, S.J. (2012). Geospatial technologies and digital geomorphological mapping: Concepts, issues and research. Geomorphology, 137, 5–26. Björnsson, H., & Pálsson, F. (2008). Icelandic glaciers. Jökull, 58, 365–386. Blaschke, T. (2010). Object based image analysis for remote sensing. ISPRS Journal of Photogrammetry and Remote Sensing, 65, 2–16. Blaschke, T., Lang, S., & Hay, G. (2008). Object-based image analysis: Spatial concepts for knowledge-driven remote sensing applications. Springer Science & Business Media. Bradwell, T., Stoker, M.S., Golledge, N.R., Wilson, C.K., Merritt, J.W., Long, D., et al. (2008). The northern sector of the last British Ice Sheet: Maximum extent and demise. EarthScience Reviews, 88, 207–226. Clark, C.D., Evans, D.J.A., Khatwa, A., Bradwell, T., Jordan, C.J., Marsh, S.H., et al. (2004). Map and GIS database of glacial landforms and features related to the last British Ice Sheet. Boreas, 33, 359–375. Clark, C.D., Hughes, A.L., Greenwood, S.L., Jordan, C., & Sejrup, H.P. (2012). Pattern and timing of retreat of the last British–Irish Ice Sheet. Quaternary Science Reviews, 44, 112–146. Comaniciu, D., & Meer, P. (2002). Mean shift: A robust approach toward feature space analysis. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 24, 603–619. Comaniciu, D., Ramesh, V., & Meer, P. (2001). The variable bandwidth mean shift and data-driven scale selection. Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on, vol. 431. (pp. 438–445). Congalton, R.G. (1991). A review of assessing the accuracy of classifications of remotely sensed data. Remote Sensing of Environment, 37, 35–46.

Derivaux, S., Forestier, G., Wemmert, C., & Lefèvre, S. (2010). Supervised image segmentation using watershed transform, fuzzy classification and evolutionary computation. Pattern Recognition Letters, 31, 2364–2374. Dilts, T. (2009). Stream gradient and sinuosity toolbox. University of Nevada Reno. Drăguţ, L., & Blaschke, T. (2006). Automated classification of landform elements using object-based image analysis. Geomorphology, 81, 330–344. Drăguţ, L., Eisank, C., & Strasser, T. (2011). Local variance for multi-scale analysis in geomorphometry. Geomorphology, 130, 162–172. Eisank, C., Smith, M., & Hillier, J. (2014). Assessment of multiresolution segmentation for delimiting drumlins in digital elevation models. Geomorphology, 214, 452–464. Evans, D.J.A., Clark, C.D., & Mitchell, W.A. (2005). The last British Ice Sheet: A review of the evidence utilised in the compilation of the Glacial Map of Britain. Earth-Science Reviews, 70, 253–312. Evans, D.J.A., & Twigg, D.R. (2002). The active temperate glacial land system: A model based on Breiðamerkurjökull and Fjallsjökull, Iceland. Quaternary Science Reviews, 21, 2143–2177. Gilvear, D., Tyler, A., & Davids, C. (2004). Detection of estuarine and tidal river hydromorphology using hyper-spectral and lidar data: Forth estuary, Scotland. Estuarine, Coastal and Shelf Science, 61, 379–392. Gonzalez, R.C., Woods, R.E., & Eddins, S.L. (2009). Digital image processing using MATLAB. Prentice-Hall, Inc. Hengl, T., & Reuter, H.I. (2009). Geomorphometry: Concepts, software, applications. Elsevier. Hillier, J.K., & Smith, M.J. (2012). Testing 3D landform quantification methods with synthetic drumlins in a real digital elevation model. Geomorphology, 153–154, 61–73. Institute of Photogrammetry and Remote Sensing (2002). SCOP++. Vienna University of Technology. Iwahashi, J., & Pike, R.J. (2007). Automated classifications of topography from DEMs by an unsupervised nested-means algorithm and a three-part geometric signature. Geomorphology, 86, 409–440. Jóhannesson, T., Björnsson, H., Pálsson, F., Sigurðsson, O., & Þorsteinn (2011). Lidar mapping of the Snæfellsjökull ice cap, western Iceland. Jökull, 61. Lillesand, T.M., Kiefer, R.W., & Chipman, J.W. (2008). Remote sensing and image interpretation. New York: Wiley & Sons. Lindsay, J. (2014). Whitebox GIS. University of Guelph. MacMillan, R.A., Martin, T.C., Earle, T.J., & McNabb, D.H. (2003). Automated analysis and classification of landforms using high-resolution digital elevation data: Applications and issues. Canadian Journal of Remote Sensing, 29, 592–606. Marcus, W.A., Legleiter, C.J., Aspinall, R.J., Boardman, J.W., & Crabtree, R.L. (2003). High spatial resolution hyperspectral mapping of in-stream habitats, depths, and woody debris in mountain streams. Geomorphology, 55, 363–380. Napieralski, J., Harbor, J., & Li, Y. (2007). Glacial geomorphology and geographic information systems. Earth-Science Reviews, 85, 1–22. Ottesen, D., Dowdeswell, J.A., Benn, D.I., Kristensen, L., Christiansen, H.H., Christensen, O., et al. (2008). Submarine landforms characteristic of glacier surges in two Spitsbergen fjords. Quaternary Science Reviews, 27, 1583–1599. Pedrycz, W. (1990). Fuzzy sets in pattern recognition: Methodology and methods. Pattern Recognition, 23, 121–146. Pennock, D.J., Zebarth, B.J., & De Jong, E. (1987). Landform classification and soil distribution in hummocky terrain, Saskatchewan, Canada. Geoderma, 40, 297–315. Peucker, T.K., & Douglas, D.H. (1975). Detection of surface-specific points by local parallel processing of discrete terrain elevation data. Computer Graphics and Image Processing, 4, 375–387. Reuter, H.I., Wendroth, O., & Kersebaum, K.C. (2006). Optimisation of relief classification for different levels of generalisation. Geomorphology, 77, 79–89. Saha, K., Wells, N.A., & Munro-Stasiuk, M. (2011). An object-oriented approach to automated landform mapping: A case study of drumlins. Computers & Geosciences, 37, 1324–1336. Salcher, B.C., Hinsch, R., & Wagreich, M. (2010). High-resolution mapping of glacial landforms in the North Alpine Foreland, Austria. Geomorphology, 122, 283–293. Schneevoigt, N.J., van der Linden, S., Thamm, H.P., & Schrott, L. (2008). Detecting Alpine landforms from remotely sensed imagery. A pilot study in the Bavarian Alps. Geomorphology, 93, 104–119. Stepinski, T.F., & Bagaria, C. (2009). Segmentation-based unsupervised terrain classification for generation of physiographic maps. Geoscience and Remote Sensing Letters, IEEE, 6, 733–737. Tagil, S., & Jennes, J.S. (2008). GIS-based automated landform classification and topographic, landcover and geologic attributes of landforms around the Yazoren Polje, Turkey. Journal of Applied Sciences, 8, 910–921. Tarboton, D.G. (2001). TauDEM (Terrain Analysis Using Digital Elevation Models). Utah State University. Tarboton, D.G., Bras, R.L., & Rodriguez-Iturbe, I. (1991). On the extraction of channel networks from digital elevation data. Hydrological Processes, 5, 81–100. van Asselen, S., & Seijmonsbergen, A.C. (2006). Expert-driven semi-automated geomorphological mapping for a mountainous area using a laser DTM. Geomorphology, 78, 309–320. Woodcock, C.E., & Strahler, A.H. (1987). The factor of scale in remote sensing. Remote Sensing of Environment, 21, 311–332.