Exercise 29d--Working with Digital Elevation Models

Going from here

... to here

Many of you will need digital elevation models or digital terrain elevation datasets (DTEDs) in your work.  Using them is straightforward: each cell contains an elevation from which you can compute all kinds of surface properties.  The main issues with DEMs concern their processing: where to get them, how to read the files, how to massage the files into the needed units of measurement (horizontal and vertical), how to reproject them if necessary, how to mosaic them (put together adjoining tiles), and how to visualize them.  All this needs to be done efficiently without losing too much accuracy in the process.

For more examples and discussion, go to the "Links" pages and refer to the previous course's class notes.

Sources of digital elevation datasets

The United States is nicely covered by three series of DEMs, at three different scales, by maps from the United States Geological Survey.  The USGS offers theme sources (ArcInfo export format) showing the outlines of each map (along with its name) at ftp://edcftp.cr.usgs.gov/pub/data/DEM.  Subdirectories offer DEMs, organized by map name, arranged by extent and scale: 7.5 minutes, 15 minutes, and 1:250,000.

The 7.5 minute series contains grids at a resolution of 30 meters (typical) or 10 meters (rare, but increasing in coverage).  Since one degree of latitude is always 111 kilometers, 7.5 minutes covers about 14 kilometers, which is about 460 grid rows at 30 meter resolution, or 1400 rows at 10 meter resolution.  The number of columns will vary for projected DEMs.
The elevation data for the 7.5 minute DEMs was obtained from the USGS 1:24,000 scale series of topographic maps.

The 1:250,000 series contains grids at a resolution of three arc-seconds (1/1200 degree).  Each grid covers one square degree, thereby requiring 1201 rows and 1201 columns (there is a one-cell overlap with neighboring grids).  These elevation data were obtained from the 1:250,000 scale series of topographic maps.

Some versions of these DEMs contain elevations in feet, others in meters, depending on how they have been processed.  It is very helpful to know, from independent sources, what typical elevations are within a region, just to make sure no error in units of measurement is made.

There are many sources of data for the rest of world, but little of it is yet covered at the level of detail of the 7.5 minute DEMs.  A fairly uniform series of DEMs of the world, at 30 arc-second intervals (roughly one kilometer), is available as the GTOPO30 data set at http://edcdaac.usgs.gov/gtopo30/gtopo30.html.  Within a few years this will be replaced by data from the recent Shuttle Radar Topography Mission, which should result in worldwide coverage (at least between -55 and 60 degrees latitude, roughly) at a 100 meter resolution and much better vertical accuracy.  This is one hundred times the amount of data available in the year 2000.

Many private organizations redistribute DEMs, sometimes in their original form, sometimes reprocessed.  The best place to start usually is with a focused web search, because that provides a good chance of uncovering recent additions to data sets.  Use several key words, like "USGS", "DEM", 'DTED", "DTM", "7.5", and so on, along with names related to the place in which you are interested.

A large and growing amount of data is available at http://www.gisdatadepot.com.  It's a good place to look if you can tolerate too many blatant, ugly, in-your-face animated ads.

The examples on these pages were obtained from the GTOPO30 dataset of Pakistan covering east longitude 70 - 80 degrees and north latitude 30 - 40 degrees.  This dataset comprises 100 tiles (one degree by one degree each) at 30 arc-second resolution.  Therefore, each tile is 121 rows by 121 columns--small enough to practice with, yet so small that many tiles are needed to produce interesting DEMs.

These tiles were imported into Spatial Analyst with the File|Import data source item.  The format is "US DMA DTED".  In the "file name" box we specified "*.dt0" and then selected these four adjoining tiles from 100 listed.  (They are the tiles at 74 and 75 degrees longitude, 35 and 36 degrees latitude.)  Spatial Analyst chose the color schemes randomly.

A larger portion of this dataset--42 tiles--forms the logo appearing at the top of the left-hand frame for these pages.

Digital elevation file formats

The official USGS DEM format, although destined to become an antique, is still widespread.  A lot of GIS software reads some variant of this format directly.  But therein lies the rub: the format has several versions ('A', 'B', and 'C') and the GIS software does not necessarily read all versions correctly.  Spatial Analyst's grid import feature, for example, often goes through the motion of reading a DEM, but produces absolutely nothing.  Other people have wrestled with problems ranging from variations of DEM header data to trouble with how the DEM files are organized into blocks.  These problems are not insurmountable, but they are enough of a nuisance to make one select an alternate format by preference.

