Spatial Analysis -An Antarctic Example (for ArcGIS 9.x)

I. Goal

To learn to use tools available in the ArcMap Spatial and 3D Analyst extensions, as they apply to digital elevation models and rasters derived from them, to answer some simple questions and produce attractive maps.

II. Problem

What does Antarctica look like beneath the ice?  A continent of mountain ranges, deep valleys, plains, inland seas, offshore islands and the like exists there, for the most part invisible but for a few features that protrude above the ice.  Wouldn’t it be nice to have a topographic map in shaded relief of Antarctica without the ice and with oceans filling areas that are below sea level?  Wouldn’t it be even nicer to have such a map that accounted for the isostatic rise of the land surface that would occur after the weight of the ice was removed?  What would the continent look like if sea level rose by an amount equal to the volume of the water locked up in ice?  How much ice is there?  Digital data are available to make such maps and answers these questions, as is software to do so.  Let’s have a crack at it!

III. Data

Data for today’s lab came from two web sources:

1. The Antarctic BEDMAP project
2. The Scientific Committee on Antarctic Research (SCAR) Antarctic Digital Database, version 3 (ADD v.3).

BEDMAP provides rasters of Antarctic bedrock elevation and adjacent ocean floor bathymetry, of Antarctic ice thickness, and of Antarctic surface elevations. These rasters are provided in ESRI “Grid” format, ready for immediate use!

ADD v.3 is the source of 1:5,000,000 scale vector files that define the coastline of the continent and its ice shelves, areas of rock outcrop, glacial flow lines and point features with names. These data are stored in ESRI coverages, also ready for immediate use.   Metadata describing how the data were created are available at the web sites and within metadata files viewable in ArcCatalog.

N.B.  The SCAR and BEDMAP data used in this lab are not the most current versions.  Downloads from the links above will provide such, but not at scale of 1:5,000,000.  This lab is written to use the older data sets, which are provided through a link below.

Data are stored in folders and sub-folders with the following formats and file names:

• 5 km x 5 km Raster Grid Files:
  • bed_elev – a raster of orthometric elevations (see the metadata for a description of the vertical datum) for Antarctica bedrock and the surrounding ocean floor.
  • ice_thick – a raster of ice thickness for Antarctic ice sheets and ice shelves.
  • surf_elev – a raster of orthometric surface elevations of Antarctica.
     
• Coverages (scale 1:5,000,000):
  • coast05 – a polygon coverage of Antarctica and the permanent ice shelves.
  • cont05 – an arc coverage showing 500 meter contours of elevation.
  • gflow05 – an arc coverage of flow lines for ice streams and glaciers.
  • rock05 – a polygon coverage of rock outcroppings
     
• Shapefile:
  • SouthPole.shp – a point shapefile that contains the location of the South Pole.
     
• Layer file:
  • For your convenience, vector data have been grouped, symbolized, and stored as a layer file named ADD_v3_vectors_layers.lyr.  Adding this layer file to your map saves having to individually load and create symbology for each of the vector layer.

The data sets for this exercises can be downloaded at https://webspace.utexas.edu/wg3486/geo-web-folders/helper/www/Cutting_Edge.  Copy the entire "Antarctic_data" folder (zipped file, ~28 Mb) to your storage space and unzip it.  Do not move individual files and subfolders - the integrity of the raster data will be destroyed if you do so.  If individual raster files or coverages must be copied, do so with ArcCatalog, not with Windows Explorer.

A separate MS Word document with the questions for this exercise is contained within the zipped folder.

A. Spatial Reference

A glance at the metadata shows that all data are stored in a Projected Coordinate System (PCS) with the following parameters:

PCS Type: Stereographic - South Pole
Units: meters
Latitude of Origin: -90.000000 (South Pole)
Central Meridian: 0.00000 (Prime Meridian)
Standard Parallel: -71.00000
False Easting: 0.00000
False Northing: 0.00000
Datum: WGS84

There is a predefined Projected Coordinate System (PCS) that matches this in ArcMap (i.e. "WGS 1984 Antarctic Polar Sterographic"); one could also be created by modifying an existing PCS. As practice, do this now:

1. Open a new map in ArcMap;
2. Right-click the “Layers” Data Frame heading in the Map’s table of contents (TOC);
3. Select the Coordinate System tab;
4. From the area showing folder icons select the coordinate system on the path Predefined>Projected Coordinate Systems>Polar>South Pole Stereographic.

Note that the parameters of this PCS do not quite match those given above. With the South Pole Stereographic projection selected, modify this PCS by:

5. Clicking the Modify button;
6. Enter a new name (e.g. SCAR Antarctic Projection)
7. Select Stereographic_South_Pole from the drop-down menu of the Projection name;
8. Edit the Standard_Parallel_1 to read -71.000000.
9. Click OK and then the “Add to Favorites” button to make this PCS readily available for later use.

If you can’t get this to work, simply loading any of the Antarctic data into ArcMap will also set the PCS, as the Data Frame will take on the coordinate system of the first file loaded and project everything else with a different PCS on-the-fly (in this case all data have the same PCS so on-the-fly projection is not needed). The PCS could also be imported from one of the data files using the “Import...” button in the Coordinate System tab of the Data Frame Properties.

IV. Procedure

The procedure we will follow involves the following general steps:

1. Create color, shaded relief maps of Antarctica and of Antarctic bedrock elevations. The bedrock elevation data we have includes offshore bathymetry that we would like to remove (or mask, i.e. hide), and both data sets need to be rendered with color ramps to show elevation.
2. Render Antarctic bedrock regions presently below sea level in blue. Use the raster calculator to construct a binary raster which can be symbolized as blue or transparent and overlain on the bedrock elevation map. We’ll also create a zero elevation contour line from the bedrock elevation map to outline shorelines of the regions above sea level.
3. Calculate bedrock topography after it has rebounded from ice removal. Using the raster calculator, we will make a new bedrock elevation raster that accounts for isostatic rise. A new water raster will be constructed and rendered as above.
4. Calculate the volume of water locked up in Antarctic ice sheets and shelves. Use published sea level rises that result from melting of the ice to make a bedrock elevation map showing water in areas that would be below sea level.

We begin first by exploring the data.

A. Explore the Data

1. Open up ArcMap and add the layer file vectors_Layers.lyr, which is on the path //SCARv3/scale1-5M.  Permanent ice shelves and ice tongues are in blue, land is outlined in black with no color, outcropping rock is shown in yellow and South Pole is a red dot, as shown below. 

Figure 1. Antarctic vector layers.

2. Add the "surf_elev" grid.  The default is a display that "stretches" the elevations along a continuous grayscale "Color Ramp", showing sea level (O m) as black and the highest elevations as white, with the remaining 254 shades of gray representing intermediate elevations.  The way in which the intermediate elevations are matched to one of the 254 shades of gray is given by the "Stretch Type" which can be changed through the Layer Properties Symbology tab.  Experiment with the Stretch Type to get a feel for how the data can be displayed differently in grayscale.

Figure 2. Antarctic surface elevation raster, displayed with default symbology.

3. Change the Symbology of the elevation raster to "Classified", click the Classify... button in in the Symbology tab, set the "Classification Method" to Defined Interval with an "Interval Size" of 200 (this will group elevations by 200m intervals), as in the Figure 3 below.  Click OK and select a "Color Ramp" from the Symbology tab (tip: Color Ramp names are accessible by right clicking on the color ramp and turning off the check mark next to the "Graphic View" option.)  The Symbology tab should look like the one shown below.

Figure 3. The classification of elevations by 23 equal intervals of 200 meters.

Figure 4. The symbology tab, with elevations classified at 200 m intervals and symbolized with the Elevation #2 color ramp.  Clicking on the "Label" bar will allow formatting of the labels to show no decimal places, as seen in the Label field above.

The resulting symbolized map should look like that Figure 5 below.

Figure 5. Antarctic elevations, classified at 200 m intervals and symbolized with the Elevation #2 color ramp.  Vector layers overlie the elevation grid.

