Creating Topographic Maps from Online Data Using Mathematica

Instructions are provided for downloading topographic data for any region of the contiguous 48 states from the United States Geographical Survey's  World Wide Web site.  Mathematica code to reduce the original 9megabytes of data to a manageable size is included.  The process is illustrated by creating maps of the Lake Superior region.

by Allan A. Struthers

A calculus sequence at Michigan Technological University (MTU) introduces students to the myriad applications of calculus and technology in a Mathematica-based computer laboratory.

The fourth quarter extends the one dimensional ideas covered earlier to three dimensions.  Students normally find different graphical representations (contour, density, and perspective) of functions of two variables confusing.  The United States Geographical Survey (USGS) data available on the World Wide Web let us introduce these representations with topographic maps of our local area.  Favourable student reaction to this lab in Fall 95 indicates the maps helped students understand 3D graphs and visualize  the geometrical meaning of partial derivatives and gradients.

The maps are constructed from "One Degree Digital Elevation Models" (DEM).  Each DEM consists of heights on a regular 1200 by 1200 grid and covers a one degree latitude by one degree longitude spherical rectangle.  At  47 degrees North, the latitude of Michigan Tech,  this gives a height roughly every 100 yards North-South and 71 yards East-West over a 70 by 50 mile rectangle.

The original datasets need to be reduced from 9megabytes to a more manageable size;  ListPlot3D, running overnight on a Sparc 5 workstation with 32 megabytes of Ram, produced a 120 megabyte PostScript image of an entire DEM!  Two functions are defined to conserve memory by reading and simultaneously reducing the size of a DEM:  compress decreases the resolution while maintaining coverage;  crop decreases the coverage while maintaining resolution.

Data Acquisition

The United States Geographical Survey gives links to a wealth of geological data  from their World Wide Web home page at the Uniform Reference Location (URL) 

http://www.usgs.gov

One of the many formats accessible from this page is an interactive, graphical interface for selecting and downloading topographical data in the "One Degree DEM" format.  The procedure for selecting and downloading data is illustrated by describing how to download the DEM which covers the area surrounding Michigan Technological University (MTU).

Open the URL
	
http://sun1.cr.usgs.gov/glis/hyper/guide/1_dgr_demfig/index1m.html
	
for the graphical interface;  this gives a map of the U.S. divided into large regions.  MTU is on the Upper Peninsula of Michigan in the region labelled "LAKE SUPERIOR";  click on this region to zoom in on this portion of the map.  You will see a map of the Upper Peninsula of Michigan divided into sections.  MTU is on the Keweenaw peninsula projecting into Lake Superior in the eastern portion of the section labeled "HANCOCK"; click on this section to bring up a page offering you a choice of the data for the East or West portion of "HANCOCK" in compressed or uncompressed format.

Click on "HANCOCK EAST COMPRESSED" to down load the file. The files are large (9Mb uncompressed, and between  100 Kb and  2Mb compressed) and the server is not particularly fast (averaging around  8 Kb/sec in the middle of the day);  it will probably take around four minutes to download the compressed file which is roughly 400Kb.

Details of downloading and uncompressing will vary with your browser and computer platfrom. For example, if you are using Mosaic the uncompressed file will be saved under the name "hancock-e.gz" and must be uncompressed with the standard gunzip utility.

The data sets are also available on CD-ROM from the US Geological Survey. Information can be obtained from http://www.usgs.gov/pubprod/index.html.

Data Format and Reduction

The format of a "One Degree DEM" is detailed at the URL.
		
	http://sun1.cr.usgs.gov/glis/hyper/guide/1_dgr_dem

DEMs consist of three types of records:  types A, B, and C.  The general characteristics of the DEM are described by the type A record at the start of the file.  In particular the type A record gives the number of type B records included in the DEM.  Each type B record contains equally spaced heights at the same longitude and covers one degree of latitude;  it also contains some header information.  Type C records describe the accuracy of the data.

