DEM Processing for Hydrologic Modeling Studies

Matthew M. McPherson and Heather E. Henneman

 

Abstract

The US Army Corps of Engineers' Hydrologic Engineering Center and Sacramento District prepared USGS 30-meter DEMs for a study of the 59,000 square mile Sacramento/San Joaquin basin. The authors' experiences show that the quality and accessibility of GIS data for hydrologic modeling are improving, as are the hardware and software necessary to process it. However, the day of bringing terrain-based modeling to the great, unspecialized masses hasn't yet arrived. Problems to overcome include downloading and processing around 3000 7.5-minute SDTS DEMs, with disparate units and vertical datums, and manually screening and editing them into shape for drainage analysis.

Table of Contents

  1. Introduction
  2. Data Acquisition
  3. Source DEM Characteristics
  4. Target DEM Characteristics
  5. Intermediate Grid Processing and Assembly
  6. Special Processing
  7. Final Assembly
  8. DEM Inspection and Correction
  9. National Elevation Dataset (NED)
  10. Conclusions

Acknowledgements

Appendices

Appendix A - Downloading DEM Data

Appendix B - sdtsdem.aml

Appendix C - Useful Arc/Info GRID Commands

References

Author Information

 

Introduction

The Sacramento River watershed in Northern California covers approximately 27,000 square miles, while the San Joaquin/Tulare River watershed has a drainage area of 32,000 square miles. The US Army Corps of Engineers, Sacramento District, and the Hydrologic Engineering Center (HEC) conducted a Comprehensive Study of these basins that included the development of hydrologic models for the watersheds (Dunn and Doan, 2000). Detailed hydrologic models were created for approximately 30,000 square miles of the project area shown in Figure 1.

Figure 1. Comprehensive Study Project Area

Methods to delineate watersheds and extract topographic, topologic and hydrologic information from digital terrain data have been developed for use with desktop GIS technology (Djokic and Ye, 1999; Olivera and Maidment, 1999). An ArcView extension called HEC-GeoHMS, developed by Esri and HEC, incorporates such methods to provide schematics, basin geometry, routing information and precipitation grid references for the HEC software program HEC-HMS, Hydrologic Modeling System (Doan, 2000). GeoHMS requires a digital terrain grid as input; therefore, a digital elevation model (DEM) covering the 59,000 square mile project area was created.

This paper details the creation of the digital elevation model used in the Comprehensive Study. The acquisition of the data is described, including instructions for downloading US Geological Survey (USGS) 30-meter DEM data. Characteristics of the source data and the target DEM are outlined, and the procedures for processing, assembling and editing the data are given.

[Back to Table of Contents]

 

Data Acquisition

DEM Background

Project time constraints dictated that the DEM data be available immediately. The size of the study area, combined with a preference that data and methods used throughout the study area be as homogeneous as possible, further restricted the potential data sources.

In the United States the most widely available digital elevation models (DEMs) are those distributed by the USGS. The USGS 7.5-minute DEM data have a spatial resolution of 30 meters, i.e., the grid cells are 30-by-30 meters. Recently, the USGS began producing 7.5-minute DEM data with a spatial resolution of 10 meters; however, although the 10-meter DEM coverage continues to expand, the product still only sparsely covers the United States. The USGS 30-meter DEM Data , downloadable for free, satisfied the project requirements.

The USGS has DEM data available in multiple "resolutions", "versions" and "levels" for each quadrangle. The resolution specifies the grid cell size, e.g., 10 or 30 meters on a side. The version refers to how recently the DEM data for a particular quadrangle were generated; the higher the version number the more recently the DEM was created. The level refers to the analytic technique used to create the DEM from the published topographic map. The quality increases with the level; however, the availability of the highest level DEMs is limited.

The 30-meter resolution DEM data were used for the Comprehensive Study because they were available for nearly the entire study area. (Of the approximately 3000 quadrangle datasets downloaded, only seven of the 30-meter 7.5-minute DEMs were unavailable. The data used to fill in these areas is described in Special Processing, number 3.)

 

Download DEMs

Each of the 30-meter and 10-meter DEMs corresponds to a standard USGS 7.5-minute quadrangle, and comes in the form of a Spatial Data Transfer Standard (SDTS) package. Each package consists of a compressed tar archive containing various files with consistent names and contents.

The USGS posted 30-meter DEMs for almost all of the country for free download. Unfortunately, its interface provides no convenient way to batch download multiple DEMs, based on drainage basin, state, county, et cetera. You can download them interactively (one at a time) from the USGS Earth Resources Observation Systems (EROS) Data Center (EDC). You can click through online index maps, or if you already know which quadrangles you need, select them from lists arranged by state or quadrangle name.

The step-by-step procedure for downloading multiple 7.5-minute DEMs is outlined in Appendix A, Downloading DEM Data.

[Back to Table of Contents]

 

Source DEM Characteristics

Despite appearing as a homogeneous data set coming from the same agency, the 2905 DEMs used in the study displayed a number of inconsistencies. The USGS derived the DEM data from its published 7.5-minute topographical maps, and carried many of the characteristics from the source quad through to the DEM product.

All the data came projected in Universal Transverse Mercator (UTM), but for the project study area about half fell in UTM zone 10 and half in UTM zone 11. The majority of the data used the North American Datum of 1927 (NAD27) for the horizontal datum. However, ten of the DEMs, all in UTM zone 10, had a horizontal datum of NAD83. All used metric units for horizontal distances.