This is just one of many possible ways to display the elevation raster.  For appearance sake it would be nice to eliminate the irregular box that defines the boundary of the data regions (the area displayed in teal color above, at zero elevation).  We could selective symbolized, by classification, all areas of zero elevation with no color (e.g. using the Exclusion button in the Classification window and excluding 0), but there is a better, more permanent way!

B. Spatial Analyst Extension - setting Analysis Options

1. Turn on the Spatial Analyst extension (Tools>Extensions...; check the box) and display the Spatial Analyst toolbar (Tools>Customize...; check the box).
2. Before beginning any Spatial Analyst procedures always set the Options.   Most functions within this extension create new grids (rasters) by performing operations on old grids.  It is imperative that you set a Working directory, so that these new grids are saved in a place where you can retrieve them.  Lacking a Working directory, they are written to the Windows temp folder, and will be erased when you close your project unless you remember to "Make permanent..." your temporary results!  Do not to forget to do this before beginning any Spatial Analyst session.
    Make a folder in your Antarctic_data folder to store your new elevation grids.  The folder name should be less than 13 characters long, with no spaces or special characters; e.g. use the name "New_Grids".
   From the Spatial Analyst toolbar drop-down menu, select "Options...".and set the Working directory to your newly created folder by browsing to it (See figure 6).
3. An Analysis mask defines a region where an analysis will be performed - any raster cells outside of the analysis mask will be ignored during the analysis and, upon creation of the new raster, will have "no data" values.  An analysis mask can be created from an existing raster (see Desktop Help, analysis mask), or a vector polygon file can be used instead.  We wish to restrict our analysis to the region within the coastline of Antarctica, so we'll use the coast05 polygon coverage as our mask.
    In the Spatial Analyst Option window, use the drop down menu in the Analyst mask line to set the Analyst mask to the 5m coast polygon coverage.
4. The Analysis extent also defines the area of an analysis, but when performing overlay analyses (involving more than one raster) can be set so that results are restricted the region where rasters overlap ("Intersection of inputs"; the default) or the entire region of rasters ("Union of Inputs").  We will accept the default Analysis extent, visible on the Extent Tab of the Spatial Analyst Options window.
5. The Cell size option specifies the raster cell size (resolution) of any new raster created during analysis.  With analyses using two or more rasters, it is always best to set this to "Maximum of Inputs" (the default) so that the new raster does not have cells smaller than any of the input rasters (this is the conservative approach; it does not require resampling of one or more of the input rasters).  Leave this set at "Maximum of Inputs".

Figure 6. Spatial Analyst Options window.

C. Spatial Analyst Extension - Clipping a grid file

To restrict the symbology of the elevation grid to the outline of Antarctica, we will clip the existing grid to the Analysis mask set above.  This will create a new grid that has elevation values only within the mask: the coast line of Antarctica and its ice shelves/tongues.

1. Open the Raster Calculator by clicking the drop-down menu on the Spatial Analyst toolbar and selecting "Raster Calculator..."

Figure 7. The Raster Calculator, with the surf_elev raster loaded.

2. Load the surf_elev grid by double clicking on it in the "Layers" box of the Calculator.  Click "Evaluate" and wait while a new grid is made.  Once completed, the new grid file is loaded into the table of contents and given the file name "Calculation".  Note the absence of zero values beyond the coastline (these cells are transparent).  These are now "no data" cells.  To see them (they are presently transparent)  go to the Symbology tab of this new layer, click the "Display NoData as" button and select a color.  The result, when displayed as 9 equal intervals of elevation, will look something like Figure 8.

Figure 8. The new elevation raster, showing "no data" cells in gray and elevations classified into 9 equal intervals.

4. Undo what you just did (i.e. display nodata as No Color) and change the symbology to show a Standard Deviations Stretched Color  Ramp that is "Cyan-Light to Blue-Dark" and check the "Invert" box to set the lightest colors to the highest elevations. (The color ramp names - e.g. cyan-light to blue-dark - are visible by right-clicking on the color ramp in the symbology tab and selecting "display names".)
5. Rename the new raster "Clipped_elev", right-click on it in the table of contents, select, go to "Data" then "Make permanent..." and give the file a name less than 13 character and save it to your New_Grids folder.  This step assures that the grid has a meaningful file name and is in a place where you can retrieve it.  Before this step, the new grid existed in your New_Grids folder with the file name "Calc1".  By "Making Permanent" you stored this file with a meaningful name.