Standard DEMs consist of 1201 type B records each containing 1201 heights;  DEMs for northern regions of Alaska contain fewer type B records.  With the header a standard format DEM occupies 9Mb of disk space.  The matrix of data for a complete DEM at 1201 by 1201 will strain the memory of most instructional machines and is much too large to use routinely  in Mathematica plotting routines.  The X windows version of Mathematica on a Sparc5 workstation with 32 Mb of RAM can not ListPlot a complete DEM due to memory constraints.  The command line version of Mathematica running overnight on the Sparc station did ListPlot a complete DEM producing an unwieldy 120Mb postscript file.

It is clearly neccessary to reduce the size of the data.  The data is extracted from a standard DEM and reduced to a manageable size by the functions compress and crop defined in the Appendices.

Compressing Data

The function compress defined in the file topoMap.m extracts the data from a standard form DEM and reduces its size by averaging the entries in blocks.  For example, the command

hanE = compress[ "hancock-e", 60 ];

extracts the data from  the DEM "hancock-e" and gives a manageable 60 by 60 dataset called hanE by averaging the heights of the original DEM in 20 by 20 blocks.  On a Unix machine the command looks for the DEM in the Mathematica directory.  On a Macintosh the command looks for the file in the folders in $Path;  you add folders to this list using the command AppendTo[$Path, "folder:subfolder:subsubfolder"].

The messages indicate:  the size of the new data set; the scale of the new dataset; and how to convert coordinates in the new dataset to coordinates for the original DEM.  This operation takes nearly 3 minutes on the Sparc and 8 minutes on a Macintosh Quadra 360.

Having constructed a smaller dataset you can rapidly plot maps as contour plots. A ListDensityPlot is useful as a previewer (it is faster and uses less memory) however the map does not look as clean.  The southern edge of the image is 47 degrees North so the scaling
 
AspectRatio -> 1/Cos[ 47 Pi/180]
		
sets the North-South and East-West scales equal;  the option
		
ColorFunction-> Function[c, Hue[0.54 - c]] 
		
makes Lake Superior blue. The built-in color function Hue gives a red lake.

ListContourPlot[ hanE,
	AspectRatio->1/Cos[47.0 Pi/180],
	PlotRange->All,
	ColorFunction->Function[c, Hue[0.54 - c]] ];

ListPlot3D produces a perspective plot of the same dataset:  The additional vertical scaling in BoxRatios is chosen to produce a reasonable image;  the mesh obscures the plot unless turned off.

ListPlot3D[ hanE,
	BoxRatios->{1, 1/Cos[47.0 Pi/180], 0.15},
	PlotRange->All,
	Mesh->False ];

Cropping Data

Isle Royale National Park occupies a large island in Lake Superior.  The rectangle with corners at  {0, 60} and {18, 52} in the North-West corner of the image contains the central portion of this island.  The function crop defined in topoMap.m  extracts the data from a standard DEM and reduces its size by cropping the coverage region of the data.  For example, the command

central = crop[ "hancock-e", 20 {0, 60}, 20 {18, 52} ];

extracts the data for the rectangle  {0, 1200} and {360, 1040} from the original DEM and names this new dataset central.  The first two messages indicate that central is large and significantly non-square.  The last message indicates how to translate the coordinate system to the original DEM.   At 360 by 160 data points central contains 16 times as many data points as the 60 by 60 dataset hanE.

ListContourPlot[ central,
	PlotRange->All,
	ColorFunction->Function[c, Hue[0.54 - c]] ];

Assembling Mosaics

Isle Royale straddles the Hancock and Thunder Bay DEMs.  The additional  DEMs "hancock-w", "thunder_bay-e" and "thunder_bay-w" were downloaded from the Web site.  All four DEMS were compressed to the same size and plotted.

thuW = compress[ "thunder_bay-w", 40];
thuE = compress[ "thunder_bay-e", 40];
hanW = compress[ "hancock-w", 40];
hanE = compress[ "hancock-e", 40];

thuWpic = ListContourPlot[ thuW,
	AspectRatio->1/Cos[ 47 Pi/180],
	PlotRange->All,
	ColorFunction->Function[c, Hue[0.54 - c]],
	DisplayFunction->Identity,
	FrameTicks->None ];
	