Inconsistencies in the vertical datum and units proved more vexing. About 85 percent of the data came in meters, with the rest in feet. Whether in feet or meters, USGS rounds the values to the nearest whole number. This often caused trouble along the edge lines when combining adjacent DEMs with different units. Converting the units of one DEM to that of the other created artificial discontinuities along the seam. This only caused problems in very flat areas that straddle quad boundaries, such as lake surfaces. However, representing one half of a lake even a little higher than the other half results in erroneous drainage analyses. Figure 2 helps visualize the situation. In this case, three 7.5-minute DEMs with elevations in meters, and one with elevations in feet, intersect in the middle of Folsom Lake. The elevation data show the water surface as 142 meters throughout most of the lake, but 466 feet (142.04 meters) in the southwest quadrant where the outlet is located. Combined with inadequate representation of a side dam, this four centimeter difference causes Folsom Lake to drain the wrong way. Approximately a dozen similar situations within the project study area required manual editing to resolve.

 Figure 2. Folsom Lake, California

Most of the 7.5-minute DEMs used NGVD (MSL 1927) for the elevation reference, while 70 listed the vertical datum as Local Mean Sea Level (LMSL). These 70 were located predominantly in the Sacramento River basin, but showed no apparent pattern to their distribution. In order to merge the DEMs, the grid elevations required adjustment to a common datum. Initially, the solution seemed simple, because whether the DEM listed its elevations as NGVD or LMSL, the USGS always included a conversion to NAVD (MSL 1988). For reasons described in Target DEM Characteristics, NAVD proved less desirable than NGVD as the common vertical datum, and the final DEM used NGVD.

The analytic technique for creating DEMs continues to evolve while the USGS assembles and updates its 30-meter coverage, resulting in some quality variations within the dataset. About 97 percent of the quadrangles in the study area featured a level 2 DEM, while the remaining data were only level 1. The manually profiled (level 1) DEMs tend to exhibit a "washboard" effect, so that the data seem to show a series of wavy horizontal ridges running through it. Such horizontal striping may result in poorly defined drainage feature outlines, drainage paths biased in the east-west direction, and drainage blockages in the north-south flow component (Garbrecht and Martz, 1999). Techniques exist for "smoothing" the DEM to minimize the effect (Russell, Kumler, and Ochis, 1995), but were not applied in this study, since the drainage analysis seemed to yield acceptable results without such intervention.

Of the 2905 DEMs used 2852 were version one, 52 were version two, while one was version three.

Age of the source topographic maps varied considerably among the DEMs: seven percent dated from the 1990s, 27 percent from the 1980s, 16 percent from the 1970s, ten percent from the 1960s, nine percent from the 1950s, and three percent from the 1940s. The remaining 28 percent failed to show the date of the source map. The oldest and newest dates listed were 1942 and 1998. It seems reasonable to expect that the DEMs included revision data (i.e., the purple features on many topographic maps), but the datasets contained no confirmation.

Figure 3 demonstrates that even if the DEM data included the latest revisions to the source topographic map, it may still leave much to be desired. Built in 1979, the New Melones Reservoir straddles the boundaries of four quads. DEMs for the two southern quads both list dates of 1987, and show the water surface as 332 meters NGVD (dark blue on the figure). The northeast quad fails to report a date; it shows a water surface of 319 meters NGVD (light blue). The northwest quad lists the source topographic map date as 1962, and fails to show the reservoir at all. The latest photo-revision of this map occurred in 1973.

Figure 3. Effect of Source Map Age Differences on New Melones Reservoir

 

[Back to Table of Contents]

 

Target DEM Characteristics

Projection

Inconsistencies in projection and horizontal datum required converting the digital elevation data to a common set of these characteristics before merging them together. The UTM zone 10/NAD27 data formed one intermediate digital elevation model; the UTM zone 11/NAD27 data made up a second DEM, while the ten UTM zone 10/NAD83 DEMs assembled into a small third DEM. The three intermediate DEMs were converted to a common projection and combined into a single digital elevation model for the project study area.

Due to the extent of the study area, the projection used for the master DEM needed to cover almost the entire state with minimal distortion. A projection advocated by California's Teale Data Center fit the requirements. It also turned out to be convenient for working with other GIS coverages. Digital Raster Graphics (DRGs) freely available online use the Teale Albers projection. [Teale Albers projection parameters]

Vertical Datum

Standardizing on a vertical datum and units proved more difficult. Elevations for the initial digital elevation model produced for the study were converted to NAVD, in the (mistaken) belief that this was necessary to provide a common datum for 7.5-minute DEM data using NGVD or LMSL. This first master DEM covered the entire state, taking up about one gigabyte of disk space. However, the initial DEM turned out to be unsatisfactory for the study for a pair of reasons.

First, the datum conversion technique introduced artificial discontinuities in the elevations along the seams between adjacent quadrangles. The vertical datum conversion varied with location, adding for instance 0.94 meters to one DEM, while adding 0.96 meters to its neighbor, creating an effect similar to the feet versus meters issue described in Source DEM Characteristics. Although these differences amounted to only a few centimeters, they still presented difficulties in the drainage analysis. Since the datum adjustment varied for practically every quadrangle, the conversion to NAVD created drainage analysis problems on every water body that straddled a quadrangle boundary. Second, the Sacramento District Corps of Engineers expressed a preference for elevations in NGVD.

For these reasons, the DEM data were reassembled without converting to NAVD. Also, instead of creating one master digital elevation model covering the entire state, two DEMs were made: one for the Sacramento River basin and one for the San Joaquin River basin and Tulare Lake region. This kept the master DEMs to a more manageable size, focusing only on the areas required for the study.