D. Spatial Analyst Extension - Creating a Hillshade raster

A "Hillshade" is a grayscale rendering that shows shadows and highlights to produce a "shaded relief" map.  Placing a hillshade behind a grid that is partially transparent makes the grid look three dimensional.  To create a hillshade:

1. From the Spatial Analyst drop-down menu select "Surface Analysis", then "Hillshade...".  The Input surface will be your new Clipped_elev grid.  The other parameters in this window can be left alone (most are self explanatory, but see the Help file on Hillshade for details) except for the "Output raster" line.  Before clicking OK, the new raster we are about to create should have a permanent name and location - call it hs_clip_elev and browse to your New_Grids folder.  Now click OK.
2. Move the new hs_clip_elev raster to the bottom of the Table of Contents, symbolize it with a Stretch Type of 2 (or less) Standard Deviations, turn it on, and make the Clipped_elev raster 50% transparent.  The resulting map should now look something like Figure 9.

Figure 9. Antarctic topography, rendered with a blue color ramp that is 50% transparent, with a grayscale hillshade beneath.  The solid light blue areas are ice shelf vector polygons with white outlines that lie on top of the other layers.  Yellow similarly shows rock outcroppings.

E. Spatial Analyst Extension - Antarctic Bedrock Elevations

1. Add the "bedrock_elev" raster.  This raster shows topography beneath the ice and bathymetry of the sea floor to 60 degrees South latitude.
2. As above, clip this raster to the "coast05" polygon coverage.  "Make permanent" this new clipped raster by saving it with the name clip_bed_elev to your New_Grids folder.
3. Create a Hillshade for this new raster, as above, and save it as "hs_bed_elev" into the New Grids folder.
4. Symbolize the new raster with a color ramp; the "Precipitation" ramp, inverted, works well for bright colors.  Set the display to 50% transparent to allow the bedrock elevation hillshade created in step 3 to show below.  You should end up with something that resembles Figure 10.

Figure 10. Antarctic sub-ice bedrock topography, including regions beneath the permanent ice shelves.  Cooler colors are lower elevations, warmer higher. The solid black lines mark the shoreline of the continent and edges of ice shelves.  Small bright yellow spots mark areas of outcropping rock.

F. Spatial Analyst Extension - Creating a binary raster

What parts of the above map are below sea level, what parts above?  We can, of course, symbolize the raster to show this, but it can also be done another way.  From the clip_bed_elev raster we'll compute a binary raster; cells with values of 1 will substitute for cells that are below sea level, those above sea level will be give  a value of "nodata".  We can then render the cells with values of 1 as blue (water) and nodata (land) will, by default, be transparent, making the underlying clip_bed_elev raster visible.  To enhance the appearance of "shorelines", we can also produce and overlay a zero elevation contour.

1. To produce a binary raster, we will use a Conditional Statement in the Raster Calculator.  See ArcMap Help (search on "Conditional Statement") for understanding the syntax of such statements.  The statement to use is con([grid]<=0,1), where [grid} is the clip_bed_elev raster.  Open the Raster Calculator and construct the expression con([bed_el_shlvs] <= 0,1) by a combination of typing and clicking on the calculator buttons and layer names.  Click Evaluate.
2. As above, "Make permanent" the new Calculation raster by saving it with the name "binary_water" to the New_Grids folder.
3. Render water areas blue, make the raster 50% transparent, and move it above the clip_bed_elev raster and hillshade in the table of contents.

G. Spatial Analyst Extension - Creating a contour line

1. From the drop-down menu of the Spatial Analyst toolbar, select Surface Analysis>Contour...  Set the "Input surface" to clip_bed_elev, the "Contour interval" to 4400 (this will produce a single contour, because the highest elevation is 4364 m), the "Base contour" to 0, "Z factor" (useful when x, y units are different from z units) to 1, and "Output features" to a new file name, "bed_elev_zero_contour", in your Bedrock_elevation folder.  Click OK.
2. Symbolize the new contour line in black with a 0.5 line width.  Order the table of contents so the contour line is visible.  Your result should resemble Figure 11.

