My DEM Processing Flow on Linux

Morgan Brown Research Talks Topo Maps Example Maps 2 DEM Processing Adventures Images Extra... TopoZone - Awesome! Free, seamless, online USGS topos. REI - REI sells topo maps - both USGS 7.5' and Trails Illustrated. USGS Mapping - Home base for USGS map info.

Part I - SDTS to normal digital image

1. Download DEM of choice from USGS. Uncompress and expand the .tar.gz file. (Visit GNU to download GNU tar and gzip.)
2. Convert binary SDTS format file to ASCII DEM format. I currently use sdts2dem, but the website appears to be down. Contact me and I'll try to help you get set up. I have a trivial perl script, sdts2dem.pl, which simplifies the use of sdts2dem.
3. Convert ASCII .dem file to gridded image. Eric Kasten distributes a conversion program, dem2ppm, in his free DEM Tools package. Unfortunately the dem2ppm converter mishandles some points on the edge of the grid. Note: I had to hack the code for use on my machine, which has an 8-bit colormap. The default colormap has 1000 grayscales by default, so I had to bump this down to 256. I have a trivial perl script, getgrid.pl, which parses the output for extreme grid and elevation values, which are then written to a file for future use.
4. I use MATLAB to process the data, but MATLAB cannot read .ppm files. So using the image conversion program of your choice, save the .ppm file as a grayscale .tif file.

Part II - Digital image to 3D maps, using MATLAB

1. The earth is not flat. Therefore, projections of "rectangular" regions of the earth (i.e., 7.5 minute maps) onto a Cartesian grid are normally not rectangular. Furthermore, as mentioned above, the edges of maps created by dem2ppm are usually corrupted. For my non-professional purposes, I want to warp the map to a Cartesian mesh, then simply trim the corrupted data. These tasks are accomplished by the MATLAB script demprep.m.

   To warp the data, I have written a small MATLAB script, affineTransf.m, which prompts the user to pick four control points (corners of original map), then warps the map via a simple Affine transformation to the corners of the map, where the coefficients of the transformation are gotten from the user picks. Mathematically, the process is defined as follows:

   	    Compute INVERSE transform coefficients, i.e.,
   	    an 8-vector p such that
   
   	      x = p(1)x' + p(2)y' + p(3)
   	          ----------------------
   	          p(7)x' + p(8)y' + 1
   
   	      y = p(4)x' + p(5)y' + p(6)
   	          ----------------------
   	          p(7)x' + p(8)y' + 1
   
   	    defines the transformation from known output
   	    to input space -- a backward map.
   	    

   where (x',y') and (x,y) are locations in the input and output coordinate systems, respectively. The script has an option to perform this entire warping procedure automatically, but some users may want details and/or more control.

   The corrupted data trimming is more brute force: prompt the user to pick regions of corrupted data, then extrapolate good data to the edge of the image via a recursive, 2-D running median.

2. To display the DEM in 3-D, I have written a MATLAB script called topo3d.m. The script reads a tif file from disk and renders it in 3-D, with many dislay options. The maps are in real earth coordinates, with no vertical exaggeration. An option exists to print the maps to a .tif file, which can be converted to the format of your choice. Click here to view a sample map, this one being the 7.5 minute quad "Stonyford, California".

   Using a simple raytracing methodology, you may compute a visibilty map from any point in 3-D in the map. More specifically, rays are traced from the starting point across many azimuth and elevation angles. When a ray hits the land's surface, the point is marked as "visible". Here is an example, raytraced from the same Stonyford, CA 7.5 minute quad.

3. Adjacent maps can be spliced together to create a mosaic of 7.5 minute maps. The scripts splicelr.m and spliceud.m do this job for left-right and up-down neighbors, respectively. Click here to see a 4x4 mosaic (30 minute) map of the central section of Yosemite National Park. From north to south, the two deep valleys are Yosemite Valley (Merced R.) and the Grand Canyon of the Tuolumne River.

Department of Geophysics Stanford University

Modified: 11/18/05, 13:52:58 PST , by morgan
Page Maintainer: morgan `AT' sep.stanford.edu

---