Last modified on Jun 07, 2024 14:12

Contributors: Jacqueline Bannwart (University of Zurich), Inés Dussailan (University of Zurich), Frank Paul (University of Zurich), Michael Zemp (University of Zurich)

Issued by: UZH / Frank Paul

Date: 02/02/2024

Ref: C3S2_312a_Lot4.WP2-FDDP-GL-v2_202312_A_ATBD-v5_i1.0

Official reference number service contract: 2021/C3S2_312a_Lot4_EODC/SC1

History of modifications

List of datasets covered by this document 

Related documents 

Acronyms

General definitions

Brokered data set: A dataset that is made available in the CDS but freely available (under given license conditions) from external sources. In the case of the glacier distribution service, the Randolph Glacier Inventory (RGI) is brokered for the CDS from https://glims.org/RGI under a CC-BY 4.0 license.

Debris-cover: Debris on a glacier is usually composed of unsorted rock fragments with highly variable grain size (from mm to several m). These might cover the ice in lines of variable width separating ice with origin in different accumulation regions of a glacier (so called medial moraines) up to a complete coverage of the ablation region. Automated mapping of glacier ice is only possible when the debris is not covering the ice completely when compared to the image pixel size.

Glacier complex: A contiguous ice mass that is the result of the binary (yes/no) glacier classification after conversion from raster to vector format. Usually, the glacier complexes are divided into individual glaciers by digital intersection with a vector layer containing hydrologic divides derived from watershed analysis of a digital elevation model (DEM).

Scope of the document

This document is the draft Algorithm Theoretical Basis Document (dATBD) for the glacier area product from the Copernicus glacier distribution service. It describes the algorithms used to create the product (glacier area), along with the satellite and auxiliary datasets required for the processing, as well as the format of the output datasets. All sections describe the data selection or processing approach as applied in the Copernicus Climate Change Service (C3S) to generate the products for the Global Land Ice Measurements from Space (GLIMS) database. For specific processing steps, we have formulated our approach as recommendations for other analysts, being aware that also other approaches are used. This document is a condensed version of the ATBD for the Glaciers_cci project, which gives full details and is available online [RD1].

The worldwide glacier area/outline dataset that is brokered as a Climate Data Record (CDR) to the Climate Data Store (CDS) is the Randolph Glacier Inventory (RGI). The RGI is a vector dataset extracted from the GLIMS database presenting for each glacier one outline (GLIMS is multi-temporal). The GLIMS database consists of several local to regional-scale glacier outline datasets provided by the scientific community (including C3S and C3S2) and has thus variable quality. For the RGI the outlines with the best quality have been selected. In this document we additionally explain how the RGI vector dataset is converted to a raster dataset.

Executive summary

Glacier outlines can be derived from a variety of optical sensors using semi-automated methods when a shortwave-infrared (SWIR) band is available, or from manual digitisation if not. Satellites, such as Landsat 5, 7 or 8 and Sentinel-2, provide a SWIR band and offer large spatial coverage combined with a sufficient to very good spatial resolution (30 to 10 m) to map this rapidly changing target precisely, and have thus been widely applied. The most important first step is selection of appropriate images from the sensors of choice. The scenes have to be cloud free over glaciers and seasonal snow should not hide the glacier perimeter. In many glacierized mountain regions this is difficult to achieve with a single scene and requires the creation of glacier outlines from a mosaic of several scenes, often acquired over several years.  

Whereas the method to map clean glacier ice is very simple (a red/SWIR band ratio with a threshold), intense manual editing is often required during post-processing, for example to add missed regions with debris cover or under clouds and in shadow, or to remove misclassification such as turbid water. The editing is done in the vector domain and usually performed against the satellite images (using different and contrast enhanced spectral band combinations) used to create the raw outlines. Hence, each glacier outline is already validated against a reference dataset, i.e. systematic errors are removed to obtain a dataset that meets the accuracy requirements of the Global Climate Observing System (GCOS). However, uncertainties (or random errors) due to a variable image interpretation and digitizing consistency remain.  

Once outlines have been corrected a Digital Elevation Model (DEM) of good quality is required for the region of interest to derive drainage divides and topographic attributes for each glacier. Both can be calculated from automated scripts, but the divides usually also require intense manual editing. Once the divides are good, they are digitally intersected with the corrected glacier extents to obtain individual glaciers from the mapped glacier complexes. The resulting vector layer is combined with the DEM (and products derived from it) to derive topographic attributes for each glacier.

