Rendering U.S.G.S. DEM Quadrangles

by
Russell Towle
February, 1997

Discussion

The ability to accurately visualize the shape of the earth's surface is only rarely enjoyed; pilots and mountain climbers are favored with good views from many directions and under varied lighting conditions, but for most of us, even the most careful analysis of a topographic contour map, aided by  our imagination, may fail to evoke a very good idea of the lay of the land.

Digital topographic maps remedy this deficiency, but the requisite software is rare, and good renderings of landscapes even rarer.  Mathematica includes a variety of plotting functions which can be applied to such files, and this notebook contains procedures to read and render "digital landscapes."

The routines in this notebook are under active development, and much room for improvement remains.  For instance, a generalized function could be devised, which would read any number of any type of Digital Elevation Model (DEM) files, or a Mathematica package might be written, containing a variety of functions and options.  I welcome suggestions!  See the Credits section for how to contact me.

Some United States Geological Survey (U.S.G.S) DEM files are available on-line (see "DEM Sources on the Internet," below).  Two compressed sample DEM files may be downloaded from MathSource.  This notebook reads and renders, using ListPlot3D and other plotting functions, the 7.5 minute DEMs, which offer the highest resolution of landscape features available.  These files span 7.5 minutes of longitude by 7.5 minutes of latitude.  The other types include 1-degree DEMs, 30 minute DEMs, 15 minute DEMs, and special forms of these for the state of Alaska.  A method of reading and rendering 1-degree DEMs is also included here.

The essential structure of a DEM is a simple list of elevations, arranged in "profiles," or columns; elevations within a profile are ordered from south to north, and the profiles themselves are ordered from east to west.  In a 7.5 minute DEM, elevations are given at intervals of 30 meters within each profile, and the profiles are 30 meters apart.  Elevations are integers, in (usually) meters.

The 7.5 minute DEMs, as with the others, are plain text (ASCII) files.  Although this is not the most compact form in which to store data, it provides the advantage of  portability to many different computer platforms and operating systems.  7.5 minute DEMs are about 1 MB in size, and consist of three record types:

1.  A-records contain header information which includes the name of the quadrangle, and a multitude of mapping parameters, most importantly the Universal Transverse Mercator (UTM) coordinates of the four corners of the quadrangle, ordered clockwise from the southwest corner; the number of profiles; the units of elevation measure, which may be meters or feet; and the minimum and maximum elevations within the quadrangle.  UTM coordinates are pairs of eastings and northings, in meters, the northings measured from the Equator, the eastings, from the central meridian of a UTM "zone," which is assigned a value of 500,000 meters to ensure that eastings are never negative.  UTM zones span six degrees of longitude.

2.  B-records begin 1024 bytes into the file, and always begin on some multiple of 1024 bytes.  Each profile has a B-record.  The records consist of a nine-number header, giving the number of elevations in the profile, the UTM coordinates of the "easting" and "northing" of the southernmost elevation in the profile, and the lowest and highest elevations within the profile.  Following this header are the elevations themselves, from south to north.  These vary in number, but once away from the western or eastern edges of the DEM, hundreds of elevations make up a single profile, and two or three blocks of 1024 bytes are required to store them.  

3.  There is only one C-record per DEM, at the very end of the file, which contains information about the accuracy of the DEM.  C-records are ignored by the routines below.

The 7.5 minute DEMs are somewhat awkward to work with, mainly because their profiles of elevations are of varying lengths, and also begin and end at various northings.  That is, all four edges of the DEM are staggered, in much the same way that a straight line at some slight angle, rendered on a computer display, exhibits "jaggies," or a stair-stepped appearance.  In order to use the elevation array in ListPlot3D, the profiles must be padded with zeros above and below, to make the array a rectangular list of lists.

On the other hand, these staggered edges of the 7.5 minute DEMs fit perfectly against the staggered edges of contiguous quadrangles, so it is possible to read whole blocks of quadrangles, joining their elevation profiles, and again padding the outer edges of the array so that all the composite profiles have the same length.  Code for joining several DEMs in this way will be found below.  However, due to the nature of the UTM projection, DEMs that lie to either side some multiple of 6 degrees of longitude will not fit perfectly together.

The worst problems presented to Mathematica users by the routines below are time and memory.  Using a Macintosh Quadra 700 with a PowerPC upgrade card and 32 MB of RAM allocated to kernel and front end alike, reading a single DEM and assembling the rectangular array of elevations requires about 40 seconds, while rendering in ListPlot3D requires around 550 seconds, or a little over nine minutes.  This is possible with a kernel memory partition of less than 20MB.  There are various work-arounds to reduce memory use:

1.  A subset of the elevation array may be taken; the code for this is commented out, in the sections which read DEMs, below.

2.  The final elevation array may be written to a file, whereupon one may kill and then revive the kernel, read the file into a variable, and then pass it to ListPlot3D.

3.  One may allocate nearly all memory to the kernel, and write the output of ListPlot3D to a (PostScript) file, then quit, allocate nearly all memory to the front end, launch the front end, open the PostScript file, select its cell bracket, and assign to it the cell style "Graphics."  However, a block of four DEMs when rendered to PostScript measures 52MB; a single DEM, rendered to PostScript, measures 13MB.

I had to allocate 60 MB of RAM to the kernel, and write the output to PostScript, in order to open and render a block of four 7.5 minute DEMs.

I welcome any comments, bug reports, or suggestions.  See the Credits section for information about how to contact me.

Read one 7.5' DEM

Example

(*Let Mathematica know where your DEMs are located;
in this case, in the folder "7.5" within the hard disk "Starlight"*)

SetDirectory["Starlight:7.5"]

(*Construct an array of elevations from the DEM*)

ReadDEM["Dutch Flat"]

(*Boxed, with axes, and special light sources*)

ListPlot3D[elevarray,
(*DisplayFunction->(Display["a landscape",#]&),*)
ViewPoint->{0, -1, 1.0},(*from south and above*)
PlotLabel->quad,
Boxed->True,
BoxRatios->Automatic,
AxesLabel->{"west to east","south to north","elevations"},
Axes->True,
Mesh->False,
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},
Background->GrayLevel[0],
LightSources->{{{0,0,10},RGBColor[0.4, 0.4, 0.80]},
{{10,-3,6},RGBColor[1, 0.75, 0.00]}}]

(*Grayscale image*)

ListPlot3D[elevarray,
(*DisplayFunction->(Display["a landscape",#]&),*)
ViewPoint->{0, -1, 0.5},(*from south and above*)
PlotLabel->quad,
Boxed->True,
BoxRatios->Automatic,
AxesLabel->{"west to east","south to north","elevations"},
Axes->True,
Mesh->False,
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},
ColorOutput->GrayLevel,
Background->GrayLevel[1]]

Initialization

(*Cast as a function, with four variables left global:
elevarray, quad (quadrangle name), q (dimensions of array)*)

ReadDEM[file_String]:=
Block[{dem,utm,pronum,south,north,capital,
pedestal,pillar,profilelength,elevflag,k},
(*Create an InputStream object, name it "dem"*)
dem=OpenRead[file];

(*Read A-record, retrieve quad name, etc.*)
SetStreamPosition[dem,0];
tmp=ReadList[dem,Word,3,NullWords->True];
mark=StreamPosition[dem];
SetStreamPosition[dem,0];
Print["Quadrangle Name:  ", 
	quad=FromCharacterCode[ReadList[dem, Byte, mark]] ];
	
SetStreamPosition[dem,145];
Print["DEM Level Code:  ", 
	Read[ dem,Number] ];
	
SetStreamPosition[dem,151];
Print["Elevation Pattern Code:  ", 
	Read[ dem,Number] ];			

SetStreamPosition[dem,157];
Print["Planimetric Ground Reference:  ", 
	Read[ dem,Number] ];

SetStreamPosition[dem,529];
u=Read[ dem,Number];
groundunits=Which[u==0,"radians",u==1,"feet",
	u==2,"meters",u==4,"arc-seconds"];
Print["Ground Planimetric Units:  ", groundunits];

SetStreamPosition[dem,535];
elevflag=Read[ dem,Number];
elevunits=Which[elevflag==1,"feet", elevflag==2,"meters"];
Print["Elevation Units:  ", elevunits];	

SetStreamPosition[dem,547];
Print["UTM Corners:  ", 
	utm=Partition[Round[ReadList[ dem,Number,8]], 2] ];

SetStreamPosition[dem,739];
Print["Minimum and maximum elevations:  ",
	Round[ReadList[ dem,Number,2] ]];

SetStreamPosition[dem,858];
Print["Number of elevation profiles:  ", 
	pronum=Read[ dem,Number] ];

(*Display outline of quadrangle*)
AppendTo[ utm, utm[[1]] ];
Show[Graphics[Line[utm] ],
	PlotLabel->quad,
	AspectRatio->Automatic];
 
south=Min[Map[Take[#, -1]&, utm]];
north=Max[Map[Take[#, -1]&, utm]];
 
(*Read B-records, create elevation array*)
SetStreamPosition[dem,1024];
k=1;elevarray={};

While[k<=pronum,(*while k is less than the number of profiles*)

b=Round[ReadList[dem,Number,9]];(*get the header*)

(*b[[3]] is number of elevations in profile*)
(*b[[6]] is the northing of southernmost elevation in profile*)
c=Round[ReadList[dem,Number,b[[3]] ] ];(*get the elevations in this profile*)
capital=Round[(north-(b[[6]]+(30*(b[[3]]-1))))/30];(*how many zeros above*)
pedestal=Round[((b[[6]])-south)/30];(*how many zeros below*)
(*join zeros below, elevations, zeros above*)
pillar=Flatten[{Table[0,{pedestal}],c,Table[0,{capital}]}];
elevarray={elevarray, pillar};k++];


profilelength=Length[pillar];(*all profiles now have equal length*)
elevarray=Flatten[elevarray];(*discard all but one set of braces*)
elevarray=Partition[elevarray,profilelength];(*restore profiles*)

(*Close InputStream*)
Close[dem];

(*If given to ListPlot3D without transposing, landscape is mirror image*)
elevarray=Transpose[elevarray];

q=Dimensions[elevarray];(*needed for MeshRange spec*)
Print["The dimensions of the final elevation array are ",q];

(*Finally, if the measure of elevations is feet, rather than meters,
we must multiply the entire array by .3048*)
If[elevflag==1,elevarray=elevarray*.3048,elevarray];
]

Read multiple 7.5' DEMs

Foreword:  The example below, which joins two 7.5 minute DEMs, required 179 seconds on a Macintosh Quadra 700 with PowerPC upgrade card, and consumed in the neighborhood of 12 MB kernel memory.  The dimensions of the final elevation array are {483,784}, and the array would measure 2MB if written to disk.  The rendering of this example would require about fifteen minutes, using ListPlot3D.

Example

(*This is one way to conserve memory use by the kernel*)
Needs["Utilities`MemoryConserve`"]

(*let Mathematica know where your DEMs are located,
in this case, in the folder "7.5" within the hard disk "Starlight"*)
SetDirectory["Starlight:7.5"]

(*Two DEMs, 174 seconds, Four DEMs, 446 seconds*)
ReadMultipleDEMs[{"Dutch Flat","Westville"}]

(*Default lighting, no box, no axes.  To write
PostScript graphics to disk rather than displaying,
de-comment the second line.  Other rendering schemes
will be found in the "Rendering Variants" section.*)


ListPlot3D[elevarray,
(*DisplayFunction->(Display["a landscape",#]&),*)
ViewPoint->{-0.0, -1.00, 1.0},(*from south and above.*)
PlotLabel->quad,
Boxed->False,
BoxRatios->Automatic,
Axes->False,
Mesh->False,
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},
Background->GrayLevel[0]]

Initialization

ReadMultipleDEMs[dora_List]:=
Block[{dem,mark,u,utm,groundunits,elevflag,
m,minmax,p,b,brec,allbrec,pronum,
easting,longprofile,pronorth,spots,a,
long,south,north,capital,pedestal,pillar,k},

k=1;dem={};
While[k<=Length[dora],
AppendTo[dem, OpenRead[ dora[[k]] ]];
k+=1];

k=1;utm={};quad={};pronum={};elevunits={};minmax={};
While[k<=Length[dem],

SetStreamPosition[dem[[k]],0];
tmp=ReadList[dem[[k]],Word,3,NullWords->True];
mark=StreamPosition[dem[[k]]];
SetStreamPosition[dem[[k]],0];
qind=FromCharacterCode[ReadList[dem[[k]], Byte, mark]];
Print["Quadrangle Name:  ",qind]; 
AppendTo[quad, qind];

SetStreamPosition[dem[[k]],145];
Print["DEM Level Code:  ", 
	Read[ dem[[k]],Number] ];
	
SetStreamPosition[dem[[k]],151];
Print["Elevation Pattern Code:  ", 
	Read[ dem[[k]],Number] ];			

SetStreamPosition[dem[[k]],157];
Print["Planimetric Ground Reference:  ", 
	Read[ dem[[k]],Number] ];	

SetStreamPosition[dem[[k]],529];
u=Read[ dem[[k]],Number];
groundunits=Which[u==0,"radians",u==1,"feet",
	u==2,"meters",u==4,"arc-seconds"];
Print["Ground Planimetric Units:  ", groundunits];

SetStreamPosition[dem[[k]],535];
elevflag=Read[ dem[[k]],Number];
either=Which[elevflag==1,"feet", elevflag==2,"meters"];
Print["Elevation Units:  ", either];
AppendTo[elevunits,elevflag];

SetStreamPosition[dem[[k]],547];
u=Partition[Round[ReadList[dem[[k]],Number,8]],2];
Print["UTM Corners:  ",u]; 
AppendTo[utm, u];

SetStreamPosition[dem[[k]],739];
Print["Minimum and maximum elevations:  ",
	m=Round[ReadList[ dem[[k]],Number,2] ]];
AppendTo[minmax, m];

SetStreamPosition[dem[[k]],858];
p=Read[ dem[[k]],Number];
Print["Number of elevation profiles:  ",p]; 
AppendTo[pronum, p];

Print[""];
k+=1];


(*Display outlines of quadrangles*)
Do[AppendTo[ utm[[i]],utm[[i,1]] ],{i, Length[utm]}];
Show[Graphics[Table[Line[utm[[k]]],{k,Length[utm]}] ],
	PlotLabel->"Quadrangle Outlines",
	AspectRatio->Automatic];
	
	
(*make lists of B-record headers, and elevation profiles*)
allbrec={};profiles={};k=1;
While[k<=Length[dem],
SetStreamPosition[dem[[k]],1024];

	j=1;brec={};elevs={};
	While[j<pronum[[k]],
	b=Round[ReadList[dem[[k]],Number,9]];
	AppendTo[brec,b];
	e=Round[ReadList[dem[[k]],Number,b[[3]] ]];
	If[elevunits[[k]]==1,e*=.3048,e];
	AppendTo[elevs,e];
	j=b[[2]] ];

AppendTo[allbrec,brec];
AppendTo[profiles,elevs];
k+=1];

k=1;
While[k<=Length[dem],
Close[dem[[k]]];
k+=1];

(*make a list of profile eastings, more than one DEM can have one easting*)
easting=Union[
	Flatten[
		Table[allbrec[[k,i,5]],{k,Length[dem]},{i,pronum[[k]]}],
	1]
];
Print["The total number of profiles is ",Length[easting]];
(*join individual profiles to make complete profiles*)
k=1;longprofile={};pronorth={};
While[k<=Length[easting],

	spots=Position[allbrec,easting[[k]]];
	a=Sort[Table[{allbrec[[spots[[i,1]],spots[[i,2]],6]],spots[[i,1]],spots[[i,2]]},
		{i,Length[spots]}]];
	AppendTo[pronorth, Min[Map[Take[#, 1]&, a]] ];
	j=1;
	
	While[j<=Length[spots],
	long=Flatten[Table[profiles[[a[[i,2]],a[[i,3]] ]],{i,Length[spots]}]];
	j+=1
	];
	
AppendTo[longprofile,long];(*add to list of combined profiles*)	
k+=1];

(*Discover maximum and minimum northings of profiles*)
south=Min[pronorth];
north=Max[Table[
		pronorth[[i]]+(30*(Length[longprofile[[i]]]-1)),
		{i,Length[longprofile]}]];
		
		
(*pad profiles with zeros, north and south*)
k=1;elevarray={};
While[k<=Length[easting],
capital=Round[(north-(pronorth[[k]]+(30*(Length[longprofile[[k]]]-1))))/30];
pedestal=Round[(pronorth[[k]]-south)/30];
pillar=Join[Join[Table[0,{pedestal}],longprofile[[k]]],Table[0,{capital}]];
AppendTo[elevarray,pillar];
k+=1];

(*Transpose padded profiles, so that in ListPlot3D x axis is west-east*)
elevarray=Transpose[elevarray];
(*elevarray=Table[Take[elevarray[[i]],{20,350}],{i,20,440}];*)
q=Dimensions[elevarray];
maxmin={Max[Flatten[minmax]],Min[Flatten[minmax]]};
On[General::spell1];
Print["The dimensions of the final elevation array are ",q];
Print["The maximum and minimum elevations in the array are",maxmin]
]

Rendering Variants

Overall note:  if one has written the elevation array to a file, and then quit, re-launched, and read that array back in from the file, the PlotLabel directive will not give the quadrangle names, and you may wish to comment it out.  It should be noted that, when axes are enabled within the plots, all dimensions are given in meters.  If one wishes to write the graphics to a PostScript file (to be able to allocate more kernel memory and also speed up the render slightly), de-comment the second line of each variant.  The x-axis goes from west (negative x) to east (positive x), the y-axis goes from south (negative y) to north (positive y), and the z axis goes from the lowest to the highest elevation (in meters above sea level) of the array.  Since a complete array is padded with zeros, the lowest elevation is almost always sea level (unless you are rendering a quadrangle in Death Valley, California, for instance, with areas below sea level).

(*No box, no axes, no mesh, default lighting.*)

ListPlot3D[elevarray,
(*DisplayFunction->(Display["a landscape",#]&),*)
ViewPoint->{1.50, -1.50, 0.75},(*from southeast*)
PlotLabel->quad,
Boxed->False,
BoxRatios->Automatic,
Axes->False,
Mesh->False,
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},
Background->GrayLevel[0]]

(*This variant adds a box and labels to the display.*)

ListPlot3D[elevarray,
(*DisplayFunction->(Display["a landscape",#]&),*)
ViewPoint->{1, 1, 1.25},(*from northeast*)
PlotLabel->quad,
Boxed->True,
BoxRatios->Automatic,
AxesLabel->{"west to east","south to north","elevations"},
Axes->True,
Mesh->False,
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},
Background->GrayLevel[0]]

(*This variant adds special light sources to the display.*)

ListPlot3D[elevarray,
(*DisplayFunction->(Display["a landscape",#]&),*)
ViewPoint->{0, -1, 1.0},(*from south and above*)
PlotLabel->quad,
Boxed->True,
BoxRatios->Automatic,
AxesLabel->{"west to east","south to north","elevations"},
Axes->True,
Mesh->False,
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},
Background->GrayLevel[0],
LightSources->{{{0,0,10},RGBColor[0.4, 0.4, 0.80]},
{{10,-3,6},RGBColor[1, 0.75, 0.00]}}]

(*This variant shades the landscape with a grayscale,
according to elevation.  If the borders of the elevation array
are first stripped off, to remove those rows and columns
containing zeros, a finer shading results.*)

(*elevarray=Table[Take[elevarray[[i]],{12,360}],{i,12,450}];
q=Dimensions[elevarray];*)

ListPlot3D[elevarray,
(*DisplayFunction->(Display["a landscape",#]&),*)
ViewPoint->{0, -1, 1.0},(*from south and above*)
PlotLabel->quad,
Boxed->True,
BoxRatios->Automatic,
AxesLabel->{"west to east","south to north","elevations"},
Axes->True,
Mesh->False,
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},
Background->GrayLevel[0],
Lighting->False,
ColorFunction->GrayLevel]

(*With appropriately chosen contour intervals, this
function will reconstitute a topographic map from a
DEM.  However, it crashes my computer when array
dimensions exceed about {200,200}.*)

ListContourPlot[elevarray,
PlotLabel->quad,
AspectRatio->Automatic,
Axes->False,
Contours->20,
(*ContourShading->False,*)(*de-comment to show just contour lines*)
(*ColorFunction->Hue,*)
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}}]

(*Shade DEM according to elevation; purely 2D.*)

ListDensityPlot[elevarray,
PlotLabel->quad,
AspectRatio->Automatic,
Axes->False,
Frame->False,
Mesh->False,
(*ColorFunction->Hue,*)
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}}]

Write and Read elevation arrays

With the routines below, one may write the elevation array to a file, and, possibly, reduce the amount of memory needed to render a DEM.

(*Tell Mathematica where to write to, or read from;
in this case, the hard disk named "Grangold"*)

SetDirectory["Grangold:"]

(*Rolf Mertig's method for writing arrays to disk*)
(*Evaluate this cell in order to use the function*)

ArrayWrite[{x_List, y__List, z_List}, fil_String] :=
Block[{}, 
OpenWrite[fil];
WriteString[fil,"{\n"];
PutAppend[x, fil];
Do[WriteString[fil,","];
	PutAppend[{y}[[i]], fil];,
{i,Length[{y}]}];
WriteString[fil,","];
PutAppend[z, fil]; WriteString[fil,"}\n"];
Close[fil]	];

(*Write the array to disk, using the name "terrain"*)

ArrayWrite[elevarray, "terrain"]

(*read a file into the variable elevarray*) 

elevarray= (<< terrain);
q=Dimensions[elevarray]

(*the dimensions in q are essential to the MeshRange spec;
now, render the DEM in ListPlot3D, etc.*)

Read and Render 1-degree DEMs

The U.S.G.S. 1-degree DEMs are available for free over the Internet, but, when uncompressed, consume 9 MB of disk space.  They have the same basic structure as other DEMs, with A-records, B-records, and C-records; but they have 1201 profiles of 1201 elevations, making for over one million elevations altogether.  This places huge demands upon memory (RAM), so the code below allows one to take a subset of the data, by assigning integer variables to the westernmost and easternmost profiles to be sampled, and the southernmost and northernmost elevations to be used from those profiles.

The landscape resolution of 1-degree DEMs is far inferior to that provided by 7.5 minute DEMs, and unfortunately this is especially apparent when only a subset of the data (a small portion of the entire quadrangle) is rendered.  Also complicating matters is that the profiles and elevations within profiles are spaced 3 arc-seconds apart; one must contrive a way to convert from arc-seconds into meters, or the output from ListPlot3D will be distorted (i.e., mountains too low or too high).  Furthermore, on average, and in temperate latitudes, an arc-second of latitude and an arc-second of longitude measure different distances.  This is because except for the equator itself the parallels of latitude are small circles, but meridians of longitude are great circles.  Moreover, an arc-second of latitude at the southern edge of a 1-degree DEM spans a different distance than an arc-second of latitude at the northern edge.

(*Optional*)
Needs["Utilities`MemoryConserve`"](*calls "Share" to reduce memory use*)

(*Tell Mathematica where your 1-degree DEMs are located*)

SetDirectory["Starlight:1degree DEMs"]

dem=OpenRead["chico-e"];

SetStreamPosition[dem,0];
Print["A-Record:"]
tmp=ReadList[dem,Word,3,NullWords->True];
mark=StreamPosition[dem];
SetStreamPosition[dem,0];
Print["Quadrangle Name:  ", quad = FromCharacterCode[ReadList[dem, Byte, mark]]]

SetStreamPosition[dem,145];
Print["DEM level code:  ", Read[ dem,Number] ]

SetStreamPosition[dem,547];
Print["UTM Corners:  ", 
	utm=Partition[Round[ReadList[ dem,Number,8]], 2] ]

SetStreamPosition[dem,739];
Print["Minimum and Maximum elevations of DEM:  ", Round[ReadList[ dem,Number,2] ]]

SetStreamPosition[dem,858];
Print["Number of profiles of DEM:  ", p=Read[ dem,Number] ]
(*1.18 seconds*)

(*Display outline of quadrangle*)

AppendTo[ utm, utm[[1]] ];
Show[Graphics[Line[utm] ],
	PlotLabel->"Quadrangle Outline",
	AspectRatio->Automatic]

(*Gather selected B-records and profiles.
A 1-degree DEM has 1201 profiles of 1201 elevations;
they are numbered from west to east, and south to north.	
Define your desired westernmost, southernmost etc. elevations
by setting the variables west, east, etc. in the line below.*)

west=150;east=300;south=150;north=300;

SetStreamPosition[dem,1024];
	j=1;brec={};elevarray={};
	
	While[j<east,
	
	b=Round[ReadList[dem,Number,9]];
	If[south>1,
	Skip[dem,Number,south];
	e=Round[ReadList[dem,Number,north-south ]];
	Skip[dem,Number,b[[3]]-north];
	];
	
	If[j>=west,
	AppendTo[brec,b];AppendTo[elevarray,e]
	];
	j=b[[2]] ];

Print["Dimensions of the array are",Dimensions[elevarray]];
Close[dem]

(*must transpose to avoid flipped image*)
elevarray=Transpose[elevarray];
q=Dimensions[elevarray];
Print["The dimensions of the final elevation array are ",q]

(*this variant adds a box and labels to the display*)
ListPlot3D[elevs,
ViewPoint->{0, -1, 1.25},(*from south*)
PlotLabel->quad,
Boxed->True,
BoxRatios->Automatic,
AxesLabel->{"west to east","south to north","elevations"},
Axes->True,
Mesh->False,
(*These values in MeshRange an approximation
and ignore the difference in distance between
arc-seconds of latitude and longitude*)
MeshRange->{{0,(70*q[[2]])},{0,(70*q[[1]])}},
Background->GrayLevel[0]]

DEM Sources, on the Internet

California

The following site provides access to 7.5 minute California DEMs, and is best used after 11:00 P.M.  One should know the name of the quadrangle one wishes to download, and then check the "DEM Codes" directory for its code.  For instance, the Dutch Flat quadrangle has code 611.CDO.  It will then be found within the sub-directory which references the letter "D."

ftp://130.166.124.202/pub/ca_dems.1/

Also, the best source for California 7.5 minute DEMS (2300 of them):
http://130.166.124.2/ca_dems.htm

A source for California 7.5 minute DEMs:
http://bard.wr.usgs.gov/htmldir/dems.html

Rocky Mountains

This site holds many 7.5 minute DEMs in Colorado and Wyoming.  They have been compressed.

http://sequoia.cnr.colostate.edu/ftp.html

 Rocky Mountain 7.5 minute DEMs:
http://corrc.cnr.colostate.edu:8000/ftp.html

Appalachia and more

This site conains an extensive archive of 7.5 minute and 1 degree DEMs, especially good for Appalachia.  I found navigation of the database somewhat difficult.

http://eoswww.essc.psu.edu/eosdb.html

Xerox Spectrum

This site forms a repository for various DEMs.

http://jei.umd.edu/jei/spectrum-index.html (links with quad names)
ftp://spectrum.xerox.com/pub/map/dem/ (the Xerox site itself)

Main U.S. Geological Survey mapping page

Follow links to download 1 degree DEMs, or learn more about DEM data structures.

http://www-nmd.usgs.gov/

For information about 1 degree DEMs, see

http://edcwww.cr.usgs.gov/glis/hyper/guide/1_dgr_dem

1-degree DEMs by state:
http://edcwww.cr.usgs.gov/glis/hyper/guide/1_dgr_demfig/states.html

Rendering Landscapes in POV-Ray

POV-Ray is a remarkable free ray-tracing application which represents the combined efforts of many individuals over several years' time.  It may be downloaded from many locations on the Internet; the POV-Ray home page is at http://www.pov.org.  It is available for all common computer platforms, including Unix, Windows, and Macintosh.

Roman Maeder created a Mathematica package called POVray.m, which converts Mathematica Graphics3D and other objects into POV-Ray format.  It is available at http://www.mathconsult.ch/showroom/ray/index.html.

Landscapes in POV-Ray are created by using height fields, and in version 2.0 these in turn were generated from bitmap files in the Targa format.  Version 3.0 of POV-Ray is now available and supports other bitmap file types besides Targa.  I have no experience in using these other types, so I will confine my remarks to using Targa files.

Why render landscapes in POV-Ray when such excellent results can be obtained directly from Mathematica?  Well, in the ray-tracing environment, mountains cast shadows, reflective surfaces reflect objects around them, and a variety of textures can be applied to surfaces.  The rendering model of POV-Ray differs from that of Mathematica in another important respect:  it makes no attempt to create the mesh of triangles which constitute an landscape in its entirety, and thus is never confronted with using RAM or disk space to store these objects.  It creates its height field polygons as it needs them, by testing to see if a ray from the camera might intersect one or more of them.  Thus, given a Targa which contains landscape information (I will get to that subject shortly), POV-Ray can easily handle height fields which contain millions of polygons.  I have rendered landscapes in POV-Ray which contained over 17 million triangles!

The POV-Ray scheme for reading Targa files to create height fields is as follows:  the 24 bits of color in a Targa are divided into 3 bytes of data, for red, green, and blue.  The red and the green bytes, together, can represent 65536 or 2^16 different colors, or, say, the integers between 0 and 65535.  Such a range of integers is much more than sufficient to handle DEM data sets, where elevations are usually given in meters, or even if in feet, we need only 29,028 of these to handle Mount Everest.

Where does Mathematica fit into this, then?  Well, I can speak only to the Macintosh platform, but something analogous is probably workable on other platforms.  Given an elevation array, such as the routines above generate, we may create a ListDensityPlot of the data; an example of this is shown above.  However, this by default uses a grayscale to color the cells of the density plot.  A special ColorFunction is needed.  I appealed for help from the Mathematica users' group, and Jeffrey Adams was kind enough to devise a ColorFunction to do the required transformation.  He also suggested using a Raster rather than ListDensityPlot to display the data.

I have adapted Jeffrey's HeightFieldColorFunction to make it a pure function.  Suppose we have used Mathematica to read one or more DEMs, using the methods in this notebook.  We have obtained, then, a list of lists named "elevarray," of dimensions named "q."  To create a Raster of this data, colored as POV-Ray will expect it, we do the following:

(*Is it possible you need a cup of coffee?  315 seconds for
one 7.5 minute DEM.*)
maxmin={Max[elevarray],Min[elevarray]};
dim = Reverse[q];
low = maxmin[[2]];
high = maxmin[[1]];
Show[ Graphics[
		Raster[elevarray,{{0,0},dim},{low,high}, 
		ColorFunction -> (RGBColor[Quotient[Floor[#*65535], 256]/255.,Mod[Floor[#*65535],256]/255,0]&)]],
AspectRatio->Automatic]

(*Copy scaling information to Clipboard for use in POV-Ray*)
vertratio=maxmin[[1]]/(30*(Min[q]-1)) //N;
horatio=Max[q]/Min[q] //N;
triple={1, vertratio, horatio};
If[Min[q]==q[[2]],triple,triple=Reverse[triple]];
Print["Use this scaling for the POV heightfield:"]
Print["scale <",triple[[1]],", ",triple[[2]],", ",triple[[3]],">"]

Having evaluated this cell to create a Raster, we obtain a strangely speckled image.  Note that the ColorFunction is applied to scaled values, between 0 and 1 inclusive; we wish the red byte to assume any of 256 different brightnesses, so we multiply by 65,535 and then divide by 256, then by 255 to return the value to its proper range; the green byte must take the elevations mod 256 and then divide by 255 to return the value to its proper range.  So we multiply the scaled value by 65535, take this multiple mod 256, then divide by 255.  In this way all the elevations are represented as a sum of the red and green bytes, which is how POV-Ray interprets a Targa.

A supplementary routine calculates the proper scaling to apply to the height field in POV-Ray.

Having created the Raster, which may take fifteen minutes or so for a block of four or six DEMs, the next step is to convert it to a Targa.  Although methods exist to do this from within Mathematica, these have been found to consume too much memory.  My method is this:  I select the Raster and convert it to object PICT format, in the Graph menu.  Object PICT will infallibly assign every cell of the raster to a single colored pixel in a PICT file.  Then, I select the converted Raster and copy it to the Clipboard.  Using the Edit menu, I choose "Convert Clipboard."  I choose the "Save in file" option and save the object PICT to an external file.

The next step is to use Adobe PhotoShop or a shareware program like Graphic Converter to open this object PICT and save it to the Targa format.  When I do this in PhotoShop, I take care to crop the white border around the Raster which was created in Mathematica.  I also take care to save it into the default directory search path of POV-Ray.

To actually render a landscape in POV-Ray, some syntax such as this will provide a starting point:

//This file renders the Dutch Flat area

#include "colors.inc"
#include "textures.inc"
#default { texture { pigment {color rgb <1, 1, 1>} finish {phong 0.001 ambient 0 diffuse 0.8}}}
#max_trace_level 5

camera
{
  location  <0,  1.0, -2> // <X Y Z>
  direction 3*z             // which way are we looking <X Y Z>
  up        y                 // which way is +up <X Y Z>
  right     (4/3)*x             // which way is +right <X Y Z> and aspect ratio
  look_at   <0,  0.0, 0> // point center of view at this point <X Y Z>
}

light_source { <-3000, 500, -300> color  rgb <.8, 0.5, 0> }
light_source { <0, 10000, 0> color rgb <0.4, 0.3, .8> }

object{
         height_field {tga "dfquad" }        
         translate <-.5, 0, -.5>  // center it on the origin
         scale < 1.0, 0.131351, 1.274932 >
}

Some other remarks might be useful:  in the POV-Ray "camera" statement, the "right" vector is given as (4/3)*x.  This defines the aspect ratio of the rendered image, and in this case permits dimensions identical to those of a typical computer display, 640 by 480 pixels, for instance.  Notice also that a default texture is declared near the beginning of the file; in this statement, the landscape (or whatver object enters the scene) is assigned a color of white, with an ambient lighting component of zero, and a diffuse lighting compnent of 0.8.  The defaults in POV-Ray are 0.2 and 0.6, respectively.  The assignment used in this example resulted from an attempt to mimic Mathematica  lighting.

Credits

These routines were made with the help of Gennady Stupakov, a theoretical physicist at the Stanford Linear Accelerator.  Many thanks, Gennady!

Rolf Mertig contributed the valuable method to write elevation arrays to disk.

Dave Wagner supplied the method to write PostScript to disk.

Paul A. Rubin and Mark Evans offered much expert advice.

Jeffrey Adams made the special ColorFunction to create Rasters which are properly colored for use in POV-Ray, as height fields.

Please send comments or suggestions to:
Russell Towle
P.O. Box 141
Dutch Flat, CA 95714
email:  rustybel@foothill.net