thuEpic = ListContourPlot[ thuE,
	AspectRatio->1/Cos[ 47 Pi/180],
	PlotRange->All,
	ColorFunction->Function[c, Hue[0.54 - c]],
	DisplayFunction->Identity,
	FrameTicks->None ];

hanWpic = ListContourPlot[ hanW,
	AspectRatio->1/Cos[ 47 Pi/180],
	PlotRange->All,
	ColorFunction->Function[c, Hue[0.54 - c]],
	DisplayFunction->Identity,
	FrameTicks->None ];

hanEpic = ListContourPlot[ hanE,
	AspectRatio->1/Cos[ 47 Pi/180],
	PlotRange->All,
	ColorFunction->Function[c, Hue[0.54 - c]],
	DisplayFunction->Identity,
	FrameTicks->None ];

The fragments can be assembled using Show  and GraphicsArray.  The options
	
FrameTicks->None 
GraphicsSpacing->0.0
	
ensure the fragments fit snugly.

Show[ GraphicsArray[ 
	{	{ thuWpic, thuEpic},
		{ hanWpic, hanEpic}	},
		GraphicsSpacing->0.0  ]];

You can assemble sections seamlessly with Show and the MeshRange option.  The data for the central portion of the island is contained in central;  this data set runs from the point {0, 1040} to {360, 1200} in the coordinate system of "hancock-e".  The southern portion of the island occupies the North-East corner of "hancock-w" from  30 {30, 30} to 30 {40, 40}.  The data for the southern part of the island is extracted and plotted.

south = crop["hancock-w", 30 {30, 30}, 30 {40, 40}];

Choosing the shared corner of the 4 DEMs as the origin south covers the rectangle with South-West corner {-300, -300} and North-East corner {0,0} while central covers the rectangle with South-West corner {0, -160} and North-East corner {360, 0}.  The following command pastes these two sections together.  It is necessary to specify the contours because the maximum elevations of the two maps are different.

contourlist = Table[ i, {i, 180, 420, (420 - 180)/8} ];

southpic = ListContourPlot[ south,
	MeshRange-> {{-300, 0}, {-300, 0}},
	PlotRange->All,
	Contours -> contourlist,
	ContourShading->False,
	DisplayFunction->Identity ];
	
centralpic = ListContourPlot[ central,
	MeshRange-> {{0, 360}, {-160, 0} },
	Contours -> contourlist,
	PlotRange->All,
	ContourShading->False,
	DisplayFunction->Identity ];
	
Show[ southpic, centralpic,
	DisplayFunction->$DisplayFunction];

Conclusions

Given sufficient memory it is possible to create high resolution maps of any locale in the contiguous 48 states.  The ability to convert rapidly between contour and perspective plots of familiar terrain will hopefully improve students understanding of graphical representations of functions of two variables.

topoMap.m

Compress

compress::usage = "compress[ DEM, newsize ] reads a standard form 
United States Geological Survey (USGS) One Degree Digital Elevation Model (DEM)
and reduces it from 1201 by 1201 to newsize by newsize.  The image is clipped 
if newsize does not divide 1200.";