For the raster dataset we determine the percentage of glacier cover per grid cell by dividing the area of the glacier cover in this cell by the area of each grid cell. For this purpose, the boundaries (ice divides) between the glaciers are dissolved and the entire global glacier area is handled as a multipart polygon with the same ID. This global vector layer is digitally intersected with another vector layer in which each polygon has a defined size of 1° by 1°, i.e. for the entire Earth it consists of 360 by 180 cells. 

1. Instruments

1.1. Optical sensors

Glacier outlines can be derived from imaging sensors operated on satellites, aircraft or drones. In order to map (debris free) glacier ice automatically (e.g. with the band ratio method described in Section 3), it is required that at least one sensor is located in the short-wave infrared (SWIR) part of the spectrum (around 1.5 m). If only visible and near-infrared (VNIR) bands are available, outlines have to be digitised manually. Manual editing is also required to correct misclassification in automatically classified images. The manual digitizing requires sensors with an appropriate spatial resolution. Previous studies have shown that the 30 m resolution of the sensors on-board Landsat is sufficient to identify and map glaciers down to a size of 0.01 km². Sensors with a higher spatial resolution, such as the 15 m Advanced Spaceborne Thermal Emission and reflection Radiometer (ASTER) (but 30 m SWIR) or the 10 m Sentinel-2 (20 m SWIR), are even better, as the interpretation of the glacier boundary for manual correction is much easier (Paul et al. 2016). For these freely available sensors, the spectral ranges and band numbers are listed in Table 1.

At even higher spatial resolutions (a few m to 10 cm) sensors can still be used to map glaciers, but some disadvantages occur: 
(1) they are usually commercially distributed and are very expensive for larger areas,
(2) they only cover comparatively small regions so that a limited number of glaciers can be mapped,
(3) they do not have a SWIR band so that glaciers can only be delineated manually.

Manual delineation strongly increases the workload compared to the automated mapping with a SWIR band, particularly in regions dominated by debris-free glaciers (that can be mapped automatically). In consequence, such very high-resolution sensors are mostly used for more detailed studies including validation (Section 2) rather than for large-scale glacier mapping. However, some high-quality inventories have been created from such sensors (Fischer et al. 2014, Linsbauer et al. 2021).

Common problems such as different spectral ranges (Table 1), sensor sensitivities and calibration coefficients (incl. ageing) of the non-commercial satellites (Landsat, Terra ASTER, Sentinel) used here do not impact on glacier mapping results, as the method applied (band ratio with a threshold, see Section 3) largely compensates for these and other (e.g. atmospheric conditions) differences. Hence, common issues related to sensor stability do not apply and the full time series of sensor data can be used to generate the climate data record from the satellite database (1984 to present).

Table 1: Spectral ranges of individual satellite image bands for several optical sensors (from Paul et al. 2016). Colours indicate spatial resolution (black: 30 m, red: 20 m, blue: 15 m, green: 10 m). Sensor names are explained in the Acronyms list, AST is ASTER.

A second point, related to stability, is the required orthorectification of the satellite images. These days this is performed automatically by the distributing space agencies. The DEMs used for it so far differed and lacked the required spatial resolution and accuracy required for correction of the panoramic distortion in steep high-mountain terrain (Toutin 2004). In consequence, outlines derived from satellite images that have been orthorectified with different DEMs can have a spatial mismatch and cannot be used directly for time-series analysis and change assessment (Kääb et al. 2016). Accordingly, scalar values of glacier area have to be derived for each individual glacier and subtracted from the value obtained for a different date to determine area changes. This has also to be considered when outlines are compared that have been created by different analysts. In the case where the outlines are (partly) manually corrected, differences in interpretation of the glacier extent can be easily larger than real area changes over a decadal period (Paul et al. 2015).  

With the new Copernicus DEM (see link in Table 2) available and usage of this dataset by USGS and ESA for orthorectification, the geolocation issue should now be solved. However, detailed tests with historic and current datasets have to be performed once the reprocessing is done.  

Detailed assessments of outline uncertainty (round-robin) have shown that the digitizing uncertainty by the same analysts is in the same order of magnitude (±½ pixel) as the results obtained by automatically derived glacier outlines for debris-free ice (Paul et al. 2013). Overall, the sensor used and method applied for glacier mapping does not play a large role for product quality. Instead, quality is largely governed by decisions of the analyst regarding debris, clouds, snow, shadow (interpretation uncertainty) and image conditions (e.g. regarding snow and cloud cover).