Inquiries to the USGS revealed that LMSL should only be found on islands and in the San Francisco Bay and delta regions. Some (or all) of the 70 DEMs labeled LMSL were apparently mislabeled, but no one could definitively explain the situation. The conversion factor supplied with the DEMs in LMSL turned out to be the conversion from NGVD to NAVD. Therefore, HEC staff decided to treat the LMSL DEMs as if they were actually in NGVD, and visually check them for serious discontinuities. The approach worked well in the upland areas processed with GeoHMS; however, the DEMs covering the Central Valley area, a flat, low-relief landscape, were not scrutinized closely enough to draw conclusions regarding the assumptions.

[Back to Table of Contents]

 

Intermediate Grid Processing and Assembly

To make the two master digital elevation models used for the study, 1063 DEMs were mosaicked for the Sacramento River basin, and 1405 were mosiacked for the San Joaquin basin. The Arc/Info GRID function "mosaic" creates one grid out of two or more adjacent grids and is designed primarily for continuous data; however, only 49 grids can be mosaicked at a time. Consequently, combining the DEMs required an intermediate stage, after which the intermediate mosaics were mosaicked. Each of the raw 7.5-minute DEM data files for UTM zone 11 were grouped into directories based on latitude. The same organization held true for the UTM zone 10 files, except for separating the handful of NAD83 DEMs into their own directory.

The Arc/Info macro sdtsdem.aml (Appendix B - sdtsdem.aml) and associated shell script sdtsinfo.s create an intermediate DEM from raw 7.5-minute DEM data in specified directories. The macro automatically de-archives each DEM package into a temporary directory and runs the script, which returns key properties of the DEM to the macro. The macro imports the SDTS files into an Arc/Info grid, and verifies that all the DEMs share the same horizontal datum, units, and projection. Then it corrects cells with the USGS null value to the Arc/Info "nodata" value. It also replaces obviously bogus elevations, i.e., 1/10 of a mile underwater or higher than Mt. Everest, with the nodata value.

If applicable, the macro next applies the conversion to NAVD that comes from the DEMs' metadata. (As stated previously, the conversion to NAVD was not performed for the final master grids in the Comprehensive Study.) The macro then converts the vertical units to decimeters. Aside from expressing any English system elevations in metric, this also produces an integer grid, without significantly rounding the DEM data. Integer grids required one-quarter to one-half the disk space of comparable floating-point grids. Discrete values also provide advantages when inspecting and editing in the final stages of grid preparation.

Decimeters were used for the vertical units instead of centimeters for two reasons. Most significant is a bug in Spatial Analyst that displays all integer grid values above a certain value (near 100,000) as a single color. Using decimeters kept the values in the tens of thousands, and avoided this bug. Also, using decimeters cut down the number of artificial discontinuities on water surfaces.

Finally, the macro mosaics up to 49 individual grids into an intermediate grid and displays it for inspection. If any of the component 7.5-minute DEMs contain visible errors, its raw DEM data can be moved out of the directory for that group and processed specially.

[Back to Table of Contents]

 

Special Processing

Considering the amount of data processed, the 30-meter 7.5-minute DEMs contained very few errors or omissions that might have adversely affected the drainage analysis in the upland areas of the study. However, the sheer scope of the coverage and the relative immaturity of the dataset practically guaranteed some exceptions to the processing steps described above. The following situations required alternative processing methods to generate intermediate grids from the raw data.

  1. The elevations for one 7.5-minute DEM (Caldwell Butte) that listed vertical units as meters were actually in feet. A copy of sdtsdem.aml, with the conversion factor inverted, processed this DEM into its own intermediate grid.
  2. Two DEMs that listed vertical units as meters were victims of some unit conversion mistake that left them at 30.48 percent of their proper values (Horse Creek, Duzel Rock). A copy of sdtsdem.aml, with the conversion factor adjusted to compensate for the error, processed the DEMs into an intermediate grid.
  3. No 30-meter DEMs were found for seven of the quadrangles, although 10-meter 7.5-minute DEMs were available for each. Sdtsdem.aml was used to produce a pair of 10-meter intermediate DEMs, one for the quads in UTM Zone 10 and the other for UTM Zone 11. These were resampled into two 30-meter intermediate grids using the Arc/Info GRID command "blockmean".
  4. For two quads, one each in UTM zones 10 and 11, no DEMs of any resolution were available for free downloading. The traditional, "native format" DEM data filled in these holes. A couple of intermediate grids were created using the Arc/Info GRID command "demgrid". Unit conversions, bogus data filtering, datum adjustments, and integer rounding operations were performed by manually entering the various Arc/Info GRID commands.

[Back to Table of Contents]

 

Final Assembly

At this stage, the intermediate DEMs all shared the same grid-cell spacing, vertical datum (NGVD), and horizontal and vertical units, although they differed slightly in UTM zone and horizontal datum. The next step mosaicked together all the intermediate DEMs with the same UTM zone and horizontal datum.

Mosaicking the intermediate grids resulted in three DEMs covering nearly the entire state of California: two massive grids for UTM Zone 10/NAD27 and UTM Zone 11/NAD27, and a small one for UTM Zone 10/NAD83. These DEMs were trimmed to include only the Sacramento River and San Joaquin River/Tulare Lake basins and eliminate unneeded grid cells outside the study area. The trimming was performed along the major HUC boundaries with a five-kilometer buffer around the perimeter.

Each of these three grids was then projected to the Teale Albers projection. This very computer-intensive process took two to three days for the biggest grids. Once the three pieces were projected into a common coordinate system, they were mosaicked together to form the two master grids used in the Comprehensive Study: one for the Sacramento River basin and one for the San Joaquin River/Tulare Lake basin. The final step used the Arc/Info GRID command "focalmean" to estimate values for any gaps between grids, based on the elevations in neighboring cells.