compress[ file_String, newsize_:60] := 
Module[ 
{tempstream, oldsize, frac, 
templatebig, temp2, temp3, temp4, tempdata},

tempstream = OpenRead[file];
Skip[ tempstream, Word, 
WordSeparators->{"-"}];

Skip[ tempstream, Word, 39 ];

(* Opens the DEM for reading and skips throught the header. *)

oldsize = Read[ tempstream, Number];
If[ oldsize =!= 1201, 
	Return[Print["Non-standard sized DEM."]]];
If[ newsize > 300, 
	Print[
	"Caution large data set"]
	];
	
(* Checks to see if the DEM is standard sized. *)
(* Quits and prints an appropriate message if it is not. *)
(* Prints a warning message if dataset is large. *)

frac = Floor[ oldsize/newsize ];

Skip[ tempstream, Word, 9];
templatebig = Table[ Number, {frac}, { 1210}];
tempdata = {};

(* Defines the reduction factor frac *)
(* Skips to the start of the first row. *)
(* Sets up templatebig for frac lines of data. *)
(* Initializes tempdata for storage. *)

Do[
temp2 = Read[ tempstream, templatebig];
temp3 = Sum[ Part[temp2, i], {i, frac}];
temp4 = Transpose[Partition[
	Drop[ temp3, - ( 1210 - frac newsize) ], frac]];
AppendTo[tempdata, Sum[ Part[ temp4, i], {i, frac}]],
							{i, newsize}];
							
(* Reads the data in groups of frac lines.*)
(* Writes the sum of frac rows into temp3. *)
(* The nine extra characters and any left overs are discarded *)
(* temp4 is partitioned and summed. *)
(* The average is not computed to   *)
(* stay in integer arithmetic.  *)
(* Transposing etc.takes advantage of storage structure.*)
(* Data is accumulated in tempdata.*)

Close[tempstream];
Print[ 
"Data reduced from 1201 by 1201 to ", 
newsize, " by ", newsize,"."];
Print[newsize, " ticks to a degree."];
Print[ 
"Multiply coordinates by ", frac," to get USGS coords."];



1/frac^2 Transpose[tempdata]

(* Closes the data stream. *)
(* Outputs the normalized *)
(* smaller data set Transposed so that North is *)
(* vertical on the screen. *)

];

Crop

crop::usage = "crop[DEM, {x1, y1}, {x2, y2} ] extracts 
a rectangular section from a standard United States 
Geophysical Survey (USGS) Digital Elevation Model (DEM) 
file in the current directory.  
The size of a standard DEM is 1201 by 1201.  The coordinates 
{x1,y1} and {x2, y2} are the corners of the section in the 
coordinates of the original DEM." ;
 

ClearAll[crop]

crop[ file_String, {x1_, y1_}, {x2_, y2_}] := 

Module[ 
{lowx, lowy, dimx, dimy,
tempstream, templine, oldsize, tempdata, newline},

{lowx, lowy} = {Min[x1, x2], Min[y1, y2]};
{dimx, dimy } = {Abs[x1-x2], Abs[y1- y2]};

(* Finds the SouthWest corner and size of the section. *)

If[ dimx dimy > 201^2, 
Print[ "Caution large dataset."] ];
 
If[ Max[ dimx/dimy, dimy/dimx] > 2.1, 
Print[ "Caution highly non-square dataset."] ];


linetemplate = Table[ Number, {dimy}];

tempdata = {};

(* Opens the data stream and initializes tempdata *)
(* Sets up a template newline for every line *)

tempstream = OpenRead[file];

Skip[ tempstream, Word, 
WordSeparators->{"-"}];

Skip[ tempstream, Word, 39 ];

(* Opens the input stream and skips through the header. *)

oldsize = Read[ tempstream, Number];

If[ oldsize =!= 1201, 
	Return[Print["Non-Standard Size DEM"]] ];
	
If[ Max[x1, x2, y1, y2] > 1200,
	Return[
	  Print[
	   "Requested data outside the coverage of DEM"]]];
	  
(* Reads the data size of the DEM.  Prints a message *)
(* and quits if the DEM is not the standard 1201 by 1201. *)
(* Prints a message and quits if a portion of the requested *)
(* coverage lies outside the covergae of the DEM in the file. *)


Skip[tempstream, Word, 9];
Skip[tempstream, Number, 1210*lowx];
Skip[tempstream, Number, lowy];
(* Skips: the first header; *)
(* the first lowx lines of 1210 numbers; *)
(* and the first lowy entries of the DEM. *)

Do[
	(tempcol = Read[ tempstream, linetemplate];
	AppendTo[ tempdata, tempcol ];
	Skip[ tempstream, Number, 1210 - dimy]),
		{dimx} ];
		
Close[tempstream];
Print[ "Rectangular data set with corners ",
	{x1, y1}, " and ", {x2, y2} ];
Print[ "Add ", {lowx, lowy},
" to translate coordinates to the DEM's."];
Transpose[ tempdata ]
	
]

Acknowledgments

References

About the Author

Allan A. Struthers
Mathematical Sciences
Michigan Technological University
Houghton, MI 49931.

Electronic Subscriptions