1.2. Microwave sensors

Satellite Synthetic Aperture Radar (SAR) sensors are not directly used in glacier mapping but can provide supportive information. Their most valuable application is in helping to interpret the extent of debris-covered glaciers from coherence images (e.g. Atwood et al. 2010, Frey et al. 2012). Combined with optical images and glacier outlines, the maps of low coherence between SAR image pairs can be interpreted as regions of active movement and thus belonging to the glacier. The method has been successfully tested with Advanced Land Observing Satellite (ALOS)-1 Phased Array Type L-Band Synthetic Aperture Radar (PALSAR) pairs (48 days repeat) but also works well with recent image pairs from TanDEM-X (11 days), Sentinel-1 (12 days), and ALOS-2 PALSAR-2 (14 days).

2. Input and auxiliary data 

2.1. Input data

The decision on the scenes to be selected for processing is with the analyst. In general, a compromise has to be made between clouds, remaining snow cover and shadow extent. Experience has shown that product quality is best when only scenes are used with glacier perimeters not covered by seasonal snow. This means that in some regions appropriate data are only available late in the autumn when shadows can already be too deep to identify all glaciers. In this case scenes taken earlier in the year might help but often clouds, snow conditions and shadow can compromise scene selection for several years. The high revisit frequency from Sentinel-2A/B and Landsat 8/9 of 5 and 8 days at the equator (more often towards the poles), respectively, helps to overcome the acquisition date issue in some cases. Combining scenes from more than one year might also be a solution.

The satellite data covering large regions with appropriate spatial resolution and being freely available are limited to Landsat (sensors Thematic Mapper (TM), Enhanced Thematic Mapper plus (ETM+) and Operational Land Imager (OLI)) and Sentinel-2 Multi Spectral Imager (MSI). Additionally, Terra ASTER scenes might be used but they cover only 1/9 of a Landsat scene and processing workload is thus much higher. is providing an overview on spectral band ranges / numbers and spatial resolution of these sensors. Orthorectified raw data are available from several sources, for example http://earthexplorer.usgs.gov for Landsat data and the Sentinel science hub for all Sentinel data (http://scihub.copernicus.eu). They are provided in geoTIFF or JPEG 2000 (JP2) formats that can be handled by a wide range of Geographic Information Systems (GIS) and digital image processing software. As a note of caution, the databases are constantly extended (e.g. Wulder et al. 2016) and sometimes datasets are reprocessed. It might, thus, be possible that at some point in time, scenes are available that have not been in the archive before. The archive should thus always be carefully checked before processing in a certain region is performed.

2.2. Auxiliary data

The three main auxiliary datasets used when creating glacier outlines and inventories are:

 (1) DEMs to create drainage divides and derive topographic information for each glacier entity,
 (2) high-resolution datasets for quality control during manual correction (available from Google Earth or web map services), and
 (3) SAR-coherence images in case of debris-covered glaciers.

The DEMs typically used for (1) are listed in Table 2 along with their resolution, sensor, coverage and data source. 

Table 2: Overview of DEMs typically used when creating glacier inventories. The first and the last DEMs in the list have been generated from SAR sensors, all others are from optical sensors (see acronym list). All URL resources on this table last viewed on 02^nd February 2024.

¹⁾ WORLDVIEW-1, 2, 3, (GEOEYE-1)
²⁾ GEOEYE-1, QUICKBIRD-2, WORLDVIEW-1, 2, 3

Deriving (1) hydrologic drainage divides from a DEM requires the use of a product that is most accurate at high elevations (e.g. along terrain ridges), whereas for topographic parameters a good temporal match with the satellite scene is more important in order to correctly represent the elevation of the glacier surface. Ridges in Shuttle Radar Topography Mission (SRTM) are sometimes wrong or missing due to data voids, so the ASTER Global DEM (GDEM) or the ALOS World 3D 30 (AW3D30) DEM might be a better choice worth checking. However, both might suffer from artefacts and the position of the related divides has, thus, a higher uncertainty. We do thus check several DEMs for the region of interest and select the best one on a case-by-case basis. Apart from data voids, the ArcticDEM is now often the best choice in Arctic regions, closely followed by the somewhat coarser resolution Copernicus DEM that has been derived from SAR data rather than optical sensors. For topographic parameters and scenes acquired around 2000 the SRTM is the best, while after 2005 the GDEM becomes superior. Despite the differences in the methodologies underlying their creation as well as their temporal differences, all DEMs are equally suitable to derive topographic parameters over glaciers (Frey and Paul 2012). For outlines derived from Landsat 8/9 and Sentinel-2, the TanDEM-X DEM would be the best choice, but it is currently 'only' freely available at 30 m resolution (10 m in Europe). A detailed inter-comparison (ice divides and topographic information) of DEMs available in Arctic regions has been performed by the Glaciers_cci project revealing that all DEMs are suitable despite their (mostly local) artefacts.