As noted in Target DEM Characteristics, Vertical Datum, initially a DEM covering the entire state was created as a single grid in elevation NAVD. For reasons pointed out previously, the final master grids were left in elevation NGVD.

[Back to Table of Contents]

 

DEM Inspection and Correction

The first step in GeoHMS performs a drainage analysis over the entire digital elevation model. Once GeoHMS establishes the drainage characteristics of the basin, it allows tremendous flexibility in creating and re-working the delineation of sub-basins. The current version makes no provision, however, for applying corrections in the original digital elevation model to all the various derived data layers. If defects in the DEM data come to light during the GeoHMS process, the modeler must either account for the inaccuracies later in HEC-HMS, or fix the initial grid external to GeoHMS and repeat all the GeoHMS processing.

Consequently, HEC staff spent considerable time manually inspecting, editing, and testing the master DEMs before using them in GeoHMS. Although Arcview's Spatial Analyst program proved useful for viewing the DEMs, commands from the Arc/Info GRID program seemed best suited for most of the work. Despite their sometimes-arcane nature, the GRID commands lent themselves to a very consistent ("cookbook") approach. Once the appropriate methodology had been worked out, students and engineers with little GIS background accomplished most of the DEM editing.

Appendix C contains some of the Arc/Info GRID commands used in the inspection and editing process.

 

Examples of DEM Errors

Figures 4 and 5 illustrate specific examples of inadequately or inappropriately represented features.

Figure 4. Flow Accumulation Grid near Huntington Lake, California

The circles in Figure 4 mark the saddle dams on Huntington Lake that required editing so that water flowed out of the correct channel. This situation seemed common in the 30-meter data. The false outlets tended to join the proper channel a short distance downstream, adding extra area to the lake's local drainage. Engineering judgement considered the significance of this additional area compared to the lake's real sub-basin size, as well the reservoir's modeling priority, in determining whether to spend time correcting the DEM data.

 Figure 5. Causeway over Goose Lake

In Figure 5, the Westside Road Causeway across Goose Lake appears as a dam. How much of a problem this actually posed to the drainage analysis remains unknown, since the Goose Lake region was dropped from the study. The example serves only to demonstrate potential pitfalls.

 Figure 6. Central Valley near Sutter Buttes, California

Figure 6 illustrates some of the reasons that HEC did not develop GeoHMS models based on the 30-meter DEM data in the flat, low-relief Central Valley region of the project area. Although initially difficult to look at, close inspection of the figure reveals several potential problems for a drainage analysis. It covers a couple of dozen quads centered on Sutter Buttes, California. The area shares many characteristics with the rest of the Central Valley grid. These characteristics include a mixture of English and metric elevation units; both level 1 and level 2 DEMs; extensive areas of very low relief; a 40-year range in the age of the source maps; and a very heavily-engineered system of canals, levees, and other works.

The DEMs tended to blend well in most areas, making it difficult to pick out the boundaries of individual quads. In the flat Central Valley, however, the lack of relief magnifies the effects of minor discontinuities and causes the borders to stand out. The figure shows each unique elevation as one of sixteen different colors. Several quads contain large areas of a single elevation. Even the small ridges introduced by disparate vertical units could play out into square miles of misdirected drainage. Additionally, for the level 2 DEMs, the root mean square error (RMSE) of the DEM elevations may be up to half a contour interval on the published quad map (USGS Mapping Applications Center, 1998). The flat terrain substantially amplifies this error in hydrologic parameters based on slope.

The DEMs in this vicinity suffer from more problems than mixing elevations in feet and meters. The lateral 'striping' typical of level 1 data in flat areas appears near the top of the figure. Drainage patterns based on these DEMs would likely contain substantial error. In the DEMs immediately below Sutter Buttes, the contours break abruptly along quad boundaries, or fail to match up at all. The Sutter Bypass, seen clearly in two quads along the lower right edge, disappears completely in the quad closest to the Buttes. A bust in the vertical datums possibly explains some of the edge discontinuities around this quad, which claimed to use LMSL.

[Back to Table of Contents]

 

National Elevation Dataset (NED)

 

The USGS has recently made available the National Elevation Dataset (NED). It is intended to be a single, homogenous, seamless, continuously updated source of the best topographical data available for the United States. The data carries a modest price. Customers indicate their region of interest by specifying a bounding-box, choose among several formats, and get the data shipped to them on CD-ROM some weeks later.

The NED shows great promise for simplifying the digital terrain preparation. The GIS-intensive work of importing DEMs, normalizing vertical units and datums, mosaicking data together, and filling in holes has already been performed by the USGS. HEC staff recently purchased and began evaluating NED data for some the problem areas described in this paper (Figures 3-6). The preliminary results are mixed:

 

 

 

[Back to Table of Contents]

 

Conclusions

Nearly 3000 USGS 7.5-minute DEMs were downloaded, imported to Arc/Info, mosaicked and edited for use in HEC-GeoHMS. Because of time constraints, the editing process was limited to approximately one to two weeks. Editing was focused on water bodies, especially those water bodies located at quadrangle boundaries.

Subsequent experience revealed that keying in on discontinuities in water bodies along quad boundaries failed to catch many of the inadequacies of the grid. The large, obvious discrepancies seemed to cause little difficulty by themselves, although they tended to exacerbate problems due to defects in the DEMs. The inspection, editing, and correction that takes place prior to working with DEM data in future GeoHMS activities should focus more on ensuring appropriate representation of man-made features, such as dams and levees. Grids for reservoirs and areas with flood control works should be tested with drainage analyses prior to running GeoHMS. The 30-meter data lacked the resolution to adequately define dikes, saddle dams, small levees, and the like. Perhaps the 10-meter DEMs minimize such issues, but modelers should exercise due caution until the datasets earn our faith.