Figure 11. Antarctic sub-ice topography, showing regions below present sea level in blue.  A suggestion of water depth is provided by the selected color ramp and hillshade raster beneath.  Looks like a terrific place to go fishing.

H. 3D Analyst Extension - Calculating ice volume and area

1. Load the "ice_thick" raster.

2.  Turn on the 3D Analyst extension (Tools>Extensions... check the 3D Analyst box) and display the 3D Analyst toolbar (Tools>Customize... check the 3D Analyst box).
3. Set the Layer on the 3D Analyst toolbar to ice_thick
4. From the 3D Analyst drop-down menu select Surface Analysis>Area and Volume...
5. In the Area and Volume Statistics window, set the Input Surface to the ice_thick raster, the Height of plane to 0, and the the Z factor to 1.
6. Click the Calculate statistics button.

The resulting statistics (fig. 12) give the 2D area (plan view area), the surface area (area of the irregular surface defined by the top of the ice-thickness raster when the base of the ice is assumed to be level, and the ice volume (2D area x sum of all cell values).  The units for the results are in the units of the spatial reference, in this case meters.

Figure 12. Results of ice volume calculation in 3D Analyst.

I. Spatial Analyst Extension - Antarctic topography after ice removal and isostatic rebound

Like a ship lightened of its load, on some time scale melting of the south polar ice cap will result in the rise ("rebound") of the underlying continent. The total amount of rise can be modeled as being directly proportional to the thickness of the ice and the ratio of the density of the underlying mantle to that of the ice. Specifically, for individual raster cells:

(Density of ice / Density of mantle) x (Ice thickness) = Elevation change

Taking:
average density of ice = 0.98 g/cm³
average density of mantle = 3.34 g/cm³
- the density ratio of ice to mantle is thus about 0.2825.

To obtain an elevation raster for the continent that includes this elevation difference, we will:

1. Multiply an ice thickness raster by 0.2825 to obtain elevation change;
2. Add the resulting raster to the bedrock elevation raster to obtain isostatically compensated elevations for the Antarctic continent.

Step 1 multiplies a integer raster (ice thickness) by a floating point value, resulting in a floating point raster.

1. Open the raster calculator, load the ice_thick raster, enter the expression [ice_thick]*0.2825, and click Evaluate.
2. Open the raster calculator, load the newly calculated elevation change raster and add to it the bedrock_elev raster; click Evaluate.
3. "Make permanent" this new raster by saving it with the name bed_el_isost in your New_Grids folder.  You do not need to save the result of step 1.
4. Symbolize the new bed_el_isost raster, create and save a Hillshade, create a zero elevation contour, and make a map like Figure 13.

Figure 13. DEM of Antarctica without ice after isostatic rebound.  Black line is present Mean Sea Level contour.

J. Spatial Analyst - A map of "Greenhouse Antarctica", showing the effects of isostatic rebound and sea level rise

Melting of the south polar ice cap, which contains for about 91.5% of the ice in the world, would raise sea level by about 73 meters; melting all of the ice on the planet would raise sea level by about 80.5 meters (see the literature in the SL_Rise folder in the Antarctic_data folder).  To produce a map like Figure 11 that shows higher sea level we must:

1. Subtract 80.5 meters from the cells of the bed_el_isost raster so that elevations are relative to this higher sea level (raster calculator can do this; c.f. section I)
2. Make a hillshade of this new raster.
3. Create a binary raster of regions above and below sea level (c.f. section D).
4. Create a shoreline (i.e. zero elevation) contour (c.f. section E)
5. Symbolize the results.

Your final product will resemble Figure 14.

Figure 14. DEM of Antarctica after isostatic rebound and a sea level rise of 80.5 meters.  Blue areas are below sea level.  Black line is coastline/zero elevation contour relative to a mean sea level that is 80.5 meters above the present level.  Bright yellow shows areas of present rock outcroppings; red dot is south pole.

You're Done!

M. Helper, Fall 2004; Updated Spring, 2010