The high-resolution datasets available in Google Earth or from web-map services for (2) sometimes suffer from adverse snow or cloud conditions so that they cannot always be used to guide the interpretation. Moreover, the image acquisition date might be unknown and/or very different from the image used to create the outlines. This has to be considered when interpreting the data. Other high-resolution reference datasets distributed commercially are expensive, so quality control with these datasets is seldom performed and other (internal) accuracy measures are applied instead. SAR sensors that can be employed to (3) compute coherence images of high contrast on debris-covered glaciers) include ALOS-1/2 PALSAR (2007-2011/2014-), TerraSAR-X (2003-) and Sentinel-1 (2014-). Also, the variable and likely different acquisition dates have to be considered when terminus positions are corrected with such images. 

A key problem with all of the above datasets is the possible spatial mismatch of orthorectified satellite scenes with the selected DEM, if a different DEM has been used for orthorectification. The related shifts in geolocation in steep high-mountain terrain might result in drainage divides (which follow the mountain ridges) cutting through glacier boundaries after digital intersection. This creates numerous very small polygons inside the original glacier boundary, named sliver polygons. In principle, all of them have to be manually corrected, but often they are later on just removed with a size threshold. This reduces the size of the original glacier by a few percent. The problem is difficult to overcome, in particular when working at large regional scales. With the new Copernicus DEM being available for both orthorectification and calculation of ice divides this problem might have been solved, at least after the related scenes have been reprocessed.

3. Algorithms

As a general note, we here describe how data products are created by C3S. As similar methods are applied by the community, we also include specific recommendations. We focus here on the most important aspects of the method; further details can be found in the Glaciers_cci ATBD [RD1].

3.1. Scientific background

Glaciers result from the metamorphosis and compression of snow. Accordingly, the spectral properties of (clean) glacier ice is very similar to snow. Additional to snow and ice, the glacier surface might be covered with other material exhibiting different spectral properties, e.g. rock (sizes from mm to m), soot, dust, other pollution and liquid water (Figure 1). Depending on the spatial density of these materials, the spectral information in a satellite image pixel (in the 10-30 m range) is in most cases a mixed signal with the respective deviations from a pure (laboratory) signal. In Figure 2a and b a comparison of the spectral reflectance for snow of varying grain size from theoretical considerations (Dozier 1989) with field-based measurements from Qunzhu et al. (1983) is shown. Besides the high reflectance of snow in the visible part (VIS) of the spectrum (independent of grain size), the strong reflectance drop in the near infrared (NIR) and a much stronger dependence on grain size can be seen (Figure 2a). In the SWIR the reflectance is very low and still depending on grain size. By dividing the high digital numbers (DNs) over ice and snow in the red band by the very low ones in the SWIR, a strongly contrast enhanced image is created that is converted with a threshold to a first raw glacier map (see Section 3.3). The principle steps of the pre, main and post processing are shown in Figure 3 and described in the following.

Figure 1: Oberaarglacier in Switzerland. The picture illustrates the differing reflectance of ice, debris-cover, snow, rock, vegetation (V) and turbid water in the visible part of the spectrum. A spectral discrimination of the debris cover on the glacier and from the lateral moraine (in the lower left of the image) is not possible.

Figure 2: Modelled spectral reflectance curves of snow with three different grain sizes and position of TM spectral bands (left). Spectral reflectance curves of snow, firn, ice and dirty glacier ice according to field measurements (right). The data for the left figure are taken from the ASTER spectral library (JPL 2002), the right figure is adapted from Hall et al. (1988).

Figure 3: Principle data processing workflow, copied from the ATBD of Glaciers_cci [RD1]. The pre-processing step ‘Scene selection’ is described in Section 2.1.

3.2. Pre-processing