A common alternative--this appears to be where the USGS has been going for several years--is the SDTS (Spatial Data Transfer Standard) format.  This format consists of about twenty separate (binary) files containing everything you would possibly want to know about the data, including their projection, how they have been processed, and so on.  The main problem with SDTS is that (to date) a lot of GIS software cannot read this complex format.  You need a translator.

The late Sol Katz wrote translators to convert SDTS DEMs into the USGS DEM format and ArcInfo grid format (the format used by Spatial Analyst).  These translators, sdts2dem.exe and sdts2arc.exe, respectively, are available at many web sites.  They are relatively easy to use: you answer a series of questions, mainly to indicate where the data sources are and to name the output data set.  In the process the sdts2arc program displays the extent of the DEM and its projection in a DOS window.  A nice trick is to copy the text from the DOS window while the conversion is running and paste it into a little metadata file for later reference.

Converting horizontal and vertical units of measure

Converting horizontal units

Unless they have been reprocessed, DEMs are georeferenced by seconds of arc.  (There are 3600 arc-seconds per degree.)  If your software does not support this unit of measurement, you have to perform the conversion yourself.

In Spatial Analyst, the most direct method is to use the Warp request.  This is a complex request, but with patience it can be executed in the "Map calculator" interface.  At its heart is a list of "links".  Each link specifies a point followed by the coordinates it should have.  The conversion from arc-seconds to degrees should simply divide each coordinate by 3600.  This is an instance of an affine, or order-1, coordinate transformation.  It can be completely specified by three starting points that are not all on the same line.  A good choice, then, is to use these three links:

(0, 0) --> (0, 0).    This fixes the origin.
(3600, 0) --> (1, 0).    This divides the first coordinates by 3600.
(0, 3600) --> (0, 1).    This divides the second coordinates by 3600.

The full details of the Warp request are available in the ArcView help under "Warp".

Suppose the imported DEM, in arc-seconds, is shown with a theme named "My DEM" in a View.  Here is the full Map Calculator expression for converting it to decimal degrees:

[My DEM].Warp({Line.Make(0@0, 0@0), Line.Make(3600@0, 1@0), Line.Make(0@3600,0@1)}, 1, #GRID_RESTYPE_NEAREST, NIL)

This expression can easily be reused by copying it and pasting it into as many Map Calculator dialogs as necessary, one for each source DEM.  Simply change the theme name appearing between square brackets [].

Because the resulting grid themes use different coordinates, it is a good idea to cut them out of the view and paste them into their own view.

Converting vertical units

Changing the elevation units is much simpler, because it simply rescales each cell value.  For example, to convert elevations from meters to feet, multiply each value by 39.37/12.  Almost all raster-based GISes provide a capability for this kind of operation.  In Spatial Analyst, use the Map Calculator to perform the multiplication.  You could even combine the elevation change with the Warp request to reduce the number of intermediate grid datasets you produce.

Reprojecting grids

If your DEM is in decimal degrees, then it probably needs to be projected.  Otherwise, many values derived from it, such as slope, aspect, and curvature, will be incorrect.

High-end GIS software, such as ESRI's ArcInfo, will reproject gridded data sets.  There also exist third-party commercial products, typically costing around $500 (US), to do this.

Spatial Analyst does not reproject grids.  However, its Warp request can be pressed into service to achieve an excellent (polynomial) approximation in almost every practical case.  (It will have problems for grids spanning 90 degrees or more in any direction.)  This is available as the Reproject Grids extension at ESRI's Arcscripts pages.

When using approximate projections, as with the Reproject Grids extension, it is probably more accurate to reproject the tiles before mosaicing them.  However, by mosaicing them first you may avoid the need to repair gaps later on (see below).

These are the separately reprojected DEMs.  The projection is UTM-1983, Zone 42.  The projection was approximated using affine (order-1) transformations, which introduce errors of up to 300 meters in the X direction: about one-third of a grid cell.  (Order-2 transformations would reduce the maximum error to one meter; order-3 transformations are accurate to a few millimeters.)

The NoData values have been made transparent to keep the borders of any grid from obscuring data in its neighbors.

Mosaicing grids

Most of any map is close to an edge.  This perverse rule is borne out by the graph showing how much of a map's area (X-axis) lies a given percentage of the map's dimensions from an edge. 

This graph is based on a simple geometric calculation.

If the map dimensions are l and w, it has area lw.  The area of the map that is not within a fraction y of the edge is the area of a rectangle with dimensions l(1-2y) by w(1-2y), which is lw(1-2y)2.  Therefore the area within a fraction y of an edge, relative to the original area, is 1 - (1-2y)2.

The graph, therefore, is given by border area = x = 1 - (1-2y)2.  Equivalently, distance from edge = y = (1-sqrt(1-x))/2.

For example, fully one-third of a map lies within 10% of an edge.  If the map is 24 inches wide by 36 inches tall, this means that one-third of the map will be within 2.4 inches of the left or right border or within 3.6 inches of the top or bottom border (or both).  If you need to work with features extending more than a few inches across at this scale, then you will frequently need data from more than one map sheet to show them.

It follows that you will often need more than one DEM to cover your region of interest.  For some purposes, such as visualization, it is ok to maintain the DEMs separately, provided you use the same classifications and symbols to visualize each theme.  However, it rapidly becomes inconvenient when those DEMs are being used for further analysis, because then every operation has to be separately performed for each DEM.  Worse, many grid operations (such as hillshading) cannot be performed on the very borders of the grid.  This is rarely a problem, but when the grid is effectively tiled into adjoining pieces, these missing borders become lines of missing values criss-crossing the map.

The solution is to mosaic the DEMs into a single grid.  Mosaicing fits the tiles together into a whole, typically averaging the values in regions where they overlap.  So long as the DEMs involved use the same elevation units, have similar cell sizes, are of comparable quality, and overlap only slightly around their edges, this is a simple and trouble-free operation.

The Spatial Analyst syntax for mosaicing is straightforward, albeit a little strange.  It requires a Map Calculator expression.  That expression must reference each grid theme to be mosaiced.  It does so by putting the grid names (enclosed in the usual square brackets  [] ) into a list.  A list is constructed simply by enclosing its elements in curly braces  {}  and separating them by commas, as in {1, 2, 3} or {"a", "b", "c"}.

The strange part is that you need to identify one grid theme as the starting theme for the operation.  It does not matter at all which one you use.  This will be the grid to which the Mosaic "request" is applied.  The list of the remaining grids is the "argument", or parameters, for this request.

In our case we have four grid themes, "grid1", "grid2", "grid3", and "grid4".  One possible mosaic expression for the Map Calculator is therefore

[grid1].Mosaic({[grid2], [grid3], [grid4]})

By paying attention to all the brackets, and remembering why each is there, you will get this right every time.

Sometimes you have to mosaic a lot of grids.  To construct a DEM of Pakistan using the GTOPO30 data, for example, you would need around 300 tiles from four datasets.  The mosaicing will have to be done in stages, successively building larger pieces out of smaller, because operations involving more than a few dozen grids at a time are (by a corollary to Murphy's Law) destined to go awry.  (Spatial Analyst limits the number of grids in one operation to 50, anyway.)

It will still be very painful to type syntactically correct expressions that mosaic dozens of grids.  Here's a Spatial Analyst trick that will work with other GISes (although the details will be different): write a tiny program to generate the mosaic expression for you.

Here is how it's done.  First, you need to provide systematic names for the tiles, like "E074N35" (as in the GTOPO30 data) or "Grid1", "Grid2", and so on.  (If you're not already doing that, you shouldn't be working with such large data sets!)

Now write the program.  In Avenue (ArcView's language) it might look like this:

s = "[Grid1].Mosaic({"
sComma = ""
for each i in 2..20
    s = s + sComma + "[Grid" + i.AsString + "]"
    sComma = ", "
end
s = s + "})"
MsgBox.Report(s, "")

In this script, "+" is the string concatenation operator and "MsgBox.Report" displays a string in a window from which the string can be copied.

Every such script has five characteristic features:

  1. It creates the initial part of the expression.
  2. A loop constructs the repeated parts of the expression.
  3. There has to be some kind of trick to prevent the common separator (a comma in this case) from appearing before the first item or after the last.  This script illustrates one method: the separator starts out empty and then is reset to the correct value after each iteration of the loop.
  4. The script creates the final part of the expression.
  5. It displays the expression.

The output of this particular script is

[Grid1].Mosaic({[Grid2], [Grid3], [Grid4], [Grid5], [Grid6], [Grid7], [Grid8], [Grid9], [Grid10], [Grid11], [Grid12], [Grid13], [Grid14], [Grid15], [Grid16], [Grid17], [Grid18], [Grid19], [Grid20]})

That is a valid Map Calculator expression to mosaic the 20 grid themes "Grid1", ..., "Grid20".  In ArcView, you can use a variation of this technique to automate the entire mosaicing process, thereby avoiding all that painful typing.

 

This is the mosaic of the four reprojected grids.  Now every cell is automatically shown using a consistent method of symbolization, which is already an advantage.

The four reprojections were four separate transformations.  Each one provided the best approximation for its tile.  Because of this, there is a slight "tearing" or gap between two tiles visible in the middle bottom as a vertical dashed white line.

This shows a five-fold enlargement of the gap.

Patching gaps

"Tears" or gaps in a mosaiced image result when the tiles fail to overlap completely at their edges.  This usually occurs because all projections distort distances differentially over large areas.  Regions that are perfectly tiled in decimal degree coordinates therefore might not overlap in places after they have been separately projected.

There are two good ways to patch over these gaps.

Estimate the gap values from the neighboring elevations.
Estimate the gap values from an alternative DEM, usually of coarser resolution and greater spatial coverage.

A combination of the two techniques can be used as well.  The simplest approach, which is very effective for gaps only one or a few cells thick, is to use a neighborhood average.  This is computed in two steps using Spatial Analyst.

  1. Compute the neighborhood average of the mosaic.  First, inspect the gaps to determine the widest part.  In the example above, the gap never gets wider than one cell.  Therefore, we will use a three-cell neighborhood, which will include one line of cells to the right and one line of cells to the left of each gap cell.  The neighborhood average is computed using the Analysis|Neighborhood statistics item.

  2. Patch the gap.  Use the Map Calculator for this.  The best method is with the Con (conditional) request.   Assuming the name of the grid theme to patch is "Mosaic" and the name of the neighborhood average is "Average", then the appropriate expression is

[Mosaic].IsNull.Con([Average], [Mosaic])

This expression first tests for gaps (IsNull).  The Con request will replace any gap cell with the value in the corresponding cell of [Average] and otherwise will keep the value from the [Mosaic] grid.

This is the patched mosaic, shown using 256 classifications.  (The color scheme was created with the ColorRamp extension on ESRI's ArcScripts pages.)  The grid is 240 rows and 200 columns.  There are fewer columns because the reprojection compressed east-west directions by the cosine of 35 degrees (about 0.82).  This would suggest only 196 columns are needed, but the slight rotation increases the width.

The same approach works when using an alternative DEM: it plays the role of the "Average" grid in the patching.  A typical Map Calculator expression would be

[Mosaic].IsNull.Con([DEM], [Mosaic])

Visualization using hillshading

A digital elevation model represents a continuous surface.  Our eyes receive many visual cues from real surfaces according to how they are lit.  To provide comparable cues in a map, hillshading estimates the intensity of light impinging on the surface at each grid cell.  To be a fast computation, hillshading does not attempt to account for shadows cast by other parts of the surface: the estimate is purely local to each cell, depending only on the angle the surface makes to an imaginary light.

Computing a hillshade requires a bit of trigonometry, but is otherwise fast and easy.  (Spatial Analyst offers a menu item, Surface|Compute Hillshade, that makes this simple.  You need only specify the position of the light.)

Hillshading helps when the software can adjust the lightness of the map according to the hillshade value.  Most raster-based GISes will do this.

Spatial Analyst provides the interface in the "advanced" button of the legend editor.  After computing the hillshade for a grid, specify this hillshade in the dialog and provide some scaling parameters for the relative lightness.

I always specify a maximum of 100 percent because anything less unnecessarily reduces contrast.  Usually, you need as much contrast as possible when visualizing elevations.  If you need a uniformly darker map, darken the colors in the legend instead.

The thing to watch out for with hillshading is the unit of measurement.  To compute the light intensity correctly, the software needs to relate elevations to distances (it is computing slopes, among other things).  If your elevations are not in the same units as the horizontal distances, you are likely to get strange or unusable results.  The solution to this problem is temporarily to rescale all elevations using a Map Calculator expression, compute the hillshade of the rescaled grid, and apply that hillshade to the original grid.

The hillshaded, patched, mosaiced, reprojected DEM, ready for use.

The classifications and colors are exactly the same as before.  The difference is due solely to the hillshade.

Update history:

20 July 2001: Added example of patching a gap with an alternative DEM.  Added an explanation of the area-distance graph in the Mosaicing grids section.  External links now open new windows.  Minor improvements made in the text.
22 February 2001: Typographical corrections.
9 April 200:

Created.