The immediately visible problems in the flat Central Valley region discouraged HEC staff from modeling the area with GeoHMS. Consequently, several concerns regarding the DEMs suitability in these areas remain untested. The adequacy of the representation of the complex flood control structures in the valley remains in doubt. Questions still exist about the effects of the "fill" algorithm across uneven, very flat areas. Concerns over the age of source maps in areas of rapid land subsidence have not been addressed.

The NED substantially reduces the GIS skill necessary to assemble a terrain model, but does not seem to reduce the need for manual corrections. With regard to conducting drainage analyses, it does not appear to address the inadequacies of the source data, and may introduce a few new difficulties.

[Back to Table of Contents]

 

Acknowledgements

 

The U.S. Geological Survey provided all of the DEM data used in the Comprehensive Study. The majority of the data were obtained free of charge from the USGS Earth Resources Observation Systems (EROS) Data Center (EDC). Messrs. George Miller and Mike Price of the USGS Rolla office and Rocky Mountain Mapping Center, respectively, provided assistance in obtaining and interpreting data.

[Back to Table of Contents]

 

Appendices

 

Appendix A - Downloading DEM Data

Appendix B - sdtsdem.aml

Appendix C - Useful Arc/Info GRID Commands

[Back to Table of Contents]

 

Appendix A - Downloading DEM Data


The USGS posted 30-meter DEMs for almost all of the country for free download. Unfortunately, its interface provides no convenient way to batch download multiple DEMs based on drainage basin, state, county, et cetera. You can download them interactively (one at a time) from the USGS EROS Data Center at

http://edcwww.cr.usgs.gov/Webglis/glisbin/searchchoiceftp.pl?dataset_name=7_MIN_DEM.

Each of the 30-meter and 10-meter DEMs corresponds to a standard USGS 7.5-minute quadrangle and comes in the form of a Spatial Data Transfer Standard (SDTS) package. Each package consists of a compressed tar archive containing various files with consistent names and contents.

You can click through online index maps, or if you already know which quadrangles you need, select them from lists arranged by state or quadrangle name.

If you need help determining which quadrangles cover your study area, anonymously ftp the file "selquads.exe" from hec63.hec.usace.army.mil/pub. Run this self-extracting archive from an empty directory and follow the instructions in the readme.txt file. It walks through the process for selecting quadrangles based on Hydrologic Unit Codes (HUCs). It also describes a procedure for automating batch downloads of DEMs using ftp scripts.

As described in Data Acquisition, the USGS has DEM data available in multiple "resolutions", "versions" and "levels" for each quadrangle. The dataset filenames reflect these characteristics. You must choose the most appropriate for your needs. Filename 30.2.1.1008352.tar.gz indicates that the resolution is 30 meters, the level is 2, and the version is 1. The "1008352" part seems to be only a unique identifier; the only pattern detected is that the values are closest among DEMs with similar creation dates.

If DEM datasets for a quadrangle are available in both 30-meter and 10-meter resolutions, choose the one appropriate to the grid you are making. In some cases the 10-meter data may be superior to the 30-meter data in more than just finer resolution, i.e., newer source data, or better production technique. Therefore, even if you want a 30-meter grid, you may want to resample some of it from 10-meter DEMs. Unfortunately, the USGS does not describe the characteristics of each DEM online, so you have to download both versions to see the differences.

Level 2 procedures yield grids much better for hydrologic analysis than level 1, which tends to induce a horizontal "striping" effect in flat areas. Always choose level 2 data over level 1. The highest version number indicates the most recently generated DEM. For a given level, always choose the most recent version.

At first, put the downloaded SDTS packages all into a single directory. The directory name can be anything you like. For the example contained in Appendix B, the datasets will be located in the directory /usr4/shenando/download.

[Back to Table of Contents]

 

Appendix B - sdtsdem.aml


 

Appendix B Steps:

  1. Get sdtsdem.aml and Associated Files
  2. Compile DEM Data Attribute List
  3. Organize DEM Data
  4. Import DEM Data and Create Intermediate Grids

[Back to Table of Contents]

 

 

Once the SDTS packages have been downloaded, the DEM data need to be imported, processed and prepared for use. The importing and intermediate processing procedures contained herein use Arc/Info version 7.2.1.

Determine whether Esri's SDTSImport facility is installed by typing "sdtsimport" at the "Arc:" prompt. If it displays a usage message, then you have it installed. If the command goes unrecognized, download it from http://www.Esri.com/software/arcinfo/sdts.html. This command is reportedly built into Arc/Info version 8.0.

 

Get sdtsdem.aml and Associated Files

Anonymously ftp the file "sdtsdem.tar.Z" from hec63.hec.usace.army.mil/pub into a separate directory. Again, the directory name is up to you. This description will use the directory containing the DEM files (/usr4/shenando/download).

Now extract the processing programs into the parent of the download directory:

cd /usr4/shenando

compress -cd download/sdtsdem.tar.Z | tar -xvf -