Glacier outlines created within C3S2 are derived from the well-established band ratio method (Section 3.3). This requires downloading a VNIR and the SWIR band from the satellite under consideration (see Table 1 for band numbers of the typically used sensors). As the red band gives often better results than the near infrared (NIR) band, we use the former. For the later manual correction (post-processing stage) of ice and snow in shadow, debris cover, water and clouds, RGB composites with the red, green, blue (natural colours) and NIR, red, green (false colour infrared) bands should be created. As the algorithm works best with the raw digital numbers (DN), raw data rather than already converted reflectance products are obtained from the respective sensor. Raw data comes in geoTIFF or JPEG 2000 (JP2) formats that can be handled by a wide range of Geographic Information Systems (GIS) and digital image processing software. We code at least sensor, path/row (or tile identifier) and date in the file name (e.g. l8_195027_150830 or s2_32tms_150829) and append further processing details with letters and numbers.

The main steps in pre-processing of satellite scenes to derive glacier outlines are:

I. Selection of the region, analysis of the scenes in the archive and selection of the best ones. Depending on snow and cloud conditions, scenes from different dates might need to be considered.
II. Download of all bands and metadata for all scenes, and possibly format conversion.
III. Creation of various RGB composites, possibly in an image stack.
IV. Careful analysis of image conditions, e.g. clouds requiring correction during post-processing.

3.3. Main Processing

The main processing steps are shown in Figure 4 and include:

(a) The calculation of a band ratio (red/SWIR or NIR/SWIR) using the raw data (Figure 4c)
(b) Application of a threshold t₁ to separate ice and snow from other terrain
(c) Possible application of a second threshold t₂ on the blue band (or green if blue not available) to improve classification in shadow or exclude dark ocean water
(d) Possible application of a spatial filter (e.g. 3 by 3 majority) to reduce noise
(e) In the final binary yes/no glacier map (Figure 4d), the 'other' pixels are set to 'no data' before the raster file is converted to vector outlines. These are then further corrected in the post processing step as described below.

Step (a) is the main classification step to map glaciers. In pseudo code it can be written as:

IF (red/SWIR) > th1 [AND blue > th2] THEN ‘glacier’ ELSE ‘other’

A recent study found that the Landsat 8 OLI panchromatic band can be used instead of the red band for the band ratio (Paul et al. 2016). When the SWIR band is resampled (using bilinear interpolation) from 30 to 15 m (for Landsat, from 20 to 10 m for Sentinel-2) before processing, the spatial resolution of the outline product can be doubled. Steps (b) and (c) are iterative. The glacier maps resulting from different thresholds are overlaid with the natural colour RGB image to check if glaciers in shadow are accurately mapped. The threshold value th₁ (or th₂) needs to be lowered/increased if too little/too much ice in shadow is mapped. The method works reliably because snow and ice have a very different (much higher) reflectance in the VNIR (Figure 4a) than in the SWIR (Figure 4b). Step (d) is particularly useful when numerous isolated small snow patches are present. A majority filter (replacing the central pixel with the majority of the surrounding values) removes them efficiently without changing the extent of the larger glaciers too much, i.e. only a few pixels at corners are removed (e.g. Paul et al. 2002).

3.4. Post Processing

During post-processing the derived raw glacier outlines (Figure 4e) are quality checked and manually corrected as required (Figure 4f). This is performed with appropriate RGB images (e.g. various colour composites or SAR coherence images) in the background. The major steps are:

(1) Removing wrongly mapped sea ice, ice bergs and water surfaces (lakes, rivers)
(2) Adding missed debris cover and possibly lakes on glacier surfaces
(3) Correcting outlines under clouds (maybe using further images)
(4) Correcting regions in shadow (to the extent possible)
(5) Intersection with drainage divides to obtain individual glaciers

Figure 4: Main processing workflow to derive glacier outlines from a red/SWIR band ratio shown for Oberaar and Unteraar glaciers in the Swiss Alps. a) The region in Landsat Thematic Mapper (TM) band 3 (red), b) the region in TM band 5 (SWIR), c) band ratio (red/SWIR), d) glacier mask after applying a threshold to c), e) outlines from d) after raster-vector conversion, f) overlay of the automatically created glacier outlines (black) and the manually added ones (yellow) on a false colour composite (bands TM 5, 4, 3 as RGB). Numbers in the image refer to the corrections (1) to (4) listed above. All corrections have been applied except for the cloud gap (3). North is at top, image width is 9.5 km.

