Processing Plan for GLIMS

Bruce H. Raup

NOTE: The following document is old and obsolete in places. In particular, the form of the GLIMS Glacier ID has changed (although the concepts behind it have not changed). Please see the data dictionary for the correct information. Also, only lat/lon/WGS84 coordinates are stored in the GLIMS Glacier Database.

Contents

• Introduction
• Single-Image Processing
  • Human-intensive tasks
  • Computer-intensive tasks
    • Creation of Glacier ID
  • First steps upon image acquisition
• Multi-Image Processing
  • Image Registration
  • Displacement Field Measurement
• About this document ...

Introduction

This document outlines the plan for data processing within the Global Land Ice Measurements from Space (GLIMS) project. The primary source of data for this project will be the ASTER instrument aboard NASA's EOS AM-1 platform, due to be launched in June of 1998.

There will be two types of processing within GLIMS: single image processing, and multi-image processing using time-sequential imagery. These are addressed separately below. The guiding principle behind the processing steps is speed of processing. Our experience to date indicates that a trained human is required to monitor or guide the processing during at least one stage. Thus, we want to use human input not only for tasks that the computer can't do, but also for tasks that a human can do faster and more accurately than the computer.

Single-Image Processing

Human-intensive tasks

Of the many tasks to be performed on an image, the initial ones of classification (finding the glacier(s) in the image), discrimination between bare ice and snow, and deciding on the flow direction of the glacier(s) involve pattern recognition -- something humans do with ease and computers do with difficulty. These initial tasks are:

1. Classify image regions into ice, snow, rock, soil, and vegetation (supervised by a human operator).
2. Draw a polygon around each glacier in the image. (supervised by a human operator).
3. Draw an arc along the centerline of the glacier, from upstream to downstream (probably done by a human operator, or automatically if elevation information is present (e.g. from ASTER Band 3 stereo)).
4. Draw a polygon around regions unsuitable for finding registration tiepoints. (Tiepoints will then be found in all regions except these regions and the regions marked as glaciers.) (This is again is supervised by a human operator.)

Computer-intensive tasks

After these steps have been carried out, the computer can automatically calculate:

1. a unique ID for glacier from average L/L position,
2. total glacier area,
3. area of accumulation zone,
4. glacier length (along centerline arc),
5. width as function of position along centerline,
6. orientation of glacier,
7. multipoint arc describing the position of the glacier terminus (in three dimensions if a DEM is available; otherwise only in the map plane),
8. multipoint arc describing the position of the snow line (in three dimensions if a DEM is available; otherwise only in the map plane),

The steps for single-image processing could therefore go like this:

1. Image is received from data source.
2. Database record is created for image, containing ID, location information, etc (see Section 4 for a full description of the fields).
3. Some combination of operator and computer effort yields polygons surrounding all glaciers, arcs down the centerlines of all the glaciers, polygons defining ablation zones and source snowfields, and polygons surrounding regions unsuitable for tiepoints. This step could largely be done using supervised computer classification if the image is multi-band.
4. Computer then calculates quantities listed above, and stores the resulting data in the GLIMS database.

Creation of Glacier ID

The ID associated with each glacier will be a key that links the data in several parts of the database. The method used for assigning IDs should be fast and unambiguous. Use of an ID based on the World Glacier Monitoring Service (WGMS) system of counting drainages was considered. However, this would not work well in all areas, such as Greenland. An ID based on the L/L position of the center of the glacier is readily understandable by humans, well-defined over the entire earth, and is easily generated. (WGMS IDs, where defined, will be incorporated into the database.)

The uncertainty in the image position (for ASTER), based on orbit geometry, is 455 m per axis (3-sigma) (150 m for satellite position, plus 315 m for uncertainty in pointing) (from the ASTER document End-to-End Data System Concept). This is sufficient to allow unique IDs based on L/L since glaciers that are less than 455 m apart (center-to-center) are rare. (I have not found any.) For example, the small glaciers in the Pyrenees are on the order of 1 km apart.

The computer could find the L/L coordinates of the center of the polygon circumscribing the glacier by averaging the coordinates of the polygon, to come up with a single L/L position for the glacier. An ID could then be created that has the form:

G_Edddddd[N|S]ddddd

where G marks the ID as a glacier ID, E stands for East longitude, and [N|S] means either N or S is to be used, meaning North or South latitude, respectively. The values represented by the "d"s are in decimal degrees, with it understood that there are three digits to the right of the implied decimal point. For example, if there were a glacier in Flagstaff, it would have the ID G_E248380N35189 (since Flagstaff is at 35.18900 degrees N, 111.61976 degrees W). Note that we could use W as well for west longitudes, but it probably will be more efficient and less confusing to require the use of east longitudes ranging from 0 to 359.999. Leading zeroes are included, so that all IDs are the same length.

First steps upon image acquisition

There are two possible scenarios after acquiring an image that contains a glacier. The glacier could be new and not be in the database, or it could be that the glacier has already been entered and assigned an ID. The steps involved for each of these are described below.

New image, New glacier:

1. Get new image. We know rough L/L of center, plus orientation, so we can calculate rough L/L for each pixel.
2. Draw polygon around a glacier in the image. We now know L/S coordinates of the polygon.
3. Average the L/S coordinates of the polygon to get a rough position of the glacier. Create ID [glacierID] based on the corresponding L/L position.
4. Do we want to use this point, or some other, as the ground reference point? If some other, find that point in the image to subpixel accuracy (if possible). Save offsets between ID and reference point in L/S coordinates (temporarily, not in the database).
5. Create the transform (``Transform 1'') between L/S and N/E, mapping the L/S location of the reference point into (0,0) N/E. Store this transform in the database (in the Transform 1 table).
6. Using this transform, convert offsets between ID and reference point to N/E coordinates. Also convert L/S coordinates of the polygon(s) to N/E. Store all these in the database only in N/E coordinates.
7. Store L/L coordinates (as well as they are known) of the reference point in the database (in the Glacier table). These coordinates can be improved as outside information allows. (This L/L pair, together with standard mapping conversions, constitutes ``Transform 2.'')

New image, Old glacier:

1. Get new image. We know rough L/L of center, plus orientation, so we can calculate rough L/L for each pixel in the image.
2. Select glacier; have computer get L/S coordinates of a point near the center of the glacier. Convert this to a rough L/L pair.
3. Search database for glacier IDs that are near this L/L location. (``Near'' means within some settable tolerance.)
4. If no IDs are found, this is a new glacier. (Follow that procedure.) If more than one ID is found, display these images and have operator find the correct one, or simply pick the one closest to the selected point.
5. Now we have the ID for the glacier. Look up its reference point in the database.
6. Look up images in the database that contain the glacier; look up images in the database that contain the reference point. (These might be the same image.)
7. To locate the reference point in the new image, we must register the new image to an old image containing the glacier. This registration defines a transform between L/S coordinates of the two images.
8. Using the ``Transform 1'' relating L/S coordinates of the old image and the glacier's reference point, find the L/S location of the reference point in the old image. We can now find the L/S location of the reference point in the new image.
9. Create a ``Transform 1'' relating L/S coordinates in the new image to the reference point and store in the database.

There is one other case: Two old images

1. Find the same reference point in each image by steps 2-5 above.
2. Look up ``Transform 1'' for each image-reference pair.
3. Link these two ``Transform 1''s together to form a Transform from one image to the other:

   (L/S)1 <==> (N/E)ref <==> (L/S)2

Multi-Image Processing

The primary tasks for mult-image processing are measurement of the velocity field of the glacier, and measurement of changes in position and areal extent of the ice. Subtasks include:

1. Registration of two time-sequential images of the same geographical area.
2. Automatic or semi-automatic determination of velocity field and mean terminus position.

Image Registration

From the orbit data we will know the geolocation of the images to within 150 m, or about 10 pixels. Thus at worst we will know the image-to-image registration to within 20 pixels. The computer will use cross-correlation of several subscenes to find matches between the two images, then generate transformation coefficients.

The first time we acquire an image, the image regions will be classified into ice, rock, water, etc as described in Section 2.1. Upon acquisition of a second image of the same region, the computer will look up from the database regions suitable for finding tiepoints. If such are found, then steps 1 through 3 below can be skipped. In some regions, such as the West Antarctic Ice Streams, there is no exposed bedrock, and more sophisticated and labor-intensive co-registration procedures will probably be required.

The steps for the registration process are:

1. Display both images on the computer screen.
2. User clicks on a point in the first image, and possibly on the corresponding point in the second image (to save computation time).
3. Repeat step 2 at least once.
4. Computer uses cross-correlation to find subpixel location of point selected in the first image (or found automatically) in the second image.
5. Computer uses these matches to create transform that maps one image into the other.

Displacement Field Measurement

The general steps for measuring the displacement field are:

1. Computer checks database for previously-obtained velocity information for this glacier. If it exists, one or more ``seed displacement vectors'' are retrieved from the database, and the next step is skipped.
2. Operator clicks on a surface feature of the glacier in the first image, then the same feature (which has moved with the glacier) in the second image. This serves as a ``seed displacement vector'' on which the computer can base an automatic search. (This step could be skipped if the GLIMS database contains seed points already.) Alternatively, one or both images are analyzed automatically to find areas on the glacier which will probably yield good correlation matches.
3. The computer uses cross-correlation to determine the precise displacement vector for the seed(s).
4. The computer expands out from each seed in some fashion (see below), finding as many displacement vectors as possible.
5. Computer stops this search when at least one criterion is met, such as arriving at the boundaries of the glacier, or expanding into a region where cross-correlation is impossible or of dubious quality.
6. Computer weeds out bad vectors based on consistency checks.
7. Operator may need to check or further edit the found vectors.
8. Computer adds these displacement vectors and estimates of their uncertainty to the GLIMS database.

Three desires for the search are 1) that the computer uses nearby already-computed vectors as ``guesses'' when working on the current vector, 2) retain the possibility of enhancements such as warping search chips according to the strain-rate field, and 3) allow the computer to track automatically areas of strong cross-correlation.

The method currently implemented spreads out from each seed vector horizontally, filling in each row as completely as it can. After the search has spread until it fails to correlate in all directions, a polygon is created to outline the region of good matches. The next seed vector in the list is then tried if it falls outside all previously-calculated polygons. For each new point, the average of nearby vectors which were previously calculated is used as a guess for the displacement of the current point. Details of this algorithm are described in a separate document.

About this document...

Processing Plan for GLIMS

This document was generated using the LaTeX2HTML translator Version 96.1 (Feb 5, 1996) Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.

[...and considerable editing was required afterward. -BR]

The command line arguments were:
latex2html -no_navigation -split 0 plan.

The translation was initiated by Bruce Raup on Tue Oct 15 15:45:11 MST 1996