Six files should result:

  1. "gzip" is an archiving utility needed to decompress the SDTS packages ("GNU Zip"). It is very handy even outside of this context, and you likely already have it installed. If not, you should ask your system administrator to move it into an appropriate directory in your $PATH, such as /usr/local/bin.
  2. "makeindex.s" is a Unix shell script that combs through the SDTS packages and generates a text file containing important characteristics of each DEM.
  3. "sdtsdem.aml" is the Arc/Info macro that processes all the raw DEM packages in a given directory, and mosaics them into a single grid.
  4. "sdtsinfo.s" is another Unix shell script called internally by "sdtsdem.aml".
  5. "msconfirm.aml" and "msconfirm.menu" are Arctool files that are needed by "sdtsdem.aml".

[Back to Appendix B Steps]

 

Compile DEM Data Attribute List

From the directory containing the DEMs, run the script "makeindex.s" to compile a list of the attributes of the data. The commands would be:

cd /usr4/shenando/download

sh ../makeindex.s

 

This creates a file named "demindex.txt" that lists the quadrangle, downloaded filename, date of source quadrangle, latitude and longitude of the southeast corner, projection, horizontal datum, vertical datum, vertical units, and conversion to NAVD (MSL 1988). Lines from an example demindex.txt file are shown below.

FRONT ROYAL, VA - 24000 30.1.1.1116391.tar.gz Unsp LAT:: 38.875 LONG:: -78.125 UTMZ17 NAS NGVD meters -0.190000
MIDDLETOWN, VA - 24000 30.1.1.1116547.tar.gz Unsp LAT:: 39 LONG:: -78.25 UTMZ17 NAS NGVD meters -0.170000
MIDDLEWAY, WV - 24000 30.1.1.1124891.tar.gz Unsp LAT:: 39.25 LONG:: -77.875 UTMZ18 NAS NGVD meters -0.170000
HARPERS FERRY, WV-VA-MD - 24000 30.1.1.1176147.tar.gz 1986 LAT:: 39.25 LONG:: -77.625 UTMZ18 NAS NGVD meters -0.200000
BLUEMONT, VA - 24000 30.2.1.1115558.tar.gz 1970 LAT:: 39 LONG:: -77.75 UTMZ18 NAS NGVD meters -0.210000
CHESTER GAP, VA - 24000 30.2.1.1116393.tar.gz 1967 LAT:: 38.75 LONG:: -78.125 UTMZ17 NAS NGVD feet -0.610000
ASHBY GAP, VA - 24000 30.2.1.1118389.tar.gz 1970 LAT:: 39 LONG:: -77.875 UTMZ18 NAS NGVD meters -0.200000
STEPHENSON, VA-WV - 24000 30.2.1.1118508.tar.gz Unsp LAT:: 39.125 LONG:: -78 UTMZ17 NAS NGVD feet -0.580000
FLINT HILL, VA - 24000 30.2.1.1124707.tar.gz 1966 LAT:: 38.75 LONG:: -78 UTMZ17 NAS NGVD meters -0.200000
CHARLES TOWN, WV-VA-MD - 24000 30.2.1.1168763.tar.gz 1984 LAT:: 39.25 LONG:: -77.75 UTMZ18 NAS NGVD feet -0.620000
LINDEN, VA - 24000 30.2.1.1179938.tar.gz 1994 LAT:: 38.875 LONG:: -78 UTMZ17 NAS NGVD meters -0.200000
ROUND HILL, VA-WV - 24000 30.2.1.1179952.tar.gz 1998 LAT:: 39.125 LONG:: -77.75 UTMZ18 NAS NGVD meters -0.200000
BERRYVILLE, VA-WV - 24000 30.2.1.1237216.tar.gz Unsp LAT:: 39.125 LONG:: -77.875 UTMZ18 NAS NGVD feet -0.630000
UPPERVILLE, VA - 24000 30.2.1.1237218.tar.gz Unsp LAT:: 38.875 LONG:: -77.875 UTMZ18 NAS NGVD feet -0.670000
BOYCE, VA - 24000 30.2.1.1237230.tar.gz Unsp LAT:: 39 LONG:: -78 UTMZ17 NAS NGVD feet -0.630000
STEPHENS CITY, VA - 24000 30.2.1.1237240.tar.gz Unsp LAT:: 39 LONG:: -78.125 UTMZ17 NAS NGVD feet -0.580000
STRASBURG, VA - 24000 30.2.1.1237262.tar.gz Unsp LAT:: 38.875 LONG:: -78.25 UTMZ17 NAS NGVD meters -0.190000

 

Look through the list for blank or garbled fields. Also look for vertical datums other than NGVD, such as LMSL or NAVD. These DEMs require special handling.

[Back to Appendix B Steps]

 

Organize DEM Data

All the 10-meter and 30-meter DEM data come with horizontal units of meters expressed in the Universal Transverse Mercator (UTM) projection. However, the DEMs for your study area may fall across more than one UTM zone or use differing horizontal datums. Grids must be in the same UTM zone and use the same horizontal datum for them to line up properly when merged. Using the index file just created in Compile DEM Data Attribute List, organize the DEMs into directories for each combination of UTM zone and horizontal datum. The directories can be named whatever you want, although it makes sense to include the zone and datum in the name. For example, put the DEM files for all the quadrangles in UTM zone 10 using NAD27 (i.e., NAS) into a directory named z10nas, and all the files for UTM zone 10 using NAD83 (i.e., NAX) into one named z10nax.

[Link to Universal Transverse Mercator grid map]

Symbolic links work very well for this purpose. Instead of copying a file to the destination directory, the link merely points to the actual location of the file. The link makes the file appear to exist in the second directory, but does not actually double the disk space used. For example, the following command creates a link named boyce_va.tar.gz in the directory z17nas, which references downloaded file 30.2.1.1237320:

ln -s /usr4/shenando/download/30.2.1.1237230.tar.gz

/usr4/shenando/z17nas/boyce_va.tar.gz