In regions with abundant heavily debris-covered glaciers and high solar elevations (reducing topographic shading), such as in the Himalaya, SAR coherence images are often the only possibility to identify glacier parts under debris cover (Figure 5). Due to movement and/or alteration of the backscatter properties of the surface, coherence over glaciers (including their debris-covered parts) is usually quickly lost, providing a strong contrast to the much brighter surrounding ice-free regions that are usually not moving. However, water (lakes/rivers) near the terminus and unstable slopes are also subject to a loss of coherence so that a fully automated mapping (using a threshold) is not possible. If not impacted by radar shadow, coherence images can also be used to map clean glacier ice under frequent (orographic) clouds. Although several other methods to identify debris cover semi-automatically have been developed, we think the manual interpretation of coherence images from SAR sensors provides the best results and will thus use them when required and available. The downsides are the larger processing effort for the SAR data, data access restrictions and missing coverage (e.g. for ALOS PALSAR images).

The resulting corrected outlines represent glacier complexes that have to be intersected with drainage divides derived by watershed analysis from a DEM (e.g. Bolch et al. 2010, Kienholz et al. 2013) to obtain individual glaciers step (5). Finally, so-called zone statistics are calculated using the DEM and products derived from it (slope, aspect) as a layer providing the values and glacier extents (polygons) defining the zone over which statistics are calculated, for example topographic information over each glacier. The statistics calculated include minimum, maximum, mean and median values (for slope and aspect only mean) following the methods described by Paul et al. (2009). This calculation is not mandatory for GLIMS but allows data providers to further analyse their dataset. For separated glacier parts that belong together (e.g. regenerated glaciers), the automated method of parameter determination provides values for each entity rather than the combined glacier. The related polygons have to be merged beforehand if topographic data should relate to both parts combined. The resulting multi-part polygons might, however, not fit to the data model of the GLIMS glacier database which needs single part polygons as an input. This issue has yet to be solved.

Figure 5: Example of debris cover mapping with coherence images for glaciers in the Himalaya (centre at 33° 38' N, 76° 12' E). The top image shows clean glacier ice and snow in light blue and clouds in white. The red outlines are mapped automatically, the yellow ones were digitised manually using the coherence image (bottom) as a guide. The degree of coherence is depicted here from white (very high) to black (very low). Dark blue regions in this image refer to radar shadow. North is at top, image width is 42 km.

3.5. Creating a raster dataset from the RGI using a fishnet

For combination with climate model data and global-scale applications, it is important to also have the RGI in raster format. We decided to provide the percentage of area covered for each grid cell as this is the information usually provided for land surface schemes. Creation of such a dataset requires knowledge about the area of glaciers compared to the area of a grid cell for each grid cell. All glacier-specific information is irrelevant, i.e. the attributes and drainage divides can be removed (i.e. all glacier polygons are dissolved). Further, glacier outlines in the RGI and the raw raster dataset are projected in geographic coordinates with WGS1984 datum. These had to be converted to a metric system to get areas in km² (we use Eckert IV projection). The RGI areas per raster cell are calculated after digitally intersecting the dissolved RGI with polygons representing the raster cells. The steps of the processing are described below and are shown in Figure 6 on the example of a 0.1° raster.

(A) Pre-processing of the RGI
1. Dissolve all outlines in each of the 19 regions and combine them to a global dataset (Fig. 6a)
2. Project the outlines to a global-scale metric system that is area preserving (we use Eckert IV)
3. The output file is stored as rgi6_dis_proj.shp

Figure 6:. Steps of the algorithm for glaciers in southern Norway. a) Dissolved glacier outlines from the RGI, b) overlay of a 0.1° by 0.1° fishnet, c) glacier area per grid cell (%) is assigned to the glacier polygons after digital intersection, d) colour-coded glacier area percentage assigned to the 0.1° grid cell (glaciers as an overlay).

(B) Creation of a fishnet dataset with polygon topology
1. Create a 'fishnet' (Figure 6b) covering the region from longitude -180° to +180° and latitude -85° to +85° (net_1d.shp).
2. Project the fishnet to a global-scale metric system that is area preserving (we use Eckert IV).
3. Calculate a running ID (org_id) for each polygon of the fishnet and the area of each polygon cell ({}n_area km2), the resulting file is named net_1d_proj.shp.

