ASF Logo
North American Space Agency (NASA) Logo

Alaska Satellite Facility - Distributed Active Archive Center

How to Create Glacier Velocity Maps with Sentinel-1 Toolbox

Rink Glacier velocity vectors. Contains modifed Copernicus Sentinel data 2016, processed by ESA.

Adapted from the European Space Agency’s STEP community platform

●Easier

Background

The goal of this tutorial is to provide novice and experienced remote sensing users with step-by-step instructions on the use of Offset Tracking tools in generating glacier velocity maps with Sentinel-1 Level-1 Ground Range Detected (GRD) products. Offset Tracking is a technique that measures feature motion between two images using patch intensity cross-correlation optimization. It is widely used in glacier motion estimation.

Level-1 GRD products consist of focused SAR data that has been detected, multi-looked and projected to ground range using an Earth ellipsoid model such as WGS84. The ellipsoid projection of the GRD products is corrected using the terrain height specified in the product general annotation. The terrain height used varies in azimuth but is constant in range.

This tutorial will examine the movement of the Rink glacier. Rink glacier is a large glacier located on the west coast of Greenland. It drains an area of 30,182 km2 (11,653 sq mi) of the Greenland ice sheet with a flux (quantity of ice moved from the land to the sea) of 12.1 km3 (2.9 cu mi) per year, as measured for 1996. It is also the swiftest moving and highest surface ice in the world.

Rink Glacier in Google Maps
Rink Glacier calving; Credit: Faezeh M. Nick

Prerequisites

Materials List

Note: The input products should be two GRD products over the same area acquired at different times. The time interval should be as short as possible. In this tutorial, we will use two GRD products acquired 12 days apart.

System Requirements

Many of the steps within this data recipe may take some time to process. We recommend the following:

  • At least 16 GB memory (RAM)
  • Close other applications if possible while using the Sentinel-1 Toolbox
  • Do not use the computer during processing to avoid crashes

Steps

Open the Products in Sentinel-1 Toolbox

    1. Open the Sentinel-1 Toolbox software
    2. Use the Open Product button in the top toolbar and browse to the location of the two GRD products
    3. Press and hold the Ctrl button to select both .zip files and press Open

Note: Do not extract the files, Sentinel-1 Toolbox will do so automatically when the files are loaded into the program. If you unzipped the files, select the manifest.safe file from each Sentinel-1 product folder.

View the Products

    1. In the Product Explorer, you will see the opened products
    2. Double-click on the product to expand
    3. Expand the Bands folder and you will find two bands: an amplitude band, which is actually the product, and a virtual intensity band, which is there to assist you in working with the GRD data.
    4. To view the data, double-click on the Amplitude_HH band.
    5. Zoom in using the mouse wheel and pan by clicking and dragging the left mouse button.

Comparing the Amplitude_HH Band of Product image to the Products Explorer with World View image, you will find that the image is flipped upside-down. This is because the image was acquired with an ascending pass and right-pointing antenna. Therefore, the bottom part of the image area was first sensed and in Sentinel-1 Toolbox, the first sensed line is always displayed on top of the image.

Product Explorer with World View
Amplitude_HH Band of Product

Apply Orbit File

    1. Select GRD product in the Product Explorer window
    2. From the top menu navigate to Radar > Apply Orbit File
    3. In the Apply-Orbit-File window, specify the output folder and the target product name. The products will automatically be appended with the suffix _Orb if you choose the default name.
    4. Leave all default parameters and click Run
    5. Repeat steps 1-4 for the second GRD product
    6. When finished, close the Apply-Orbit-File window
Navigate to Apply Orbit File
Apply Orbit File dialog

Coregister the Images into a Stack using DEM

From the top menu, navigate to Radar > Coregistration > DEM Assisted Coregistration and select DEM Assisted Coregistration with XCorr.

    1. In the ProductSet-Reader tab, press the “+” button to add the products, ensuring that the older product is selected first.
    2. In the DEM-Assisted-Coregistration tab, select the Digital Elevation Model (DEM) to use, the DEM resampling method, and image resampling method.