If any DEMs use vertical datums other than NGVD, be sure to copy or link them into separate directories for special handling later.

Finally, Arc/Info GRID can only mosaic 49 grids at one time. Any directories with more than 49 DEMs will have to be broken up into two or more directories. If you need to recreate a link in another directory, first remove it from the overcrowded directory using the "rm" command:

rm /usr4/shenando/z17nas/boyce_va.tar.gz

Similarly, separate any DEMs of a resolution different than your primary cell spacing. These must also be separated according to any variations in UTM zone, horizontal datum, and vertical datum. For example, say you are building a 30-meter DEM that overlaps UTM zones 18 and 17, and all the DEMs use NAD27 and NGVD. No 30-meter DEM is available for one of the quadrangles in zone 18, but a 10-meter DEM is available. This situation would require distributing into three directories, with names something like z17nas, z18nas, and z18nas10.

[Back to Appendix B Steps]

 

Import DEM Data and Create Intermediate Grids

Start with a directory containing one or more DEMs with consistent properties. The cell size, projection (UTM zone), horizontal datum, horizontal units, and vertical datum must be identical for the DEMs in each directory created in Organize DEM Data. Differences in vertical units are common and are handled during the import process. The downloaded DEMs in each directory will be mosaicked together into intermediate grids. These intermediate grids will later be converted to a common grid spacing and coordinate system, and mosaicked into a master grid.

Choose a directory for your intermediate grids, copy the files "sdtsdem.aml", "sdtsinfo.s" "msconfirm.aml" and "msconfirm.menu" into it, and run Arc/Info. (See Get sdtsdem.aml and Associated Files for downloading instructions.) This directory may be located and named whatever you like. This example just uses the same directory that contains the copied or linked DEMs:

cd /usr4/shenando/z17nas

cp ../sdts* .

cp ../msconfirm* .

arc

From the Arc: prompt, run "sdtsdem.aml" with two arguments. The first argument should be the full pathname of one of the directories containing DEMs. The second argument should be the name of the intermediate grid you are creating. For example:

&r sdtsdem /usr4/shenando/z17nas iz17nas

The script will begin processing the DEMs in /usr4/shenando/z17nas, printing various messages along the way. Ignore the lines saying "Warning: new location is not a workspace". The macro automatically de-archives each DEM package into a temporary directory and runs the shell script "sdtsinfo.s", which returns key properties of the DEM to the macro. The macro imports the SDTS files into a temporary Arc/Info grid and verifies that all the DEMs share the same horizontal datum, units, and projection. Then it corrects cells with the USGS null value to the Arc/Info "nodata" value. It also replaces obviously bogus elevations, i.e., 1/10 of a mile underwater or higher than Mt. Everest, with the nodata value.

If applicable, the macro next applies the conversion to NAVD that comes from the DEMs metadata. At present, users must edit "sdtsinfo.s" and uncomment the appropriate lines to activate this feature. Then it converts the vertical units to decimeters.

These operations result in a second temporary grid for each DEM.

Finally, the macro mosaics the individual DEMs into an intermediate grid and displays it for inspection. You will likely notice one-cell gaps of NODATA along some of the boundaries between individual DEMs. These gaps can be filled later using the Arc/Info GRID command "focalmean". Systematic errors in any of the component DEMs are usually visible as a quadrangle being shaded differently than its neighbors or misaligned. When such anomalies are encountered, separate the errant DEM into its own directory and run "sdtsdem.aml" to build the intermediate grid without it.

After "sdtsdem.aml" finishes, review the file "sdtsdem.log". In particular, check to see that the minimum and maximum values for each DEM are reasonable. Remember that these are expressed in decimeters.

[Back to Appendix B Steps]

 

Appendix C - Useful Arc/Info GRID Commands


 

The following are some commands available at the GRID prompt that will help you do various useful things. Remember to avoid uppercase letters in anything involving Arc/Info, and to use plenty of spaces in Arc/Info commands.

Display Commands

Grid Management Commands

Grid Editing Commands

Drainage Analysis Commands

Miscellaneous Commands

[Back to Table of Contents]

 

The following commands help you plot grids and coverages in a display window:

Display Commands:

Command:

Description:

display 9999

Brings up an x-window for your graphics.

mape <gridname or coveragename>

mapextent: The extent of the map is the extent of [gridname], and it must be set before drawing. This command is like a zoom level.

gridp <gridname>

gridpaint: Plots [gridname] using a unique-value scheme.

gridp <gridname> value linear nowrap gray

Paints <gridname> using a continuous grayscale.

polys <coverage>

Draws polygons of specified coverage.

point <coverage>

Draws points of specified coverage.

arcs <coverage>

Draws vector data.

linecolor <color>

Changes the color used for plotting vector data.

[Back to Appendix C - Grid Commands]

 

Grid Management Commands:

Command:

Description:

lg <directoryname>

Lists the grids in the specified directory. If lg command given without directory name, the current workspace grid contents are listed.

lc <directoryname>

Lists the coverages in the specified directory. If lc command given without directory name, the current workspace coverages are listed.

describe <gridname>

Gives information about a grid including if it contains floating point or integer values.

copy <in_grid> <out_grid>

Copies a grid (its directory and associated files) to the current workspace. The path of the grid to be copied (in_grid) must be specified. If the grid is not to be copied to the current workspace, a path for out_grid must also be specified. Example: copy sj88 /usr4/sjbasin/sj88.

rename <old gridname> <new gridname>

Replaces old gridname with new gridname.

kill <grid, coverage, etc.> all

Deletes the grid or coverage and all associated files.