(C) Intersection, joining and conversion
1. Intersect A3 (rgi6_dis_proj.shp) with B3 (net_1d_proj.shp) to create cell-specific glacier areas in the file int_rgi_net_1d.shp (Figure 6c).
2. Calculate the glacier area in each cell (area_km2) for the intersect file (int_rgi_net_1d.shp).
3. Calculate the percent glacier cover per cell by dividing area_km2 by _n_area km2_ (percent).
4. Join the percent item from int_rgi_net_1d.shp with the fishnet shape file (net_1d.shp).
5. Convert the fishnet vector file to a raster dataset (rgi6_1d_perc.tif) using the attribute percent as the cell value (Figure 6d).
6. Set 0 to no data and export the file to netCDF format.
7. Check the final products (attributes, correct values).

3.6. Creating a raster dataset from the RGI using a ca. 100 m raster

As the method described in Section 3.5 was very demanding, we decided to also test another method. Here, the glacier extents from each of the 19 glacier regions is converted to a 0.001° (about 100 m) grid and then aggregated to cell sizes of 0.01° (about 1 km), 0.1°(about 11 km) and 1° (about 110 km). The aggregation simply counts the number of cells from the finer resolution mesh, i.e. the 0.01° grid has a maximum of 100 cells from the 0.001° grid and the 0.1° a maximum of 10000 cells. This allows coding of glacier coverage per grid cell from 0 to 100 (with no decimals) for the 0.01° version and from 0 to 100 with two decimals for the 0.1° resolution. The latter would be better for various applications. In a next step, we will also determine the total glacier area (in km²) per grid cell which would be required to transform the glacier mass change product in Gt (which is currently provided by the glacier change service, [RD3]) to specific mass change in m water equivalent (w.e.) per m and year.

The image in Figure 7 depicts the resulting grids indicating (colour-coded) the glacier area in per cent of the grid cell area at 0.01 and 0.1° resolution.

Figure 7: Glacier area per grid cell in percent of the total grid cell area for a region in south-west Switzerland, a) at a resolution of 0.01° and b) at a resolution of 0.1°. The black lines are the vector outlines, shown here for comparison. The scale and orientation is the same for both images.

4. Output data

4.1. Vector dataset

The first output of the glacier mapping algorithm is a vector dataset in shape file format. This format consists of a series of individual files all starting with the same name and different three letter coded endings (Table 3). For exchange across platforms, all components of a shape file are usually compressed into one file (.zip or .gz). The vector outlines derived from satellite imagery have still the rasterized shape of the original images (often with millions of entries to store the coordinates of each corner point) and are smoother where they have been manually corrected. To keep this contrast, it is recommended to not further smooth or generalize the original outlines.

Table 3: The four sub-files required to characterize a shape file.

For the scene-specific processing at a regional scale, our datasets have the projection of the satellite scenes used (UTM with WGS1984 datum). If datasets from different UTM zones have to be merged for practical reasons (e.g. to cover the entire Alps), the central zone (for 3 zones) or those hosting most of the glaciers (for two zones) is chosen. For examples, the Alps cover UTM zones 31, 32 and 33 and the output dataset is thus in UTM zone 32. More than three UTM zones should not be merged to a specific UTM zone as the UTM projection does not conserve area. For the global dataset provided in the Climate Data Store (the RGI) all outlines are projected to geographic coordinates (longitude/latitude) in decimal degrees with WGS 1984 datum. For specific calculations in sub-regions, they might have to be re-projected to a metric coordinate system.

Details of the output fields that are stored in the attribute table (.dbf file) and the csv file containing the hypsometric information are described in the Product User Guide and Specification (PUGS) document [RD2] and are, thus, not repeated here.

4.2. Raster dataset

The resulting raster dataset is exported to netCDF formats. For the version with a resolution of 1° by 1°, the file has 360 by 170 cells ranging from 180° W to 180° E longitude and from 85° S to 85° N latitude. The projection is geographic with WGS 1984 datum. Cells without glacier cover are assigned to 'no data', all others give area in per cent with five decimals.

References

Atwood, D., Meyer, F. and Arendt, A. (2010): Using L-band SAR coherence to delineate glacier extent. Canadian Journal of Remote Sensing, 36, S186-S195.

Bolch, T., Menounos, B. and Wheate, R. (2010): Landsat-based glacier inventory of western Canada, 1985-2005. Remote Sensing of Environment 114, 127-137.

Dozier, J. (1989): Spectral signature of alpine snow cover from Landsat 5 TM. Remote Sensing of Environment, 28, 9-22.

Fischer, M., Huss, M., Barboux, C. and Hoelzle, M. (2014): The new Swiss Glacier Inventory SGI2010: Relevance of using high-resolution source data in areas dominated by very small glaciers. Arctic, Antarctic and Alpine Research, 46, 935-947.