Note: The default DEM is SRTM 3 Sec, which covers most area of Earth’s surface between -60° latitude and +60° latitude. However, it does not cover the high latitude area where Rink Glacier in the sample granule is located. Therefore, ASTER GDEM, GETASSE30 or ACE30 DEM could be selected. Areas outside of the DEM or in the sea may be optionally masked out.

    1. In the Write tab, specify the output folder and the target product name. The product will automatically be appended with the suffix _Stack if you choose the default name.
    2. Leave all other parameters as default and click Run
    3. When finished, close the Coregistration window
DEM Assisted Coregistration with paired products
DEM Assisted Coregistration with DEM selected
Coregistered stack product

Create a Subset Image Containing Rink Glacier

Since the image covers a large area of the west coast of Greenland and we are interested only in the Rink Glacier area, we will create a subset of the coregistered stack that contains the Rink Glacier area only.

    1. Open one of the bands from your coregistered stack
    2. Zoom in the image until the image window contains only the area of interest
    3. Right-click on the image and select Spatial Subset from View…
    4. When the Specify Product Subset window appears, click OK to create the subset
    5. Save the newly created subset product by right-clicking on the product in the Product Explorer and select Save Product
Context menu over area of interest with Spatial Subset from View selected.
Subset area in the image

Offset Tracking

The Offset Tracking operator estimates the movement of glacier surfaces between master and slave images in both slant-range and azimuth direction. It performs cross-correlation on a selected Ground Control Point (GCP) in master and slave images. Then the glacier velocities on the selected GCPs are computed based on the offsets estimated by the cross-correlation. Finally, the glacier velocity map is generated through interpolation of the velocities computed on the GCP grid. The Offset Tracking is performed in the following substeps:

    1. For each point in the user-specified GCP grid in the master image, compute its corresponding pixel position in the slave image using normalized cross-correlation.
    2. If the compute offset between master and slave GCP positions exceeds the maximum offset (computed from user-specified maximum velocity), then the GCP point is marked as an outlier.
    3. Perform local average for the offset on valid GCP points.
    4. Fill holes caused by the outliers. The offset at a hole-point will be replaced by a new offset computed by local weighted average.
    5. Compute the velocities for all points on GCP grid from their offsets.
    6. Finally, compute velocities for all pixels in the master image from the velocities on the GCP grid by interpolation.

In the Offset Tracking dialog window, the user needs to define a GCP grid by specifying the grid-point spacing in range and azimuth directions. The spacing is specified in terms of the number of pixels. Then the system will convert the spacing into corresponding spacing in meters and compute the grid dimension.

The user also needs to specify other processing parameters such as Registration Window dimension and Maximum Velocity. Before running the operator, the user is suggested to do some research on the maximum velocity of the glacier under study on the season of acquisition. This will be helpful in guiding the user selecting meaningful processing parameters. For example, the maximum velocity for Rink Glacier is around 10 meters per day. Given that the SAR image acquisition period is 12 days and range and azimuth spacing is 10 meters, we can calculate the maximum shift of a target in the glacier is about 12 pixels. We know that the default Registration Window dimension (128 pixels) is large enough to cover the target in both master and slave images.

The Spatial Average and Fill Holes processing steps can be optionally turned off by deselecting the corresponding checkboxes in the dialog window.

Generate Glacier Velocity Map

    1. Navigate to Radar > SAR Applications > Offset Tracking from the top menu
    2. Select your subset image as the source product
    3. Specify the output folder and the target product name. The product will automatically be appended with the suffix _Vel if you choose the default name.
    4. In the Processing Parameters tab, change Max Velocity (m/day): to “10”
    5. Leave all other default parameters and click Run
Offset Tracking dialog

View the Glacier Velocity Map

    1. Double-click on the velocity band in the resulting product to display the velocity map
    2. Navigate to Layer > Layer Manager
    3. From the Layer Manager window, deselect Vector data to remove the grid
    4. Click on the “+” button to open the Add Layer window
    5. In the Add Layer window, click Coregistered GCP Movement Vector and click on Finish. You will see the velocity vectors displayed on the GCP grid showing direction and speed.
Velocity Map of Rink Glacier; Contains modified Copernicus Sentinel data
Velocity Map of Rink Glacier; Contains modified Copernicus Sentinel data
Add Layer dialog

Sample Image

Rink Glacier velocity vectors. Contains modifed Copernicus Sentinel data 2016, processed by ESA.