[Back to Appendix C - Grid Commands]

 

Grid Editing Commands:

Command:

Description:

gridclip <in_grid> <out_grid> *

Used to work on a smaller area defined by a box drawn on the grid with a mouse. The * signifies that the user will draw a box to define the area to be clipped. Example: gridclip sac2729 berry *

gridclip <in_grid> <out_grid> cover <coverage>

Creates a grid file from another grid file with the boundaries of the coverage file. Example: gridclip sac8329 featherdem cover feather

cellvalue <gridname> *

Used to determine value of cells by placing cursor on grid.

gridedit edit <gridname>

Moves into editing mode for grid file.

gridedit fillvalue <value>

Example: gridedit fillvalue 1330 This sets the value to be filled in the grid to 1330.

gridedit fillregion *

Click on the region within the grid to be filled with previously defined fillvalue. All grid cells within a region (defined by adjacent cells of equal value) are filled with the fillvalue.

gridedit fillpolygon *

Create a polygon on the grid within which the fillvalue will be applied.

gridedit fillcell *

Click on individual grid cells to be given previously defined fillvalue.

gridedit save

Save changes to grid.

gridedit quit

Exits the grid-editing mode.

[Back to Appendix C - Grid Commands]

 

Drainage Analysis Commands:

Command:

Description:

fill <in_grid> <out_grid>

Fills each depression in a grid to the elevation of the lowest overflow point out of the sink. For example, fill berry berryf will create a grid berryf by filling in depressions in the grid berry.

flowdirection

Computes flow direction and creates a grid containing cells with integer values 1, 2, 4, 8, 16, 32, 64 or 128 to represent eight directions where 1 is east, 2 is southeast, 4 is south and so on. Example: berryd = flowdirection (berryf, #, #). # represents default

flowaccumulation

Computes a stream network by accumulating flow in each grid cell. Example: berrya = flowaccumulation (berryd, #)

[Back to Appendix C - Grid Commands]

 

 Miscellaneous Commands:

Command:

Description:

&sys <Unix command>

Runs the given Unix command and shows the output.

arc w <directory>

Changes workspace to the specified directory.

arc shapearc <shapefile> <coverage>

Creates a coverage for use in ARC from a shapefile created in ArcView. Example: arc shapearc febuf feather. Note: this command can be used in ARC directly without preceding shapearc with arc.

arc build <coverage>

Builds the ARC coverage

con, focalmean

The following command checks the grid tempsac for nodata, or null, values. It then averages the surrounding area, specified here as a rectangle of area 2 x 2, and fills in the null cells with that average. Example: sac8329 = con (isnull (tempsac), focalmean (tempsac, rectangle, 2, 2), tempsac)

identify <coverage> <feature class> *

Identifies and gives information about the feature selected. The * signifies that the feature will be selected interactively. Example: identify /usr4/sacsjgis/quads polys. This identifies the quad selected with the mouse.

int

Creates an integer grid. Example: sj88dm = int (sj88cm / 10 + .5) This creates an integer decimeter grid from a centimeter grid and includes rounding.

mosaic

Creates a grid composed of all of the input grids. Example: z10 = mosaic ( z1034_675, ... z1037_5). This mosaics all of the input grids and creates a new grid called z10 that contains all the data from the input grids.

project

Projects a grid to the specified projection. Example: z10prj = project ( sjz10, #, BILINEAR, 30, 0, 0). The user will then be prompted with Project: after which a series of key words are entered, e.g., datum, projection, units, parameters.

merge

Merges a series of input grids based on order of input. The first grid has priority over the second and so on. Example: midsacdem2 = ( churn, jewett, toomes, midsacdem). This merges the example edited grids with midsacdem and overwrites the cells in midsacdem with the corrected, or edited cells, from the first input grids.

[Back to Appendix C - Grid Commands]

 

References

Djokic, D. and Z. Ye, 1999, DEM Preprocessing for Efficient Watershed Delineation, Proceedings of the 19th Esri Users Conference, San Diego, CA.

Doan, J., 2000, Development and Application of HEC-GeoHMS, Proceedings of the 20th Esri Users Conference, San Diego, CA (this issue).

Dunn, C. and J. Doan, 2000, Hydrologic Modeling for the Sacramento and San Joaquin Basins using GeoHMS, Proceedings of the 20th Esri Users Conference, San Diego, CA (this issue).

Garbrecht, J. and L. W. Martz, 1999, Digital Elevation Model Issues in Water Resources Modeling, Proceedings of the 19th Esri Users Conference, San Diego, CA.

Olivera, F. and D. R. Maidment, 1999, GIS Tools for HMS Modeling Support, Proceedings of the 19th Esri Users Conference, San Diego, CA.

USGS Mapping Applications Center, 1998, Fact Sheet 102-96, Reston, VA

Russell, E., Kumler, M., and Ochis, H., 1995, Identifying and Removing Systematic Errors in USGS Digital Elevation Models (DEMs), Poster presented at the Annual GIS in the Rockies Conference, Denver, CO.

[Back to Table of Contents]

 

Author Information

 

Matthew M. McPherson

Research Hydraulic Engineer
US Army Corps of Engineers
Hydrologic Engineering Center
609 Second Street
Davis, CA 95616
Tel: 530/756-1104
Fax: 530/756-8520

mattm@superior.lre.usace.army.mil

 

Heather E. Henneman

Hydraulic Engineer
US Army Corps of Engineers
Chicago District
111 N. Canal Street
Chicago, IL 60606
Tel: 312/353-6400, ext. 3008
Fax: 312/353-2156

Heather.E.Henneman@usace.army.mil