- Email: [email protected]
The Cookie Cutter: A method for obtaining a quantitative 3D description of glacial bedforms
Geomorphology 108 (2009) 209–218
Contents lists available at ScienceDirect
Geomorphology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / g e o m o r p h
The Cookie Cutter: A method for obtaining a quantitative 3D description of glacial bedforms M.J. Smith a,⁎, J. Rose b, M.B. Gousie c a b c
School of Earth Sciences and Geography, Kingston University, Kingston-upon-Thames, Surrey, KT1 2EE, UK Department of Geography, Royal Holloway, University of London, Egham, Surrey, TW20 0EX, UK Department of Mathematics and Computer Science, Wheaton College, MA, 02766, USA
a r t i c l e
i n f o
Article history: Received 1 October 2008 Received in revised form 16 January 2009 Accepted 19 January 2009 Available online 27 January 2009 Keywords: DEM Relief Volume Drumlin Glacier bedform
a b s t r a c t Recording the position and attributes of subglacial bedforms, particularly drumlins, is paramount in determining the extent, and dynamics, of former ice sheets. This paper presents a method of deriving 3D properties of glacier bedforms (drumlins in this case) in order to quantify bedform dimensions and acquire information necessary for further investigations, such as calculating amounts of sediment moved subglacially and hence interactions between a glacier and a deformable bed. The method developed here is a semi-automated technique, called Cookie Cutter after the baking implement/process. This method is applied to 5 m spatial resolution DEM data and is based on manual geomorphological mapping of the DEM images, which forms the judgmental part of the process. The mapped bedforms are then processed individually using an automated technique which is described in detail with worked examples from western Central Scotland, which was last glaciated during the Last Glacial Maximum and the Younger Dryas. The advantages and potential sources of error are discussed and results from the sample area are used to compare the volumes of drumlins in an area of LGM glaciation with those in an area glaciated during the Younger Dryas. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Glacier bedforms, such as drumlins, are visually distinctive landforms and have had a long history of study (Close, 1867). Early research mapped drumlin distribution in order to reconstruct the regional palaeo-glaciology (Close, 1867; Charlesworth, 1924). During the latter half of the 20th century, research focussed upon the parameterization of drumlins (Smalley and Unwin, 1968) and their internal composition in an attempt to describe and classify them and provide a basis for understanding their genesis. Subsequently, research into drumlin morphology examined ideas of glacier/bed interactions and bed preservation (Rose and Letzer, 1975; Kleman, 1994) and the landform record of multiple ice flows (Rose and Letzer, 1977; Boulton and Clark, 1990). In parallel with these later developments the concepts of bedform theory (Rose, 1987) and deforming bed processes (Smalley and Unwin, 1968; Boulton, 1987) were introduced to the study of drumlins. The concepts developed from research into drumlin morphology have contributed to the reconstruction of palaeo-ice sheets (e.g. Boulton and Clark, 1990), determination of patterns of glaciation across many areas of former ice cover in temperate latitudes
⁎ Corresponding author. Tel.: +44 207 099 2817; fax: +44 870 063 3061. E-mail address: [email protected] (M.J. Smith). 0169-555X/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.geomorph.2009.01.006
(Kleman et al., 1997) and, in one case, the quantification of the amount of sediment moved during particular glacial events (Rose, 1989). Methods used for the recording of glacier bedform distribution and morphology have also developed considerably, enabling, in part, the gathering of larger data sets over shorter time scales and in greater detail (Clark, 1997). However, 3D quantification has received little attention despite the importance of this property for determining the relief and volume of bedforms. Early quantitative work focused solely on twodimensional (2D) parameters, partly because of the inadequacy of elevational data (Rose and Letzer, 1975). 3D properties are a fundamental requirement for quantifying the amount of sediment moved by glaciers beneath an active ice sheet or glacier and understanding glacier dynamics in terms of work done at a deformable glacier bed. Barriers to calculating bedform volumes are being broken down with the increasing availability of high spatial resolution surface elevation data (Smith et al., 2006), primarily in the form of digital elevation models (DEMs). This issue was addressed by Evans (1987) who noted the trend towards the production of regional, detailed, digital elevation data and illustrated how it could be used to calculate morphometric parameters of drumlins. In particular he described the potential to outline drumlins on DEMs and subsequently characterize their 3D form. Here we develop a methodology for characterizing the relief and volume of drumlins and demonstrate the application of the method with worked case studies. Results are discussed in conjunction with potential sources of error.
210
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
2. Methodology In this paper we propose a three stage procedure for the calculation of landform volume (Fig. 1), as a pre-requisite to further geodynamic studies. Stage 1 involves initial mapping to identify and outline landforms; Stage 2 defines landform relief and Stage 3 calculates individual landform volumes. Stage 1 requires scientific judgement; Stages 2 and 3 are automated. For the purposes of this paper this methodology is applied to a Late Quaternary glaciated region with a drumlin landscape. For this study, the NEXTMap Britain™ DEM (referred to as NEXTMap) is used as the data source for geomorphological mapping and subsequent calculations of drumlin relief and volume. NEXTMap images have been tested against field mapped data (Smith et al., 2006) and found to be the highest quality, readily available data source. NEXTMap is a single-pass interferometric synthetic aperture radar (IfSAR) product, with a spatial resolution of 5 m and a vertical accuracy of 0.5–1 m (Intermap, 2005). The DEM is supplied as a raster grid, with each cell within the grid assigned a value that represents the average elevation of the land surface. Fig. 1 outlines the conceptual and technical process used in the derivation of drumlin volume. All processing is performed within a geographic information system (GIS). This requires an input DEM (Fig. 1, I) which is then used to visualize the drumlized terrain (Fig. 1, II). The outlines of the drumlins are manually digitized from the DEM (Fig. 1, III). The area with the outline is then removed (the “cookie cutting” process) to produce a new DEM of “holes” or voids (Fig. 1, IV). A planar surface is then interpolated across each of the individual voids (Fig. 1, V) to leave a new “infilled” DEM. This is then subtracted from the
original DEM (Fig. 1, VI), leaving a DEM of drumlin relief. Volume is calculated using drumlin area and height (Fig. 1, VII). 2.1. Study area With the primary requirement being a study area where a high spatial resolution DEM can be validated against geomorphological ground control, we selected sites from the western part of central Scotland (Fig. 2; Smith et al., 2006). The landforms studied are drumlins formed during the Last Glacial Maximum (LGM), when ice moved west to east across the region, and during the Younger Dryas (YD) glaciation, when ice flow was predominantly from north to south (Fig. 2). Thus, although the principles outlined in this paper are applicable to all glacier bedforms (flutes, megaflutes, drumlins, megadrumlins, streamlined hills and Rogen) the only examples studied are drumlins. Henceforward the term drumlin is substituted for bedform. The full study area covers ~750 km2 comprising a total of 2076 drumlins mapped in the field at 1:10,000 scale (Rose and Smith, 2008) of which 214 are from the area of YD glaciation which covers ~104 km2 (excluding Loch Lomond and the islands). For the purpose of demonstrating the methodology of this paper, two worked examples have been taken from the LGM glaciated region. Additionally a sample of 175 drumlins has been studied from an area that was in part last glaciated during the LGM and in part during the YD [note that the size of this sample differs from that (189) in Smith et al. (2006) because that reports all of the drumlins in the area, including parts of drumlins, whereas the sample used in this study only deals with complete drumlins]. All the drumlins used for these tests are identified from the NEXTMap DEM and the area used in the large sample is the same as that used for the comparison of field mapped results with NEXTMap results in Smith et al. (2006). The worked examples do not discuss the glaciodynamic implications of volumetric calculations; as it is necessary to understand the internal composition of the bedform before this process can be carried out, and this is an entirely different issue from that considered here. All figures depicted in map coordinates are projected to the British National Grid, with north–south indicated by vertical grid lines and north implicitly to the top of the figure. 2.2. Determination of bedform area Identification and representation of the drumlin outline requires the optimal visualization of the DEM. Smith and Clark (2005) suggest using images of local contrast stretch and gradient, supplemented with two orthogonal relief shaded images. The drumlin outline is then identified based upon breaks-of-slope and recorded as a vector polygon through on-screen digitizing. It is important to use visualization techniques that do not incorporate biases to the orientation of the drumlin (Smith and Wise, 2007), although Smith and Clark (2005) recommend the use of relief shaded images for the identification of topographically subtle drumlins. This process requires scientific judgement with volumetric calculations dependent upon the quality of digitization. Quantitative validation of the DEM mapping performed here is reported in Smith et al. (2006). 2.3. The Cookie Cutter
Fig. 1. Flow chart outlining the routine for the calculation of glacier bedform volume. This process requires visualization of a DEM (I and II) and the digitizing of drumlin outlines (III; a subjective operator process) which are then used to cookie cut voids from the original DEM (IV). The voids are then infilled (V) and the resultant data set is subtracted from the original DEM (VI) to produce estimates of relief. Volume can then be calculated (VII). IV to VII are performed on individual drumlins.
Cookie Cutter refers to the process whereby a drumlin is extracted from the original DEM (Fig. 1, IV). The term is derived from a cookie cutter used to cut shapes from pastry; here we cut drumlins from the DEMs. The cookie cutter cuts a single defined shape for each drumlin individually in order that each feature is treated separately. This is necessary in order to overcome problems associated with drumlin adjacency (where two drumlins are in contact), a normal feature of glacier bedform terrain (Fig. 3a), and is achieved by adopting a buffering process. If the drumlins were processed in a single operation, the adjacent rising slopes would influence the interpolation of the basal surface.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
211
Fig. 2. Location of the field mapping area (~ 750 km2) and comparative study area (~104 km2). The inset map (top left) shows generalized ice flow during the LGM (black arrows after Sissons, 1967) and YD (red outline after Rose, 1987). Map coordinates are projected to the British National Grid, with north–south indicated by vertical grid lines and north implicitly to the top of the figure.
The buffering process is only possible when the cookie cutter is performed as a separate operation for each drumlin (Fig. 3b). The process involves conversion of the area within the drumlin outline to a raster image comprised of a null value which is ignored in further calculations. This forms the void that is filled at the interpolation stage. The buffer used to separate adjacent drumlins can be user-specified, and in this instance we have defined a zone of 20 m around the drumlin outline. A conditional
function within the GIS is used to process each drumlin individually, on a cell-by-cell basis. Any cells in the buffered region are replaced by an elevation value derived from the level of the boundary of the DEM, whilst any cells within the drumlin are replaced with the null value. The output of this process is a drumlin void surrounded by a 20 m buffer with surface elevation identical to that of the drumlin boundary. This process is performed on every drumlin to produce DEMs with drumlin voids.
Fig. 3. Illustration of drumlin adjacency. (a) Digitized outlines of two drumlins that are adjacent. (b) The cookie cutter requires a buffer to be created around each drumlin (dashed line), a process not possible with a single operation as the buffered regions overlap. This stage therefore requires each drumlin to be processed individually.
212
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Fig. 4. Examples of drumlins on different topographic surfaces. (a) Diagrammatic representation of drumlin relief (relative elevation) on a planar horizontal surface depicting the difference in elevation from surrounding terrain, (b) drumlin relief on a planar tilted surface (e.g. hillside) and (c, d) depiction of drumlins formed on curving terrain configurations.
2.4. Determination of basal surface For the purpose of this research, the relief (relative elevation) of a drumlin is defined as the difference in elevation between an individual point on a drumlin surface and a plane interpolated across a drumlin outline which is “representative” of the base of the
drumlin (Figs. 1, V, 4a). This can be a simple relationship when the representative base is planar and level (Fig. 4a), but becomes much more complex when the representative base is inclined (Fig. 4b) or complex (Fig. 4c,d). The interpolated plane creates a basal surface for the drumlin void produced during the cookie cutting stage. It therefore prescribes the datum that is used for the calculation of relief
Fig. 5. Example of basal surface interpolation for a single drumlin. (a) relief shaded DEM, with digitized outline, (b) the blank void created after the removal of the drumlin, (c) the drumlin basal surface interpolated using IDW and (d) the drumlin basal surface interpolated using a tensioned thin-plate spline.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Fig. 6. Perspective view of a drumlin overlaid on a DEM grid. The interpolated base is shaded, with relief for each individual cell calculated from the base to the drumlin surface.
on a point-by-point basis (i.e. the elevation of the void will be different at every location). Relief is then calculated for all cells on the drumlin above this surface (Fig. 5) and these values are summed to cover the planform area and thus compute volume. Note that accurate calculations of bedform relief make assumptions about the terrain prior to the formation of the drumlin. This is a fundamental geomorphological problem that will be considered in a later paper; here we are concerned only with demonstrating and applying the cookie cutter technique. Although the above is simple in concept, calculation of relief can be difficult as bedforms can be found in complex terrain configurations, including sloping terrain and across ridges and depressions (Fig. 4b,c, d; Rose and Letzer, 1975, 1977; Mitchell, 1994). The interpolated surface needs to slope in order to take these variations into account and, therefore, is likely to be complex. A number of commonly used interpolation techniques were considered: inverse distance weighting (IDW), splines and kriging (Dubrule, 1984; Oliver and Webster, 1990). Kriging was discounted as, being a geostatistical technique, it ideally needs to be optimized for each application (Isaaks and Srivastava, 1990). Therefore both inverse distance weighting and splines were tested as to their suitability for this application. Inverse distance weighting is a simple and fast method in which an interpolated point is computed, based on the weighted average of its neighbours, where the weight of a neighbour is determined by its
213
distance from the point (Heine, 1986; Watson, 1992). The problem with using such a simple interpolator is that an artefact at the edge of the drumlin void (e.g. a hedge or man-made object) will, in effect, be “pulled” into the area of interpolation. This is illustrated in Fig. 5 which shows a drumlin crossed by a hedge (Fig. 5a). The void created by the cookie cutter is shown on Fig. 5b and the effect of the hedge is evident when the IDW interpolator is extended across this void. This hedge creates error because artifacts on a drumlin contribute to the definition of the basal surface (Fig. 5c). This problem can be resolved to some degree by making the area of interpolation as planar as possible with the use of a thin-plate spline under tension (Smith and Wessel, 1990; Fig. 5d). The thin-plate spline under tension can be described, through analogy, as a membrane that is pulled taught to the edges of the drumlin void; the force can be adjusted by a tension parameter. In this case, artifacts on the edges of the drumlin are smoothed, leaving the basal surface with increased planarity (Fig. 5d). The tension parameter is defined between 0 and 1, where 0 has no tension and means the interpolated surface conforms to the minimum curvature and 1 has infinite tension and maximum curvature. Infinite tension is not physically meaningful (Smith and Wessel, 1990), but implies that any further increase of tension cannot result in a smoother surface while remaining true to the original data points. Infinite tension was chosen to minimize the effects of artifacts whilst conforming to the surrounding terrain. As can be seen in Fig. 5 the IDW method produces a relatively poor result as the hedges can be seen to “smear” across the drumlin area. The thin-plate spline performs well, defining a surface around the base of the drumlin. The thin-plate spline was therefore selected for the interpolation of surfaces across all drumlin voids. 2.5. Determination of bedform relief and volume After the interpolation of the drumlin basal surface, relief can be calculated by subtracting the interpolated basal values (excluding the buffer) from the original DEM comprised of drumlin elevations (Fig. 1, VI). Each cell in the output dataset contains an elevation value representing the height of the drumlin above its basal surface (Fig. 6). When this value is multiplied by the area of each cell (25 m2 in this instance), an estimate of volume of the cell column is produced. The summation of the volumes of all cell columns that comprise a drumlin produces a final total volume for the individual bedform.
Fig. 7. Flow diagram showing the Cookie Cutter computer processing stages (circles represent input parameters). (a) Visualisation of an input DEM is used as a data source to digitize drumlin outlines, (b) a project directory is specified prior to running the “Cookie Setup” script and (c) separated drumlin outlines and the original DEM form the primary inputs to the “Cookie Cutter” script. Four input parameters are required and, once complete, the script outputs drumlin relief and summary statistics.
214
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
3. Worked examples on small samples The previous sections have outlined the principle steps required to estimate drumlin volume. This section outlines the procedure in detail and then provides two worked examples illustrating how this method can be applied. 3.1. Procedure (Fig. 7) The procedure for carrying out this process is as follows (Fig. 7): i) Digitization of drumlin outlines is performed manually within ESRI ArcMap and saved as a polygon shapefile (Fig. 7a) ii) The cookie-cutting process is automated using two custom ArcGIS Python scripts which are loaded and run from ArcToolBox. The Python scripts are freely available for download as ancillary materials accompanying this paper. The “Cookie Setup” script (Fig. 7b) separates each digitized drumlin outline into a separate file, as well as establishing a series of directories for computer processing. Input in to the script is a shapefile of digitized drumlin outlines, with the user required to specify the project directory. The “Cookie Cutter” script (Fig. 7c) buffers the original drumlin outline in order to produce the drumlin void independent of
interference from other drumlins. The spline interpolator within ArcGIS fills the void before the production of the final output. A raster file represents the relief for each drumlin, supplemented with a database collating summary statistics for each drumlin. Input into the script is a DEM, with drumlin outlines automatically selected from the output of the “Cookie Setup” script. The user is required to specify the project folder, cell size (it is recommended to use the cell size of the input DEM), tension parameter for the spline interpolator and distance that each drumlin will be buffered (at least four cells). 3.2. Bowridge (Fig. 8) Fig. 8a depicts a well defined megadrumlin with three superimposed drumlins that are in the order of hundreds of metres long and tens of metres high. The location is shown on Fig. 2; the hill is known as Bowridge and is centred on 1 km grid square NS 7880. Fig. 8b shows a standard greyscale visualization of the DEM and demonstrates the ineffectiveness of this type of image. Fig. 8c shows the most effective visualization of the drumlins after rotating light source directions on the NEXTMap image and the boundaries of these drumlins are outlined. The drumlin outlines were buffered and then used to cookie cut the DEM to create a series of drumlin voids (Fig. 8d). A tensioned
Fig. 8. Worked example demonstrating the steps required to calculate drumlin volume —(a) the mega drumlin located at Bowridge (NS 7880) (see Fig. 2 for location), (b) raw DEM data, (c) relief shaded visualization of terrain, with mapped drumlin outlines, (d) DEM voids, (e) interpolation of drumlin basal surfaces and (f) relief shaded visualization of drumlin volumes (1.51 × 106 m3). Note the “stepping” in (e), a result of artifacts at the edges being interpolated across the basal surface.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
215
Fig. 9. Worked example demonstrating the steps required to calculate drumlin volume — the drumlin located near Milton of Campsie (NS 6677) (see Fig. 2 for location). (a) original field mapping, (b) raw DEM data, (c) relief shaded visualization of terrain, with mapped drumlin outlines, (d) DEM voids, (e) interpolation of drumlin basal surfaces and (f) relief shaded visualization of drumlin volumes (2.07 × 106 m3). Note the “stepping” in (e), a result of artifacts at the edges being interpolated across the basal surface. Note also the difference between the interpretations of drumlin form and extent derived from field mapping and analysis of NEXTMap showing the scope of the two methods. This difference does not affect the calculations but illustrates the need for ground truthing whenever possible (Rose and Smith, 2008).
spline was then used to interpolate the basal surface of the drumlins (Fig. 8e; note the effect of edge artifacts on the smoothness of the basal surfaces). This basal surface is subtracted from the original drumlins to enable a calculation of drumlin volumes (Fig. 8f; note that this final image has been relief shaded to help visualization).
smoothness of the basal surfaces) and the final drumlin volumes were calculated (Fig. 9f).
3.3. Milton of Campsie (Fig. 9)
In addition to developing a technique for the estimation of drumlin volume, or the volume of any other landform, we have applied the method to a small area from the western part of central Scotland that includes drumlins formed during both the LGM and the YD glaciations. This area covers ˜ 104 km2 and contains 175 drumlins mapped from the NEXTMap DEM, ranging from 88 to 844 m in length (see Smith et al., 2006 for further details). The drumlins formed during the LGM and the YD glaciations are identified on Fig. 10. Table 1 presents drumlin relief and drumlin volume results for the entire study area and, separately, the sub-samples for drumlins formed during the LGM and YD. The minimum drumlin relief value of
The second worked example was selected as the drumlins are larger and clearly arranged adjacent to one another, with shared outline boundaries. The field mapping (Fig. 9a) and analysis of the DEM (Fig. 9c) depict the same general arrangement of drumlins and, as at Bowridge, the greyscale image is ineffective (Fig. 9b). The drumlin outlines are depicted on the NEXTMap image following analysis of the NEXTMap visualizations (Fig. 9c). Fig. 9d shows the cookie cut drumlin outlines; basal surfaces were subsequently interpolated (Fig. 9e; again note the effect of edge artifacts on the
4. Comparing the volume of drumlins formed during LGM and YD glaciations
216
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Fig. 10. Comparative study area depicting drumlins mapped from NEXTMap Britain DEM data. The dotted line demarcates the boundary between LGM and YD drumlins. Map coordinates are projected to the British National Grid, with north–south indicated by vertical grid lines and north implicitly to the top of the figure.
−16 m (and minimum drumlin volume of −1.15 × 103 m3) indicates that the results must contain error as it is not possible to have negative relief values. Investigation of the sources of these errors indicated eight negative values that resulted from the creation of “sliver polygons” (Longley et al., 2005) (Fig. 11a) or “clutter” (such as trees or hedges) at the drumlin edges creating the effect of negative relief when the basal surface is interpolated (Fig. 11b). Sliver polygons result from incorrect digitization and can simply be deleted. Drumlins with negative relief values resulting from clutter were considered erroneous and excluded from further analysis. The summary statistics (Table 1) indicate that relative to the LGM sub-sample of drumlins the volumes for the YD sub-sample give a smaller mean size with a larger
range and larger standard deviation. A t-test comparing drumlin volumes for the two sub-samples, with equal variances not assumed, was statistically significant (0.012) indicating that the datasets are different. This fact is reported here as an outcome from this exercise,
Table 1 Volume calculations for the study area (figures in bold are corrected after the removal of negative relief values), incorporating the two sub-samples for drumlins formed during the YD and LGM. Measure
Study area (n = 175)
Relief min (m) Relief max (m) Relief mean (m) Relief standard deviation (m) Volume min (106 m3) Volume max (106 m3) Volume mean (106 m3) Volume standard deviation (m) Total volume (106 m3)
− 16.3 65 3.7 1.8 − 0.11 2.75 0.12 0.24 22.75
1.7 65 9.3 7.0 0.0016 2.75 0.13 0.24 23.00
Younger Dryas (n = 104)
Last Glacial Maximum (n = 71)
1.65 65 8.1 7.2
2.41 31 11.1 6.2
0.0016 0.70 0.12 0.29 12.31
0.0031 2.75 0.15 0.16 10.69
Fig. 11. Two common errors associated with automated calculation of drumlin volumes. (a) Mis-digitization can create erroneous “sliver” polygons. (b) The presence of edge artifacts can lead to the interpolation of basal surfaces not related to the drumlin topography. The dotted line shows the drumlin basal surface, whilst the dashed line shows the actual interpolated surface as a result of the tree increasingly the elevation at the edge of the drumlin.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
but the results may identify a critical difference between drumlin samples formed by glaciers with different dimensions, and may confirm the critical importance of this property in controlling bedform dimensions, as proposed in Rose (1989). 5. Discussion The Cookie Cutter is a straightforward approach to calculate the volume of drumlins, other glacier bedforms and well defined positiverelief landforms in general. However, the derived values presented here are only a first approximation due to several potential sources of error that introduce uncertainty into the results. In addition to the errors inherent in the method of collecting the data (quality of the DEM and quality of the geomorphological interpretation; Smith et al., 2006), further sources of uncertainty are outlined below: a) Quality of digitization of interpreted data — this can mean that terrain may be inappropriately included or excluded from the calculation. For example, pronounced furrows can occur around the stoss end of drumlins and along their sides, erroneously incorporating volume in to the calculations. b) DEM vertical accuracy — the accuracy of the measurement of surface elevation will affect the calculation of drumlin volume. The vertical accuracy of NEXTMap is ~0.5–1.0 m. c) Surface Clutter — the NEXTMap DEM used in this study quantifies visible surface elevation, including buildings and vegetation. These features are not relevant to the volume of drumlins and it is desirable that they are removed through an automated procedure (e.g. Wang et al., 2001; Sithole and Vosselman, 2004), although this may have the detrimental effect of removing subtle topographic features, including, for instance small drumlins. Removal through interactive editing of the DEM is also possible, although this can be a time consuming task. Buildings and forestry cover only small parts of the study sites, but hedges are an important element of the landscape and caused problems throughout. Of particular importance is the presence of clutter at the edges of drumlins which can substantially affect drumlin volume estimates and, in some instances, lead to calculations of negative volume. Calculations presented for the study area above are therefore considered to be underestimates, although it should be noted that the presence of clutter on drumlins will cause overestimation. d) Interpolation — the drumlin basal elevation values are estimated through the use of an interpolator and therefore a degree of uncertainty will be introduced as these represent estimates of an unknown, prior surface. In addition to the methodological problems listed above, it should be noted that in-situ landform degradation will have occurred. The drumlins examined in this study were exposed to subaerial process as the ice retreated some c. 15 ka BP and c. 11.8 ka BP following the LGM and YD glaciations respectively. Mass wastage processes, therefore, will have degraded the glaciogenic landforms to varying degrees (Rose, 1991). The combination of the above factors will contribute to error in the calculation of drumlin volume. Potential mitigation measures include the following: i) the use of LiDAR data with a higher spatial resolution, higher vertical accuracy and ability to penetrate vegetation canopies and b) a prescriptive procedure for digitization which will improve consistent mapping. 6. Conclusions This paper presents a technique for the calculation of material volumes of drumlins from sites in western central Scotland. The method can be applied to other glacial bedforms and indeed to any landform. Conclusions from this work are:
217
• Digitized drumlin outlines are used to extract (“cookie cut”) landforms from an underlying DEM leaving empty voids. • A void is “infilled” through the application of a tensioned spline interpolator thereby estimating the basal surface. • Drumlin relief is calculated by subtracting the basal surface from the original DEM and then converted to volume by multiplying by the planform area. • Landform adjacency requires drumlins to be processed individually (to avoid the creation of adjacent voids). The drumlin outlines are therefore buffered (to include surrounding terrain) prior to being “cookie cut”. • A protocol is established, using the above methodology, for the calculation of material volumes for similarly wholly concave or wholly convex landforms. • The results provide a first-order approximation that facilitates the rapid calculation of drumlin volumes and, in the example given in this paper, the differentiation of the volume of drumlins formed by an ice sheet during the LGM, and drumlins formed at the margin of the piedmont lobe of the Loch Lomond outlet glacier during the YD. • Volume calculations are subject to error from the source DEM, digitization procedure and presence of surface clutter, in addition to geomorphological problems inherent in the sample studied. • This technique can potentially be used wherever landform volume estimates are required, for example, calculations of sand dune volumes or sediment input into depositional environments.
Acknowledgements We thank Steve Booth and James O'Brien for their assistance with this work. Reviews from Ian Evans and John Shaw helped improve this manuscript. The NEXTMap DEM data was kindly supplied by the British Geological Survey. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.geomorph.2009.01.006. References Boulton, G.S., 1987. A theory of drumlin formation by sub-glacial sediment deformation. In: Menzies, J., Rose, J. (Eds.), Drumlin Symposium. Balkema, Rotterdam, pp. 25–80. Boulton, G.S., Clark, C.D., 1990. A highly mobile Laurentide Ice Sheet revealed by satellite images of glacial lineations. Nature 346, 813–817. Charlesworth, J.K., 1924. The glacial geology of the north-west of Ireland. Proceedings of the Royal Irish Academy, vol 36B, pp. 174–314. Clark, C.D., 1997. Reconstructing the evolutionary dynamics of palaeo-ice sheets using multi-temporal evidence, remote sensing and GIS. Quaternary Science Reviews 16, 1067–1092. Close, M.H., 1867. Notes on the general glaciation of Ireland. Journal of the Royal Geographical Society of London 1, 207–242. Dubrule, O., 1984. Comparing splines and kriging. Computers and Geosciences 10, 327–338. Evans, I.S., 1987. A new approach to drumlin morphometry. In: Menzies, J., Rose, J. (Eds.), Drumlin Symposium. Balkema, Rotterdam, pp. 119–130. Heine, G.W., 1986. A controlled study of some two-dimensional interpolation methods. Computer Oriented Geological Society Computer Contributions 2, 60–72. Isaaks, E.H., Srivastava, R.M., 1990. An Introduction to Applied Geostatistics. Oxford University Press, New York. Intermap, 2005. Product Handbook [online]. Available from: http://www.intermap. com/images/handbook/producthandbook.pdf [Accessed: 21st October 2005]. Kleman, J., 1994. Preservation of landforms under ice sheets and ice caps. Geomorphology 9, 19–32. Kleman, J., Hätterstrand, C., Borgström, I., Stroeven, A., 1997. Fennoscandian palaeoglaciology reconstructed using a glacial geological inversion model. Journal of Glaciology 43, 283–299. Longley, P.A., Goodchild, M.F., Maguire, D.J., Rhind, D.W., 2005. Geographical Information Systems and Science. Wiley, Chichester. Mitchell, W.A., 1994. Drumlins in ice sheet reconstructions, with reference to the western Pennines, northern England. Sedimentary Geology 91, 313–331. Oliver, M.A., Webster, R., 1990. Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Information Systems 4, 313–332.
218
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Rose, J., 1987. Drumlins as part of a glacier bedform continuum. In: Menzies, J., Rose, J. (Eds.), Drumlin Symposium. Balkema, Rotterdam, pp. 103–116. Rose, J., 1989. Glacial stress patterns and sediment transfer associated with the formation of superimposed flutes. Sedimentary Geology 62, 151–176. Rose, J., 1991. Subaerial modification of glacier bedforms immediately following ice wastage. Norsk Geografisk Tidsskrift 45, 143–153. Rose, J., Letzer, J., 1975. Drumlin measurements: a test of the reliability of data derived from 1:25,000 scale topographic maps. Geological Magazine 112, 361–371. Rose, J., Letzer, J., 1977. Superimposed drumlins. Journal of Glaciology 18, 471–480. Rose, J., Smith, M.J., 2008. Glacial geomorphological maps of the Glasgow region, western central Scotland. Journal of Maps v2008, 399–416. Sissons, J.B., 1967. The Evolution of Scotland's Scenery. Oliver and Boyd, Edinburgh. Sithole, G., Vosselman, G., 2004. Experimental comparison of filter algorithms for bareEarth extraction from airborne laser scanning point clouds. ISPRS Journal of Photogrammetry & Remote Sensing 59, 85–101. Smalley, I.J., Unwin, D.J., 1968. The formation and shape of drumlins and their distribution and orientation in drumlin fields. Journal of Glaciology 7, 377–390.
Smith, W.H.F., Wessel, P., 1990. Gridding with continuous curvature splines in tension. Geophysics 55, 293–305. Smith, M.J., Clark, C.D., 2005. Methods for the visualisation of digital elevation models for landform mapping. Earth Surface Processes and Landforms 30, 885–900. Smith, M.J., Wise, S.M., 2007. Mapping glacial lineaments from satellite imagery: an assessment of the problems and development of best procedure. International Journal of Applied Earth Observation and Geoinformation 9, 65–78. Smith, M.J., Rose, J., Booth, S., 2006. Geomorphological mapping of glacial landforms from remotely sensed data: an evaluation of the principal data sources and an assessment of their quality. Geomorphology 76, 148–165. Wang, Y., Mercer, B., Tao, V.C., Sharma, J., Crawford, S., 2001. Automatic generation of bald earth digital elevation models from digital surface models created using airborne IFSAR. Proceedings of 2001 ASPRS Annual Conference, St. Louis, Missouri, USA. Watson, D., 1992. Contouring: A Guide to the Analysis and Display of Spatial Data. Pergamon Press.
Contents lists available at ScienceDirect
Geomorphology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / g e o m o r p h
The Cookie Cutter: A method for obtaining a quantitative 3D description of glacial bedforms M.J. Smith a,⁎, J. Rose b, M.B. Gousie c a b c
School of Earth Sciences and Geography, Kingston University, Kingston-upon-Thames, Surrey, KT1 2EE, UK Department of Geography, Royal Holloway, University of London, Egham, Surrey, TW20 0EX, UK Department of Mathematics and Computer Science, Wheaton College, MA, 02766, USA
a r t i c l e
i n f o
Article history: Received 1 October 2008 Received in revised form 16 January 2009 Accepted 19 January 2009 Available online 27 January 2009 Keywords: DEM Relief Volume Drumlin Glacier bedform
a b s t r a c t Recording the position and attributes of subglacial bedforms, particularly drumlins, is paramount in determining the extent, and dynamics, of former ice sheets. This paper presents a method of deriving 3D properties of glacier bedforms (drumlins in this case) in order to quantify bedform dimensions and acquire information necessary for further investigations, such as calculating amounts of sediment moved subglacially and hence interactions between a glacier and a deformable bed. The method developed here is a semi-automated technique, called Cookie Cutter after the baking implement/process. This method is applied to 5 m spatial resolution DEM data and is based on manual geomorphological mapping of the DEM images, which forms the judgmental part of the process. The mapped bedforms are then processed individually using an automated technique which is described in detail with worked examples from western Central Scotland, which was last glaciated during the Last Glacial Maximum and the Younger Dryas. The advantages and potential sources of error are discussed and results from the sample area are used to compare the volumes of drumlins in an area of LGM glaciation with those in an area glaciated during the Younger Dryas. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Glacier bedforms, such as drumlins, are visually distinctive landforms and have had a long history of study (Close, 1867). Early research mapped drumlin distribution in order to reconstruct the regional palaeo-glaciology (Close, 1867; Charlesworth, 1924). During the latter half of the 20th century, research focussed upon the parameterization of drumlins (Smalley and Unwin, 1968) and their internal composition in an attempt to describe and classify them and provide a basis for understanding their genesis. Subsequently, research into drumlin morphology examined ideas of glacier/bed interactions and bed preservation (Rose and Letzer, 1975; Kleman, 1994) and the landform record of multiple ice flows (Rose and Letzer, 1977; Boulton and Clark, 1990). In parallel with these later developments the concepts of bedform theory (Rose, 1987) and deforming bed processes (Smalley and Unwin, 1968; Boulton, 1987) were introduced to the study of drumlins. The concepts developed from research into drumlin morphology have contributed to the reconstruction of palaeo-ice sheets (e.g. Boulton and Clark, 1990), determination of patterns of glaciation across many areas of former ice cover in temperate latitudes
⁎ Corresponding author. Tel.: +44 207 099 2817; fax: +44 870 063 3061. E-mail address: [email protected] (M.J. Smith). 0169-555X/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.geomorph.2009.01.006
(Kleman et al., 1997) and, in one case, the quantification of the amount of sediment moved during particular glacial events (Rose, 1989). Methods used for the recording of glacier bedform distribution and morphology have also developed considerably, enabling, in part, the gathering of larger data sets over shorter time scales and in greater detail (Clark, 1997). However, 3D quantification has received little attention despite the importance of this property for determining the relief and volume of bedforms. Early quantitative work focused solely on twodimensional (2D) parameters, partly because of the inadequacy of elevational data (Rose and Letzer, 1975). 3D properties are a fundamental requirement for quantifying the amount of sediment moved by glaciers beneath an active ice sheet or glacier and understanding glacier dynamics in terms of work done at a deformable glacier bed. Barriers to calculating bedform volumes are being broken down with the increasing availability of high spatial resolution surface elevation data (Smith et al., 2006), primarily in the form of digital elevation models (DEMs). This issue was addressed by Evans (1987) who noted the trend towards the production of regional, detailed, digital elevation data and illustrated how it could be used to calculate morphometric parameters of drumlins. In particular he described the potential to outline drumlins on DEMs and subsequently characterize their 3D form. Here we develop a methodology for characterizing the relief and volume of drumlins and demonstrate the application of the method with worked case studies. Results are discussed in conjunction with potential sources of error.
210
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
2. Methodology In this paper we propose a three stage procedure for the calculation of landform volume (Fig. 1), as a pre-requisite to further geodynamic studies. Stage 1 involves initial mapping to identify and outline landforms; Stage 2 defines landform relief and Stage 3 calculates individual landform volumes. Stage 1 requires scientific judgement; Stages 2 and 3 are automated. For the purposes of this paper this methodology is applied to a Late Quaternary glaciated region with a drumlin landscape. For this study, the NEXTMap Britain™ DEM (referred to as NEXTMap) is used as the data source for geomorphological mapping and subsequent calculations of drumlin relief and volume. NEXTMap images have been tested against field mapped data (Smith et al., 2006) and found to be the highest quality, readily available data source. NEXTMap is a single-pass interferometric synthetic aperture radar (IfSAR) product, with a spatial resolution of 5 m and a vertical accuracy of 0.5–1 m (Intermap, 2005). The DEM is supplied as a raster grid, with each cell within the grid assigned a value that represents the average elevation of the land surface. Fig. 1 outlines the conceptual and technical process used in the derivation of drumlin volume. All processing is performed within a geographic information system (GIS). This requires an input DEM (Fig. 1, I) which is then used to visualize the drumlized terrain (Fig. 1, II). The outlines of the drumlins are manually digitized from the DEM (Fig. 1, III). The area with the outline is then removed (the “cookie cutting” process) to produce a new DEM of “holes” or voids (Fig. 1, IV). A planar surface is then interpolated across each of the individual voids (Fig. 1, V) to leave a new “infilled” DEM. This is then subtracted from the
original DEM (Fig. 1, VI), leaving a DEM of drumlin relief. Volume is calculated using drumlin area and height (Fig. 1, VII). 2.1. Study area With the primary requirement being a study area where a high spatial resolution DEM can be validated against geomorphological ground control, we selected sites from the western part of central Scotland (Fig. 2; Smith et al., 2006). The landforms studied are drumlins formed during the Last Glacial Maximum (LGM), when ice moved west to east across the region, and during the Younger Dryas (YD) glaciation, when ice flow was predominantly from north to south (Fig. 2). Thus, although the principles outlined in this paper are applicable to all glacier bedforms (flutes, megaflutes, drumlins, megadrumlins, streamlined hills and Rogen) the only examples studied are drumlins. Henceforward the term drumlin is substituted for bedform. The full study area covers ~750 km2 comprising a total of 2076 drumlins mapped in the field at 1:10,000 scale (Rose and Smith, 2008) of which 214 are from the area of YD glaciation which covers ~104 km2 (excluding Loch Lomond and the islands). For the purpose of demonstrating the methodology of this paper, two worked examples have been taken from the LGM glaciated region. Additionally a sample of 175 drumlins has been studied from an area that was in part last glaciated during the LGM and in part during the YD [note that the size of this sample differs from that (189) in Smith et al. (2006) because that reports all of the drumlins in the area, including parts of drumlins, whereas the sample used in this study only deals with complete drumlins]. All the drumlins used for these tests are identified from the NEXTMap DEM and the area used in the large sample is the same as that used for the comparison of field mapped results with NEXTMap results in Smith et al. (2006). The worked examples do not discuss the glaciodynamic implications of volumetric calculations; as it is necessary to understand the internal composition of the bedform before this process can be carried out, and this is an entirely different issue from that considered here. All figures depicted in map coordinates are projected to the British National Grid, with north–south indicated by vertical grid lines and north implicitly to the top of the figure. 2.2. Determination of bedform area Identification and representation of the drumlin outline requires the optimal visualization of the DEM. Smith and Clark (2005) suggest using images of local contrast stretch and gradient, supplemented with two orthogonal relief shaded images. The drumlin outline is then identified based upon breaks-of-slope and recorded as a vector polygon through on-screen digitizing. It is important to use visualization techniques that do not incorporate biases to the orientation of the drumlin (Smith and Wise, 2007), although Smith and Clark (2005) recommend the use of relief shaded images for the identification of topographically subtle drumlins. This process requires scientific judgement with volumetric calculations dependent upon the quality of digitization. Quantitative validation of the DEM mapping performed here is reported in Smith et al. (2006). 2.3. The Cookie Cutter
Fig. 1. Flow chart outlining the routine for the calculation of glacier bedform volume. This process requires visualization of a DEM (I and II) and the digitizing of drumlin outlines (III; a subjective operator process) which are then used to cookie cut voids from the original DEM (IV). The voids are then infilled (V) and the resultant data set is subtracted from the original DEM (VI) to produce estimates of relief. Volume can then be calculated (VII). IV to VII are performed on individual drumlins.
Cookie Cutter refers to the process whereby a drumlin is extracted from the original DEM (Fig. 1, IV). The term is derived from a cookie cutter used to cut shapes from pastry; here we cut drumlins from the DEMs. The cookie cutter cuts a single defined shape for each drumlin individually in order that each feature is treated separately. This is necessary in order to overcome problems associated with drumlin adjacency (where two drumlins are in contact), a normal feature of glacier bedform terrain (Fig. 3a), and is achieved by adopting a buffering process. If the drumlins were processed in a single operation, the adjacent rising slopes would influence the interpolation of the basal surface.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
211
Fig. 2. Location of the field mapping area (~ 750 km2) and comparative study area (~104 km2). The inset map (top left) shows generalized ice flow during the LGM (black arrows after Sissons, 1967) and YD (red outline after Rose, 1987). Map coordinates are projected to the British National Grid, with north–south indicated by vertical grid lines and north implicitly to the top of the figure.
The buffering process is only possible when the cookie cutter is performed as a separate operation for each drumlin (Fig. 3b). The process involves conversion of the area within the drumlin outline to a raster image comprised of a null value which is ignored in further calculations. This forms the void that is filled at the interpolation stage. The buffer used to separate adjacent drumlins can be user-specified, and in this instance we have defined a zone of 20 m around the drumlin outline. A conditional
function within the GIS is used to process each drumlin individually, on a cell-by-cell basis. Any cells in the buffered region are replaced by an elevation value derived from the level of the boundary of the DEM, whilst any cells within the drumlin are replaced with the null value. The output of this process is a drumlin void surrounded by a 20 m buffer with surface elevation identical to that of the drumlin boundary. This process is performed on every drumlin to produce DEMs with drumlin voids.
Fig. 3. Illustration of drumlin adjacency. (a) Digitized outlines of two drumlins that are adjacent. (b) The cookie cutter requires a buffer to be created around each drumlin (dashed line), a process not possible with a single operation as the buffered regions overlap. This stage therefore requires each drumlin to be processed individually.
212
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Fig. 4. Examples of drumlins on different topographic surfaces. (a) Diagrammatic representation of drumlin relief (relative elevation) on a planar horizontal surface depicting the difference in elevation from surrounding terrain, (b) drumlin relief on a planar tilted surface (e.g. hillside) and (c, d) depiction of drumlins formed on curving terrain configurations.
2.4. Determination of basal surface For the purpose of this research, the relief (relative elevation) of a drumlin is defined as the difference in elevation between an individual point on a drumlin surface and a plane interpolated across a drumlin outline which is “representative” of the base of the
drumlin (Figs. 1, V, 4a). This can be a simple relationship when the representative base is planar and level (Fig. 4a), but becomes much more complex when the representative base is inclined (Fig. 4b) or complex (Fig. 4c,d). The interpolated plane creates a basal surface for the drumlin void produced during the cookie cutting stage. It therefore prescribes the datum that is used for the calculation of relief
Fig. 5. Example of basal surface interpolation for a single drumlin. (a) relief shaded DEM, with digitized outline, (b) the blank void created after the removal of the drumlin, (c) the drumlin basal surface interpolated using IDW and (d) the drumlin basal surface interpolated using a tensioned thin-plate spline.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Fig. 6. Perspective view of a drumlin overlaid on a DEM grid. The interpolated base is shaded, with relief for each individual cell calculated from the base to the drumlin surface.
on a point-by-point basis (i.e. the elevation of the void will be different at every location). Relief is then calculated for all cells on the drumlin above this surface (Fig. 5) and these values are summed to cover the planform area and thus compute volume. Note that accurate calculations of bedform relief make assumptions about the terrain prior to the formation of the drumlin. This is a fundamental geomorphological problem that will be considered in a later paper; here we are concerned only with demonstrating and applying the cookie cutter technique. Although the above is simple in concept, calculation of relief can be difficult as bedforms can be found in complex terrain configurations, including sloping terrain and across ridges and depressions (Fig. 4b,c, d; Rose and Letzer, 1975, 1977; Mitchell, 1994). The interpolated surface needs to slope in order to take these variations into account and, therefore, is likely to be complex. A number of commonly used interpolation techniques were considered: inverse distance weighting (IDW), splines and kriging (Dubrule, 1984; Oliver and Webster, 1990). Kriging was discounted as, being a geostatistical technique, it ideally needs to be optimized for each application (Isaaks and Srivastava, 1990). Therefore both inverse distance weighting and splines were tested as to their suitability for this application. Inverse distance weighting is a simple and fast method in which an interpolated point is computed, based on the weighted average of its neighbours, where the weight of a neighbour is determined by its
213
distance from the point (Heine, 1986; Watson, 1992). The problem with using such a simple interpolator is that an artefact at the edge of the drumlin void (e.g. a hedge or man-made object) will, in effect, be “pulled” into the area of interpolation. This is illustrated in Fig. 5 which shows a drumlin crossed by a hedge (Fig. 5a). The void created by the cookie cutter is shown on Fig. 5b and the effect of the hedge is evident when the IDW interpolator is extended across this void. This hedge creates error because artifacts on a drumlin contribute to the definition of the basal surface (Fig. 5c). This problem can be resolved to some degree by making the area of interpolation as planar as possible with the use of a thin-plate spline under tension (Smith and Wessel, 1990; Fig. 5d). The thin-plate spline under tension can be described, through analogy, as a membrane that is pulled taught to the edges of the drumlin void; the force can be adjusted by a tension parameter. In this case, artifacts on the edges of the drumlin are smoothed, leaving the basal surface with increased planarity (Fig. 5d). The tension parameter is defined between 0 and 1, where 0 has no tension and means the interpolated surface conforms to the minimum curvature and 1 has infinite tension and maximum curvature. Infinite tension is not physically meaningful (Smith and Wessel, 1990), but implies that any further increase of tension cannot result in a smoother surface while remaining true to the original data points. Infinite tension was chosen to minimize the effects of artifacts whilst conforming to the surrounding terrain. As can be seen in Fig. 5 the IDW method produces a relatively poor result as the hedges can be seen to “smear” across the drumlin area. The thin-plate spline performs well, defining a surface around the base of the drumlin. The thin-plate spline was therefore selected for the interpolation of surfaces across all drumlin voids. 2.5. Determination of bedform relief and volume After the interpolation of the drumlin basal surface, relief can be calculated by subtracting the interpolated basal values (excluding the buffer) from the original DEM comprised of drumlin elevations (Fig. 1, VI). Each cell in the output dataset contains an elevation value representing the height of the drumlin above its basal surface (Fig. 6). When this value is multiplied by the area of each cell (25 m2 in this instance), an estimate of volume of the cell column is produced. The summation of the volumes of all cell columns that comprise a drumlin produces a final total volume for the individual bedform.
Fig. 7. Flow diagram showing the Cookie Cutter computer processing stages (circles represent input parameters). (a) Visualisation of an input DEM is used as a data source to digitize drumlin outlines, (b) a project directory is specified prior to running the “Cookie Setup” script and (c) separated drumlin outlines and the original DEM form the primary inputs to the “Cookie Cutter” script. Four input parameters are required and, once complete, the script outputs drumlin relief and summary statistics.
214
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
3. Worked examples on small samples The previous sections have outlined the principle steps required to estimate drumlin volume. This section outlines the procedure in detail and then provides two worked examples illustrating how this method can be applied. 3.1. Procedure (Fig. 7) The procedure for carrying out this process is as follows (Fig. 7): i) Digitization of drumlin outlines is performed manually within ESRI ArcMap and saved as a polygon shapefile (Fig. 7a) ii) The cookie-cutting process is automated using two custom ArcGIS Python scripts which are loaded and run from ArcToolBox. The Python scripts are freely available for download as ancillary materials accompanying this paper. The “Cookie Setup” script (Fig. 7b) separates each digitized drumlin outline into a separate file, as well as establishing a series of directories for computer processing. Input in to the script is a shapefile of digitized drumlin outlines, with the user required to specify the project directory. The “Cookie Cutter” script (Fig. 7c) buffers the original drumlin outline in order to produce the drumlin void independent of
interference from other drumlins. The spline interpolator within ArcGIS fills the void before the production of the final output. A raster file represents the relief for each drumlin, supplemented with a database collating summary statistics for each drumlin. Input into the script is a DEM, with drumlin outlines automatically selected from the output of the “Cookie Setup” script. The user is required to specify the project folder, cell size (it is recommended to use the cell size of the input DEM), tension parameter for the spline interpolator and distance that each drumlin will be buffered (at least four cells). 3.2. Bowridge (Fig. 8) Fig. 8a depicts a well defined megadrumlin with three superimposed drumlins that are in the order of hundreds of metres long and tens of metres high. The location is shown on Fig. 2; the hill is known as Bowridge and is centred on 1 km grid square NS 7880. Fig. 8b shows a standard greyscale visualization of the DEM and demonstrates the ineffectiveness of this type of image. Fig. 8c shows the most effective visualization of the drumlins after rotating light source directions on the NEXTMap image and the boundaries of these drumlins are outlined. The drumlin outlines were buffered and then used to cookie cut the DEM to create a series of drumlin voids (Fig. 8d). A tensioned
Fig. 8. Worked example demonstrating the steps required to calculate drumlin volume —(a) the mega drumlin located at Bowridge (NS 7880) (see Fig. 2 for location), (b) raw DEM data, (c) relief shaded visualization of terrain, with mapped drumlin outlines, (d) DEM voids, (e) interpolation of drumlin basal surfaces and (f) relief shaded visualization of drumlin volumes (1.51 × 106 m3). Note the “stepping” in (e), a result of artifacts at the edges being interpolated across the basal surface.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
215
Fig. 9. Worked example demonstrating the steps required to calculate drumlin volume — the drumlin located near Milton of Campsie (NS 6677) (see Fig. 2 for location). (a) original field mapping, (b) raw DEM data, (c) relief shaded visualization of terrain, with mapped drumlin outlines, (d) DEM voids, (e) interpolation of drumlin basal surfaces and (f) relief shaded visualization of drumlin volumes (2.07 × 106 m3). Note the “stepping” in (e), a result of artifacts at the edges being interpolated across the basal surface. Note also the difference between the interpretations of drumlin form and extent derived from field mapping and analysis of NEXTMap showing the scope of the two methods. This difference does not affect the calculations but illustrates the need for ground truthing whenever possible (Rose and Smith, 2008).
spline was then used to interpolate the basal surface of the drumlins (Fig. 8e; note the effect of edge artifacts on the smoothness of the basal surfaces). This basal surface is subtracted from the original drumlins to enable a calculation of drumlin volumes (Fig. 8f; note that this final image has been relief shaded to help visualization).
smoothness of the basal surfaces) and the final drumlin volumes were calculated (Fig. 9f).
3.3. Milton of Campsie (Fig. 9)
In addition to developing a technique for the estimation of drumlin volume, or the volume of any other landform, we have applied the method to a small area from the western part of central Scotland that includes drumlins formed during both the LGM and the YD glaciations. This area covers ˜ 104 km2 and contains 175 drumlins mapped from the NEXTMap DEM, ranging from 88 to 844 m in length (see Smith et al., 2006 for further details). The drumlins formed during the LGM and the YD glaciations are identified on Fig. 10. Table 1 presents drumlin relief and drumlin volume results for the entire study area and, separately, the sub-samples for drumlins formed during the LGM and YD. The minimum drumlin relief value of
The second worked example was selected as the drumlins are larger and clearly arranged adjacent to one another, with shared outline boundaries. The field mapping (Fig. 9a) and analysis of the DEM (Fig. 9c) depict the same general arrangement of drumlins and, as at Bowridge, the greyscale image is ineffective (Fig. 9b). The drumlin outlines are depicted on the NEXTMap image following analysis of the NEXTMap visualizations (Fig. 9c). Fig. 9d shows the cookie cut drumlin outlines; basal surfaces were subsequently interpolated (Fig. 9e; again note the effect of edge artifacts on the
4. Comparing the volume of drumlins formed during LGM and YD glaciations
216
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Fig. 10. Comparative study area depicting drumlins mapped from NEXTMap Britain DEM data. The dotted line demarcates the boundary between LGM and YD drumlins. Map coordinates are projected to the British National Grid, with north–south indicated by vertical grid lines and north implicitly to the top of the figure.
−16 m (and minimum drumlin volume of −1.15 × 103 m3) indicates that the results must contain error as it is not possible to have negative relief values. Investigation of the sources of these errors indicated eight negative values that resulted from the creation of “sliver polygons” (Longley et al., 2005) (Fig. 11a) or “clutter” (such as trees or hedges) at the drumlin edges creating the effect of negative relief when the basal surface is interpolated (Fig. 11b). Sliver polygons result from incorrect digitization and can simply be deleted. Drumlins with negative relief values resulting from clutter were considered erroneous and excluded from further analysis. The summary statistics (Table 1) indicate that relative to the LGM sub-sample of drumlins the volumes for the YD sub-sample give a smaller mean size with a larger
range and larger standard deviation. A t-test comparing drumlin volumes for the two sub-samples, with equal variances not assumed, was statistically significant (0.012) indicating that the datasets are different. This fact is reported here as an outcome from this exercise,
Table 1 Volume calculations for the study area (figures in bold are corrected after the removal of negative relief values), incorporating the two sub-samples for drumlins formed during the YD and LGM. Measure
Study area (n = 175)
Relief min (m) Relief max (m) Relief mean (m) Relief standard deviation (m) Volume min (106 m3) Volume max (106 m3) Volume mean (106 m3) Volume standard deviation (m) Total volume (106 m3)
− 16.3 65 3.7 1.8 − 0.11 2.75 0.12 0.24 22.75
1.7 65 9.3 7.0 0.0016 2.75 0.13 0.24 23.00
Younger Dryas (n = 104)
Last Glacial Maximum (n = 71)
1.65 65 8.1 7.2
2.41 31 11.1 6.2
0.0016 0.70 0.12 0.29 12.31
0.0031 2.75 0.15 0.16 10.69
Fig. 11. Two common errors associated with automated calculation of drumlin volumes. (a) Mis-digitization can create erroneous “sliver” polygons. (b) The presence of edge artifacts can lead to the interpolation of basal surfaces not related to the drumlin topography. The dotted line shows the drumlin basal surface, whilst the dashed line shows the actual interpolated surface as a result of the tree increasingly the elevation at the edge of the drumlin.
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
but the results may identify a critical difference between drumlin samples formed by glaciers with different dimensions, and may confirm the critical importance of this property in controlling bedform dimensions, as proposed in Rose (1989). 5. Discussion The Cookie Cutter is a straightforward approach to calculate the volume of drumlins, other glacier bedforms and well defined positiverelief landforms in general. However, the derived values presented here are only a first approximation due to several potential sources of error that introduce uncertainty into the results. In addition to the errors inherent in the method of collecting the data (quality of the DEM and quality of the geomorphological interpretation; Smith et al., 2006), further sources of uncertainty are outlined below: a) Quality of digitization of interpreted data — this can mean that terrain may be inappropriately included or excluded from the calculation. For example, pronounced furrows can occur around the stoss end of drumlins and along their sides, erroneously incorporating volume in to the calculations. b) DEM vertical accuracy — the accuracy of the measurement of surface elevation will affect the calculation of drumlin volume. The vertical accuracy of NEXTMap is ~0.5–1.0 m. c) Surface Clutter — the NEXTMap DEM used in this study quantifies visible surface elevation, including buildings and vegetation. These features are not relevant to the volume of drumlins and it is desirable that they are removed through an automated procedure (e.g. Wang et al., 2001; Sithole and Vosselman, 2004), although this may have the detrimental effect of removing subtle topographic features, including, for instance small drumlins. Removal through interactive editing of the DEM is also possible, although this can be a time consuming task. Buildings and forestry cover only small parts of the study sites, but hedges are an important element of the landscape and caused problems throughout. Of particular importance is the presence of clutter at the edges of drumlins which can substantially affect drumlin volume estimates and, in some instances, lead to calculations of negative volume. Calculations presented for the study area above are therefore considered to be underestimates, although it should be noted that the presence of clutter on drumlins will cause overestimation. d) Interpolation — the drumlin basal elevation values are estimated through the use of an interpolator and therefore a degree of uncertainty will be introduced as these represent estimates of an unknown, prior surface. In addition to the methodological problems listed above, it should be noted that in-situ landform degradation will have occurred. The drumlins examined in this study were exposed to subaerial process as the ice retreated some c. 15 ka BP and c. 11.8 ka BP following the LGM and YD glaciations respectively. Mass wastage processes, therefore, will have degraded the glaciogenic landforms to varying degrees (Rose, 1991). The combination of the above factors will contribute to error in the calculation of drumlin volume. Potential mitigation measures include the following: i) the use of LiDAR data with a higher spatial resolution, higher vertical accuracy and ability to penetrate vegetation canopies and b) a prescriptive procedure for digitization which will improve consistent mapping. 6. Conclusions This paper presents a technique for the calculation of material volumes of drumlins from sites in western central Scotland. The method can be applied to other glacial bedforms and indeed to any landform. Conclusions from this work are:
217
• Digitized drumlin outlines are used to extract (“cookie cut”) landforms from an underlying DEM leaving empty voids. • A void is “infilled” through the application of a tensioned spline interpolator thereby estimating the basal surface. • Drumlin relief is calculated by subtracting the basal surface from the original DEM and then converted to volume by multiplying by the planform area. • Landform adjacency requires drumlins to be processed individually (to avoid the creation of adjacent voids). The drumlin outlines are therefore buffered (to include surrounding terrain) prior to being “cookie cut”. • A protocol is established, using the above methodology, for the calculation of material volumes for similarly wholly concave or wholly convex landforms. • The results provide a first-order approximation that facilitates the rapid calculation of drumlin volumes and, in the example given in this paper, the differentiation of the volume of drumlins formed by an ice sheet during the LGM, and drumlins formed at the margin of the piedmont lobe of the Loch Lomond outlet glacier during the YD. • Volume calculations are subject to error from the source DEM, digitization procedure and presence of surface clutter, in addition to geomorphological problems inherent in the sample studied. • This technique can potentially be used wherever landform volume estimates are required, for example, calculations of sand dune volumes or sediment input into depositional environments.
Acknowledgements We thank Steve Booth and James O'Brien for their assistance with this work. Reviews from Ian Evans and John Shaw helped improve this manuscript. The NEXTMap DEM data was kindly supplied by the British Geological Survey. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.geomorph.2009.01.006. References Boulton, G.S., 1987. A theory of drumlin formation by sub-glacial sediment deformation. In: Menzies, J., Rose, J. (Eds.), Drumlin Symposium. Balkema, Rotterdam, pp. 25–80. Boulton, G.S., Clark, C.D., 1990. A highly mobile Laurentide Ice Sheet revealed by satellite images of glacial lineations. Nature 346, 813–817. Charlesworth, J.K., 1924. The glacial geology of the north-west of Ireland. Proceedings of the Royal Irish Academy, vol 36B, pp. 174–314. Clark, C.D., 1997. Reconstructing the evolutionary dynamics of palaeo-ice sheets using multi-temporal evidence, remote sensing and GIS. Quaternary Science Reviews 16, 1067–1092. Close, M.H., 1867. Notes on the general glaciation of Ireland. Journal of the Royal Geographical Society of London 1, 207–242. Dubrule, O., 1984. Comparing splines and kriging. Computers and Geosciences 10, 327–338. Evans, I.S., 1987. A new approach to drumlin morphometry. In: Menzies, J., Rose, J. (Eds.), Drumlin Symposium. Balkema, Rotterdam, pp. 119–130. Heine, G.W., 1986. A controlled study of some two-dimensional interpolation methods. Computer Oriented Geological Society Computer Contributions 2, 60–72. Isaaks, E.H., Srivastava, R.M., 1990. An Introduction to Applied Geostatistics. Oxford University Press, New York. Intermap, 2005. Product Handbook [online]. Available from: http://www.intermap. com/images/handbook/producthandbook.pdf [Accessed: 21st October 2005]. Kleman, J., 1994. Preservation of landforms under ice sheets and ice caps. Geomorphology 9, 19–32. Kleman, J., Hätterstrand, C., Borgström, I., Stroeven, A., 1997. Fennoscandian palaeoglaciology reconstructed using a glacial geological inversion model. Journal of Glaciology 43, 283–299. Longley, P.A., Goodchild, M.F., Maguire, D.J., Rhind, D.W., 2005. Geographical Information Systems and Science. Wiley, Chichester. Mitchell, W.A., 1994. Drumlins in ice sheet reconstructions, with reference to the western Pennines, northern England. Sedimentary Geology 91, 313–331. Oliver, M.A., Webster, R., 1990. Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Information Systems 4, 313–332.
218
M.J. Smith et al. / Geomorphology 108 (2009) 209–218
Rose, J., 1987. Drumlins as part of a glacier bedform continuum. In: Menzies, J., Rose, J. (Eds.), Drumlin Symposium. Balkema, Rotterdam, pp. 103–116. Rose, J., 1989. Glacial stress patterns and sediment transfer associated with the formation of superimposed flutes. Sedimentary Geology 62, 151–176. Rose, J., 1991. Subaerial modification of glacier bedforms immediately following ice wastage. Norsk Geografisk Tidsskrift 45, 143–153. Rose, J., Letzer, J., 1975. Drumlin measurements: a test of the reliability of data derived from 1:25,000 scale topographic maps. Geological Magazine 112, 361–371. Rose, J., Letzer, J., 1977. Superimposed drumlins. Journal of Glaciology 18, 471–480. Rose, J., Smith, M.J., 2008. Glacial geomorphological maps of the Glasgow region, western central Scotland. Journal of Maps v2008, 399–416. Sissons, J.B., 1967. The Evolution of Scotland's Scenery. Oliver and Boyd, Edinburgh. Sithole, G., Vosselman, G., 2004. Experimental comparison of filter algorithms for bareEarth extraction from airborne laser scanning point clouds. ISPRS Journal of Photogrammetry & Remote Sensing 59, 85–101. Smalley, I.J., Unwin, D.J., 1968. The formation and shape of drumlins and their distribution and orientation in drumlin fields. Journal of Glaciology 7, 377–390.
Smith, W.H.F., Wessel, P., 1990. Gridding with continuous curvature splines in tension. Geophysics 55, 293–305. Smith, M.J., Clark, C.D., 2005. Methods for the visualisation of digital elevation models for landform mapping. Earth Surface Processes and Landforms 30, 885–900. Smith, M.J., Wise, S.M., 2007. Mapping glacial lineaments from satellite imagery: an assessment of the problems and development of best procedure. International Journal of Applied Earth Observation and Geoinformation 9, 65–78. Smith, M.J., Rose, J., Booth, S., 2006. Geomorphological mapping of glacial landforms from remotely sensed data: an evaluation of the principal data sources and an assessment of their quality. Geomorphology 76, 148–165. Wang, Y., Mercer, B., Tao, V.C., Sharma, J., Crawford, S., 2001. Automatic generation of bald earth digital elevation models from digital surface models created using airborne IFSAR. Proceedings of 2001 ASPRS Annual Conference, St. Louis, Missouri, USA. Watson, D., 1992. Contouring: A Guide to the Analysis and Display of Spatial Data. Pergamon Press.