Frey, H. and Paul, F. (2012): On the suitability of the SRTM DEM and ASTER GDEM for the compilation of topographic parameters in glacier inventories. International Journal of Applied Earth Observation and Geoinformation, 18, 480-490.

Frey, H., Paul, F. and Strozzi, T. (2012): Compilation of a glacier inventory for the western Himalayas from satellite data: Methods, challenges and results. Remote Sensing of Environment, 124, 832-843.

Hall, D.K., Chang, A.T.C. and Siddalingaiah, H. (1988): Reflectances of glaciers as calculated using Landsat 5 Thematic Mapper data. Remote Sensing of Environment, 25, 311-321.

JPL (2002): ASTER spectral library Version 2.0. Online at: http://speclib.jpl.nasa.gov/

Kääb, A., Winsvold, S.H., Altena, B., Nuth, C., Nagler, T. and Wuite, J. (2016): Glacier remote sensing using Sentinel-2. Part I: Radiometric and geometric performance, and application to ice velocity. Remote Sensing, 8 (7), doi: 10.3390/rs8070598 .

Kienholz, C., Hock, R. and Arendt, A.A. (2013): A new semi-automatic approach for dividing glacier complexes into individual glaciers. Journal of Glaciology, 59 (217), 913-925.

Linsbauer, A., Huss, M., Hodel, E., Bauder, A., Fischer, M., Weidmann, Y., Bärtschi, H. and Schmassmann, E. (2021): The New Swiss Glacier Inventory SGI2016: From a Topographical to a Glaciological Dataset. Frontiers in Earth Science, 9, 704189, doi: 10.3389/feart.2021.704189.

Paul, F., A. Kääb, M. Maisch, T.W. Kellenberger and W. Haeberli (2002): The new remote-sensing-derived Swiss glacier inventory: I. Methods. Annals of Glaciology, 34, 355-361.

Paul, F., Barry, R., Cogley, J.G., Frey, H., Haeberli, W., Ohmura, A., Ommanney, S., Raup, B., Rivera, A. and T, M. (2009): Recommendations for the compilation of glacier inventory data from digital sources. Annals of Glaciology, 50 (53), 119-126.

Paul, F., Barrand, N.E., Baumann, S., Berthier, E., Bolch, T., Casey, K., Frey, H., Joshi, S.P., Konovalov, V., Le Bris, R., Mölg, N., Nosenko, G., Nuth, C., Pope, A., Racoviteanu, A., Rastner, P., Raup, B., Scharrer, K., Steffen, S. and Winsvold, S.H. (2013): On the accuracy of glacier outlines derived from remote sensing data. Annals of Glaciology, 54 (63), 171-182.

Paul, F. and 24 others (2015): The Glaciers Climate Change Initiative: Algorithms for creating glacier area, elevation change and velocity products. Remote Sensing of Environment, 162, 408-426.

Paul, F., Winsvold, S.H., Kääb, A., Nagler, T. and Schwaizer G. (2016): Glacier remote sensing using Sentinel-2. Part II: Mapping glacier extents and surface facies, and comparison to Landsat 8. Remote Sensing , 8 (7), 575; doi:10.3390/rs8070575.

Qunzhu, Z., Meisheng, C., Xuezhi, F., Fengxian, L., Xianzhang, C. and Wenkuna, S. (1983): A study of spectral reflection characteristics for snow, ice and water in the north of China. IAHS, 145, 451-462.

Toutin, T. (2004): Geometric processing of remote sensing images: models, algorithms and methods. International Journal of Remote Sensing, 25 (10), 1893-1924.

Wulder, M.A., White, J.C., Loveland, T.R., Woodcock, C.E., Belward, A.S., Cohen, W.B., Fosnight, E.A., Shaw, J., Masek, J.G., Roy, D.P. (2016): The global Landsat archive: status, consolidation, and direction. Remote Sensing of Environment, 185, 271–283.

Related articles

• Page:
  Glacier Area document version 2024: Product Quality Assessment Report (PQAR)
• Page:
  Glacier Area document version 2023: Target Requirements and Gap Analysis Document (TRGAD)
• Page:
  Glacier Area data version 7.0: System Quality Assurance Document (SQAD)
• Page:
  Glacier Area data version 7.0: Product Quality Assurance Document (PQAD)
• Page:
  Glacier Area data version 7.0: Algorithm Theoretical Basis Document (ATBD)