(***********************************************************************
Mathematica-Compatible Notebook
This notebook can be used on any computer system with Mathematica 3.0,
MathReader 3.0, or any compatible application. The data for the notebook
starts with the line of stars above.
To get the notebook into a Mathematica-compatible application, do one of
the following:
* Save the data starting with the line of stars above into a file
with a name ending in .nb, then open the file inside the application;
* Copy the data starting with the line of stars above to the
clipboard, then use the Paste menu command inside the application.
Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode. Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).
NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing the
word CacheID, otherwise Mathematica-compatible applications may try to
use invalid cache data.
For more information on notebooks and Mathematica-compatible
applications, contact Wolfram Research:
web: http://www.wolfram.com
email: info@wolfram.com
phone: +1-217-398-0700 (U.S.)
Notebook reader applications are available free of charge from
Wolfram Research.
***********************************************************************)
(*CacheID: 232*)
(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[ 49283, 1217]*)
(*NotebookOutlinePosition[ 50381, 1253]*)
(* CellTagsIndexPosition[ 50337, 1249]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[CellGroupData[{
Cell[TextData[
"Rendering U.S.G.S. DEM Quadrangles"], "Title",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
by
Russell Towle
February, 1997\
\>", "Subsubtitle",
Evaluatable->False,
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData["Discussion"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[{
StyleBox[
"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.\n\
\nDigital topographic maps remedy this deficiency, but the requisite software \
is rare, and good renderings of landscapes even rarer. ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontSlant->"Italic"],
StyleBox[
" includes a variety of plotting functions which can be applied to such \
files, and this notebook contains procedures to read and render \"digital \
landscapes.\"\n\nThe 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 ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontSlant->"Italic"],
StyleBox[
" package might be written, containing a variety of functions and options. \
I welcome suggestions! See the Credits section for how to contact me.\n\n\
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 ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["MathSource",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontSlant->"Italic"],
StyleBox[". This notebook reads and renders, using ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["ListPlot3D",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontWeight->"Bold"],
StyleBox[
" 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.\n\nThe 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.\n\nThe 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:\n\n1. 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.\n\n2. 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. \n\n3. 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.\n\nThe 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 ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["ListPlot3D",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontWeight->"Bold"],
StyleBox[
", the profiles must be padded with zeros above and below, to make the \
array a rectangular list of lists.\n\nOn 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.\n\nThe \
worst problems presented to ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontSlant->"Italic"],
StyleBox[
" 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 ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["ListPlot3D",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontWeight->"Bold"],
StyleBox[
" 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:\n\n1. A subset of the elevation array may \
be taken; the code for this is commented out, in the sections which read \
DEMs, below.\n\n2. 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 ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["ListPlot3D",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontWeight->"Bold"],
StyleBox[
".\n\n3. One may allocate nearly all memory to the kernel, and write the \
output of ",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14],
StyleBox["ListPlot3D",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14,
FontWeight->"Bold"],
StyleBox[
" 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.\n\nI 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.\n\nI welcome any comments, bug reports, or \
suggestions. See the Credits section for information about how to contact \
me.",
Evaluatable->False,
AspectRatioFixed->True,
FontSize->14]
}], "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Read one 7.5' DEM"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData["Example"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData[{
StyleBox["(*Let ",
AspectRatioFixed->True],
StyleBox["Mathematica",
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" know where your DEMs are located;\nin this case, in the folder \"7.5\" \
within the hard disk \"Starlight\"*)\n\nSetDirectory[\"Starlight:7.5\"]",
AspectRatioFixed->True]
}], "Input",
AspectRatioFixed->True],
Cell[OutputFormData["\<\
\"Starlight:7.5\"\
\>",
"\<\
Starlight:7.5\
\>"], "Output",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[TextData[
"(*Construct an array of elevations from the DEM*)\n\nReadDEM[\"Dutch \
Flat\"]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"(*Boxed, with axes, and special light sources*)\n\nListPlot3D[elevarray,\n\
(*DisplayFunction->(Display[\"a landscape\",#]&),*)\nViewPoint->{0, -1, \
1.0},(*from south and above*)\nPlotLabel->quad,\nBoxed->True,\n\
BoxRatios->Automatic,\nAxesLabel->{\"west to east\",\"south to \
north\",\"elevations\"},\nAxes->True,\nMesh->False,\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},\nBackground->GrayLevel[0],\n\
LightSources->{{{0,0,10},RGBColor[0.4, 0.4, 0.80]},\n{{10,-3,6},RGBColor[1, \
0.75, 0.00]}}]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"(*Grayscale image*)\n\nListPlot3D[elevarray,\n\
(*DisplayFunction->(Display[\"a landscape\",#]&),*)\nViewPoint->{0, -1, \
0.5},(*from south and above*)\nPlotLabel->quad,\nBoxed->True,\n\
BoxRatios->Automatic,\nAxesLabel->{\"west to east\",\"south to \
north\",\"elevations\"},\nAxes->True,\nMesh->False,\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},\nColorOutput->GrayLevel,\n\
Background->GrayLevel[1]]"], "Input",
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Initialization"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"(*Cast as a function, with four variables left global:\nelevarray, quad \
(quadrangle name), q (dimensions of array)*)\n\nReadDEM[file_String]:=\n\
Block[{dem,utm,pronum,south,north,capital,\n\
pedestal,pillar,profilelength,elevflag,k},\n(*Create an InputStream object, \
name it \"dem\"*)\ndem=OpenRead[file];\n\n(*Read A-record, retrieve quad \
name, etc.*)\nSetStreamPosition[dem,0];\n\
tmp=ReadList[dem,Word,3,NullWords->True];\nmark=StreamPosition[dem];\n\
SetStreamPosition[dem,0];\nPrint[\"Quadrangle Name: \", \n\t\
quad=FromCharacterCode[ReadList[dem, Byte, mark]] ];\n\t\n\
SetStreamPosition[dem,145];\nPrint[\"DEM Level Code: \", \n\tRead[ \
dem,Number] ];\n\t\nSetStreamPosition[dem,151];\nPrint[\"Elevation Pattern \
Code: \", \n\tRead[ dem,Number] ];\t\t\t\n\nSetStreamPosition[dem,157];\n\
Print[\"Planimetric Ground Reference: \", \n\tRead[ dem,Number] ];\n\n\
SetStreamPosition[dem,529];\nu=Read[ dem,Number];\n\
groundunits=Which[u==0,\"radians\",u==1,\"feet\",\n\t\
u==2,\"meters\",u==4,\"arc-seconds\"];\nPrint[\"Ground Planimetric Units: \
\", groundunits];\n\nSetStreamPosition[dem,535];\nelevflag=Read[ dem,Number];\
\nelevunits=Which[elevflag==1,\"feet\", elevflag==2,\"meters\"];\n\
Print[\"Elevation Units: \", elevunits];\t\n\nSetStreamPosition[dem,547];\n\
Print[\"UTM Corners: \", \n\tutm=Partition[Round[ReadList[ dem,Number,8]], \
2] ];\n\nSetStreamPosition[dem,739];\nPrint[\"Minimum and maximum elevations: \
\",\n\tRound[ReadList[ dem,Number,2] ]];\n\nSetStreamPosition[dem,858];\n\
Print[\"Number of elevation profiles: \", \n\tpronum=Read[ dem,Number] ];\n\n\
(*Display outline of quadrangle*)\nAppendTo[ utm, utm[[1]] ];\n\
Show[Graphics[Line[utm] ],\n\tPlotLabel->quad,\n\tAspectRatio->Automatic];\n \
\nsouth=Min[Map[Take[#, -1]&, utm]];\nnorth=Max[Map[Take[#, -1]&, utm]];\n \n\
(*Read B-records, create elevation array*)\nSetStreamPosition[dem,1024];\n\
k=1;elevarray={};\n\nWhile[k<=pronum,(*while k is less than the number of \
profiles*)\n\nb=Round[ReadList[dem,Number,9]];(*get the header*)\n\n(*b[[3]] \
is number of elevations in profile*)\n(*b[[6]] is the northing of \
southernmost elevation in profile*)\nc=Round[ReadList[dem,Number,b[[3]] ] \
];(*get the elevations in this profile*)\n\
capital=Round[(north-(b[[6]]+(30*(b[[3]]-1))))/30];(*how many zeros above*)\n\
pedestal=Round[((b[[6]])-south)/30];(*how many zeros below*)\n(*join zeros \
below, elevations, zeros above*)\n\
pillar=Flatten[{Table[0,{pedestal}],c,Table[0,{capital}]}];\n\
elevarray={elevarray, pillar};k++];\n\n\nprofilelength=Length[pillar];(*all \
profiles now have equal length*)\nelevarray=Flatten[elevarray];(*discard all \
but one set of braces*)\n\
elevarray=Partition[elevarray,profilelength];(*restore profiles*)\n\n(*Close \
InputStream*)\nClose[dem];\n\n(*If given to ListPlot3D without transposing, \
landscape is mirror image*)\nelevarray=Transpose[elevarray];\n\n\
q=Dimensions[elevarray];(*needed for MeshRange spec*)\nPrint[\"The dimensions \
of the final elevation array are \",q];\n\n(*Finally, if the measure of \
elevations is feet, rather than meters,\nwe must multiply the entire array by \
.3048*)\nIf[elevflag==1,elevarray=elevarray*.3048,elevarray];\n]"], "Input",
InitializationCell->True,
AspectRatioFixed->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Read multiple 7.5' DEMs"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[{
StyleBox[
"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 ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["ListPlot3D",
Evaluatable->False,
AspectRatioFixed->True,
FontWeight->"Bold"],
StyleBox[".",
Evaluatable->False,
AspectRatioFixed->True]
}], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData["Example"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"(*This is one way to conserve memory use by the kernel*)\n\
Needs[\"Utilities`MemoryConserve`\"]"], "Input",
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData[{
StyleBox["(*let ",
AspectRatioFixed->True],
StyleBox["Mathematica",
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" know where your DEMs are located,\nin this case, in the folder \"7.5\" \
within the hard disk \"Starlight\"*)\nSetDirectory[\"Starlight:7.5\"]",
AspectRatioFixed->True]
}], "Input",
AspectRatioFixed->True],
Cell[OutputFormData["\<\
\"Starlight:7.5\"\
\>",
"\<\
Starlight:7.5\
\>"], "Output",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[TextData[
"(*Two DEMs, 174 seconds, Four DEMs, 446 seconds*)\nReadMultipleDEMs[{\"Dutch \
Flat\",\"Westville\"}]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"(*Default lighting, no box, no axes. To write\nPostScript graphics to disk \
rather than displaying,\nde-comment the second line. Other rendering schemes\
\nwill be found in the \"Rendering Variants\" section.*)\n\n\n\
ListPlot3D[elevarray,\n(*DisplayFunction->(Display[\"a landscape\",#]&),*)\n\
ViewPoint->{-0.0, -1.00, 1.0},(*from south and above.*)\nPlotLabel->quad,\n\
Boxed->False,\nBoxRatios->Automatic,\nAxes->False,\nMesh->False,\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},\nBackground->GrayLevel[0]]"],
"Input",
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Initialization"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"ReadMultipleDEMs[dora_List]:=\nBlock[{dem,mark,u,utm,groundunits,elevflag,\n\
m,minmax,p,b,brec,allbrec,pronum,\neasting,longprofile,pronorth,spots,a,\n\
long,south,north,capital,pedestal,pillar,k},\n\nk=1;dem={};\n\
While[k<=Length[dora],\nAppendTo[dem, OpenRead[ dora[[k]] ]];\nk+=1];\n\n\
k=1;utm={};quad={};pronum={};elevunits={};minmax={};\nWhile[k<=Length[dem],\n\
\nSetStreamPosition[dem[[k]],0];\n\
tmp=ReadList[dem[[k]],Word,3,NullWords->True];\n\
mark=StreamPosition[dem[[k]]];\nSetStreamPosition[dem[[k]],0];\n\
qind=FromCharacterCode[ReadList[dem[[k]], Byte, mark]];\nPrint[\"Quadrangle \
Name: \",qind]; \nAppendTo[quad, qind];\n\nSetStreamPosition[dem[[k]],145];\n\
Print[\"DEM Level Code: \", \n\tRead[ dem[[k]],Number] ];\n\t\n\
SetStreamPosition[dem[[k]],151];\nPrint[\"Elevation Pattern Code: \", \n\t\
Read[ dem[[k]],Number] ];\t\t\t\n\nSetStreamPosition[dem[[k]],157];\n\
Print[\"Planimetric Ground Reference: \", \n\tRead[ dem[[k]],Number] ];\t\n\n\
SetStreamPosition[dem[[k]],529];\nu=Read[ dem[[k]],Number];\n\
groundunits=Which[u==0,\"radians\",u==1,\"feet\",\n\t\
u==2,\"meters\",u==4,\"arc-seconds\"];\nPrint[\"Ground Planimetric Units: \
\", groundunits];\n\nSetStreamPosition[dem[[k]],535];\nelevflag=Read[ \
dem[[k]],Number];\neither=Which[elevflag==1,\"feet\", \
elevflag==2,\"meters\"];\nPrint[\"Elevation Units: \", either];\n\
AppendTo[elevunits,elevflag];\n\nSetStreamPosition[dem[[k]],547];\n\
u=Partition[Round[ReadList[dem[[k]],Number,8]],2];\nPrint[\"UTM Corners: \
\",u]; \nAppendTo[utm, u];\n\nSetStreamPosition[dem[[k]],739];\n\
Print[\"Minimum and maximum elevations: \",\n\tm=Round[ReadList[ \
dem[[k]],Number,2] ]];\nAppendTo[minmax, m];\n\n\
SetStreamPosition[dem[[k]],858];\np=Read[ dem[[k]],Number];\nPrint[\"Number \
of elevation profiles: \",p]; \nAppendTo[pronum, p];\n\nPrint[\"\"];\nk+=1];\
\n\n\n(*Display outlines of quadrangles*)\nDo[AppendTo[ utm[[i]],utm[[i,1]] \
],{i, Length[utm]}];\nShow[Graphics[Table[Line[utm[[k]]],{k,Length[utm]}] ],\n\
\tPlotLabel->\"Quadrangle Outlines\",\n\tAspectRatio->Automatic];\n\t\n\t\n\
(*make lists of B-record headers, and elevation profiles*)\n\
allbrec={};profiles={};k=1;\nWhile[k<=Length[dem],\n\
SetStreamPosition[dem[[k]],1024];\n\n\tj=1;brec={};elevs={};\n\t\
While[jTrue,
AspectRatioFixed->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Rendering Variants"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"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)."], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"(*No box, no axes, no mesh, default lighting.*)\n\nListPlot3D[elevarray,\n\
(*DisplayFunction->(Display[\"a landscape\",#]&),*)\nViewPoint->{1.50, -1.50, \
0.75},(*from southeast*)\nPlotLabel->quad,\nBoxed->False,\n\
BoxRatios->Automatic,\nAxes->False,\nMesh->False,\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},\nBackground->GrayLevel[0]]"],
"Input",
AspectRatioFixed->True],
Cell[TextData[
"(*This variant adds a box and labels to the display.*)\n\n\
ListPlot3D[elevarray,\n(*DisplayFunction->(Display[\"a landscape\",#]&),*)\n\
ViewPoint->{1, 1, 1.25},(*from northeast*)\nPlotLabel->quad,\nBoxed->True,\n\
BoxRatios->Automatic,\nAxesLabel->{\"west to east\",\"south to \
north\",\"elevations\"},\nAxes->True,\nMesh->False,\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},\nBackground->GrayLevel[0]]"],
"Input",
AspectRatioFixed->True],
Cell[TextData[
"(*This variant adds special light sources to the display.*)\n\n\
ListPlot3D[elevarray,\n(*DisplayFunction->(Display[\"a landscape\",#]&),*)\n\
ViewPoint->{0, -1, 1.0},(*from south and above*)\nPlotLabel->quad,\n\
Boxed->True,\nBoxRatios->Automatic,\nAxesLabel->{\"west to east\",\"south to \
north\",\"elevations\"},\nAxes->True,\nMesh->False,\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},\nBackground->GrayLevel[0],\n\
LightSources->{{{0,0,10},RGBColor[0.4, 0.4, 0.80]},\n{{10,-3,6},RGBColor[1, \
0.75, 0.00]}}]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"(*This variant shades the landscape with a grayscale,\naccording to \
elevation. If the borders of the elevation array\nare first stripped off, to \
remove those rows and columns\ncontaining zeros, a finer shading results.*)\n\
\n(*elevarray=Table[Take[elevarray[[i]],{12,360}],{i,12,450}];\n\
q=Dimensions[elevarray];*)\n\nListPlot3D[elevarray,\n\
(*DisplayFunction->(Display[\"a landscape\",#]&),*)\nViewPoint->{0, -1, \
1.0},(*from south and above*)\nPlotLabel->quad,\nBoxed->True,\n\
BoxRatios->Automatic,\nAxesLabel->{\"west to east\",\"south to \
north\",\"elevations\"},\nAxes->True,\nMesh->False,\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}},\nBackground->GrayLevel[0],\n\
Lighting->False,\nColorFunction->GrayLevel]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"(*With appropriately chosen contour intervals, this\nfunction will \
reconstitute a topographic map from a\nDEM. However, it crashes my computer \
when array\ndimensions exceed about {200,200}.*)\n\n\
ListContourPlot[elevarray,\nPlotLabel->quad,\nAspectRatio->Automatic,\n\
Axes->False,\nContours->20,\n(*ContourShading->False,*)(*de-comment to show \
just contour lines*)\n(*ColorFunction->Hue,*)\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}}]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"(*Shade DEM according to elevation; purely 2D.*)\n\n\
ListDensityPlot[elevarray,\nPlotLabel->quad,\nAspectRatio->Automatic,\n\
Axes->False,\nFrame->False,\nMesh->False,\n(*ColorFunction->Hue,*)\n\
MeshRange->{{0,(30*q[[2]])},{0,(30*q[[1]])}}]"], "Input",
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Write and Read elevation arrays"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"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."], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData[{
StyleBox["(*Tell ",
AspectRatioFixed->True],
StyleBox["Mathematica",
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" where to write to, or read from;\nin this case, the hard disk named \
\"Grangold\"*)\n\nSetDirectory[\"Grangold:\"]",
AspectRatioFixed->True]
}], "Input",
AspectRatioFixed->True],
Cell[OutputFormData["\<\
\"Grangold:\"\
\>", "\<\
Grangold:\
\>"], "Output",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[TextData[
"(*Rolf Mertig's method for writing arrays to disk*)\n(*Evaluate this cell in \
order to use the function*)\n\nArrayWrite[{x_List, y__List, z_List}, \
fil_String] :=\nBlock[{}, \nOpenWrite[fil];\nWriteString[fil,\"{\\n\"];\n\
PutAppend[x, fil];\nDo[WriteString[fil,\",\"];\n\tPutAppend[{y}[[i]], fil];,\n\
{i,Length[{y}]}];\nWriteString[fil,\",\"];\nPutAppend[z, fil]; \
WriteString[fil,\"}\\n\"];\nClose[fil]\t];"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"(*Write the array to disk, using the name \"terrain\"*)\n\n\
ArrayWrite[elevarray, \"terrain\"] "], "Input",
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData[
"(*read a file into the variable elevarray*) \n\nelevarray= (<< terrain);\n\
q=Dimensions[elevarray]\n\n(*the dimensions in q are essential to the \
MeshRange spec;\nnow, render the DEM in ListPlot3D, etc.*)"], "Input",
AspectRatioFixed->True],
Cell[OutputFormData["\<\
{472, 371}\
\>", "\<\
{472, 371}\
\>"], "Output",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Read and Render 1-degree DEMs"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"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.\n\nThe 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."], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"(*Optional*)\nNeeds[\"Utilities`MemoryConserve`\"](*calls \"Share\" to \
reduce memory use*)"], "Input",
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData[{
StyleBox["(*Tell ",
AspectRatioFixed->True],
StyleBox["Mathematica",
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" where your 1-degree DEMs are located*)\n\n\
SetDirectory[\"Starlight:1degree DEMs\"]",
AspectRatioFixed->True]
}], "Input",
AspectRatioFixed->True],
Cell[OutputFormData["\<\
\"Starlight:1degree DEMs\"\
\>",
"\<\
Starlight:1degree DEMs\
\>"], "Output",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData[
"dem=OpenRead[\"chico-e\"];\n\nSetStreamPosition[dem,0];\n\
Print[\"A-Record:\"]\ntmp=ReadList[dem,Word,3,NullWords->True];\n\
mark=StreamPosition[dem];\nSetStreamPosition[dem,0];\nPrint[\"Quadrangle \
Name: \", quad = FromCharacterCode[ReadList[dem, Byte, mark]]]\n\n\
SetStreamPosition[dem,145];\nPrint[\"DEM level code: \", Read[ dem,Number] ]\
\n\nSetStreamPosition[dem,547];\nPrint[\"UTM Corners: \", \n\t\
utm=Partition[Round[ReadList[ dem,Number,8]], 2] ]\n\n\
SetStreamPosition[dem,739];\nPrint[\"Minimum and Maximum elevations of DEM: \
\", Round[ReadList[ dem,Number,2] ]]\n\nSetStreamPosition[dem,858];\n\
Print[\"Number of profiles of DEM: \", p=Read[ dem,Number] ]\n(*1.18 \
seconds*) "], "Input",
AspectRatioFixed->True],
Cell[OutputFormData["\<\
InputStream[\"chico-e\", 4]\
\>",
"\<\
InputStream[chico-e, 4]\
\>"], "Output",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"A-Record\nQuadrangle Name: CHICO - E\nDEM level code: 3\nUTM Corners: \
{{-435600, 140400}, {-435600, 144000}, {-432000, 144000}, {-432000, 140400}}\n\
Minimum and Maximum elevations of DEM: {251, 2809}\nNumber of profiles of \
DEM: 1201"], "Print",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[TextData[
"(*Display outline of quadrangle*)\n\nAppendTo[ utm, utm[[1]] ];\n\
Show[Graphics[Line[utm] ],\n\tPlotLabel->\"Quadrangle Outline\",\n\t\
AspectRatio->Automatic]\n\t"], "Input",
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData[
"(*Gather selected B-records and profiles.\nA 1-degree DEM has 1201 profiles \
of 1201 elevations;\nthey are numbered from west to east, and south to north.\
\t\nDefine your desired westernmost, southernmost etc. elevations\nby setting \
the variables west, east, etc. in the line below.*)\n\n\
west=150;east=300;south=150;north=300;\n\nSetStreamPosition[dem,1024];\n\t\
j=1;brec={};elevarray={};\n\t\n\tWhile[j1,\n\tSkip[dem,Number,south];\n\t\
e=Round[ReadList[dem,Number,north-south ]];\n\tSkip[dem,Number,b[[3]]-north];\
\n\t];\n\t\n\tIf[j>=west,\n\tAppendTo[brec,b];AppendTo[elevarray,e]\n\t];\n\t\
j=b[[2]] ];\n\nPrint[\"Dimensions of the array are\",Dimensions[elevarray]];\n\
Close[dem]\n"], "Input",
AspectRatioFixed->True],
Cell[TextData["Dimensions of the array are{150, 150}"], "Print",
Evaluatable->False,
AspectRatioFixed->True],
Cell[OutputFormData["\<\
\"chico-e\"\
\>", "\<\
chico-e\
\>"], "Output",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData[
"(*must transpose to avoid flipped image*)\nelevarray=Transpose[elevarray];\n\
q=Dimensions[elevarray];\nPrint[\"The dimensions of the final elevation array \
are \",q]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"The dimensions of the final elevation array are {150, 150}"], "Print",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[TextData[
"(*this variant adds a box and labels to the display*)\nListPlot3D[elevs,\n\
ViewPoint->{0, -1, 1.25},(*from south*)\nPlotLabel->quad,\nBoxed->True,\n\
BoxRatios->Automatic,\nAxesLabel->{\"west to east\",\"south to \
north\",\"elevations\"},\nAxes->True,\nMesh->False,\n(*These values in \
MeshRange an approximation\nand ignore the difference in distance between\n\
arc-seconds of latitude and longitude*)\n\
MeshRange->{{0,(70*q[[2]])},{0,(70*q[[1]])}},\nBackground->GrayLevel[0]]"],
"Input",
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["DEM Sources, on the Internet"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData["California"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
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\
\>", "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Rocky Mountains"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
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\
\>", "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Appalachia and more"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"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.\n\nhttp://eoswww.essc.psu.edu/eosdb.html"], "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Xerox Spectrum"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
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)
\
\>", "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Main U.S. Geological Survey mapping page"], "Subsection",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
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\
\>",
"Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[TextData["Rendering Landscapes in POV-Ray"], "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[{
StyleBox[
"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.\n\nRoman Maeder created a ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[" package called POVray.m, which converts ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" Graphics3D and other objects into POV-Ray format. It is available at ",
Evaluatable->False,
AspectRatioFixed->True],
"http://www.mathconsult.ch/showroom/ray/index.html",
StyleBox[
".\n\nLandscapes 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.\n\nWhy render landscapes \
in POV-Ray when such excellent results can be obtained directly from ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
"? 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 ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" 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!\n\n\
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.\n\nWhere does ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" 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 ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" 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.\n\nI have adapted Jeffrey's \
HeightFieldColorFunction to make it a pure function. Suppose we have used ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" 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:",
Evaluatable->False,
AspectRatioFixed->True]
}], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[
"(*Is it possible you need a cup of coffee? 315 seconds for\none 7.5 minute \
DEM.*)\nmaxmin={Max[elevarray],Min[elevarray]};\ndim = Reverse[q];\nlow = \
maxmin[[2]];\nhigh = maxmin[[1]];\nShow[ Graphics[\n\t\t\
Raster[elevarray,{{0,0},dim},{low,high}, \n\t\tColorFunction -> \
(RGBColor[Quotient[Floor[#*65535], \
256]/255.,Mod[Floor[#*65535],256]/255,0]&)]],\nAspectRatio->Automatic]\n"],
"Input",
AspectRatioFixed->True],
Cell[CellGroupData[{
Cell[TextData[
"(*Copy scaling information to Clipboard for use in POV-Ray*)\n\
vertratio=maxmin[[1]]/(30*(Min[q]-1)) //N;\nhoratio=Max[q]/Min[q] //N;\n\
triple={1, vertratio, horatio};\n\
If[Min[q]==q[[2]],triple,triple=Reverse[triple]];\nPrint[\"Use this scaling \
for the POV heightfield:\"]\nPrint[\"scale <\",triple[[1]],\", \
\",triple[[2]],\", \",triple[[3]],\">\"]"], "Input",
AspectRatioFixed->True],
Cell[TextData[
"Use this scaling for the POV heightfield:\nscale <1.24538, 0.353263, 1>"],
"Print",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[TextData[{
StyleBox[
"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.\n\nA \
supplementary routine calculates the proper scaling to apply to the height \
field in POV-Ray.\n\nHaving 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 ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
", these have been found to consume too much memory. My method is this: I \
select the Raster and convert it to ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["object",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
" 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.\n\nThe 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 ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[
". I also take care to save it into the default directory search path of \
POV-Ray.\n\nTo actually render a landscape in POV-Ray, some syntax such as \
this will provide a starting point:",
Evaluatable->False,
AspectRatioFixed->True]
}], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[StyleBox[
"//This file renders the Dutch Flat area\n\n#include \"colors.inc\"\n#include \
\"textures.inc\"\n#default { texture { pigment {color rgb <1, 1, 1>} finish \
{phong 0.001 ambient 0 diffuse 0.8}}}\n#max_trace_level 5\n\ncamera\n{\n \
location <0, 1.0, -2> // \n direction 3*z // which way \
are we looking \n up y // which way is +up \n right (4/3)*x // which way is +right and \
aspect ratio\n look_at <0, 0.0, 0> // point center of view at this point \
\n}\n\nlight_source { <-3000, 500, -300> color rgb <.8, 0.5, 0> }\n\
light_source { <0, 10000, 0> color rgb <0.4, 0.3, .8> }\n\nobject{\n \
height_field {tga \"dfquad\" } \n translate <-.5, 0, -.5> // \
center it on the origin\n scale < 1.0, 0.131351, 1.274932 >\n}",
Evaluatable->False,
AspectRatioFixed->True,
FontWeight->"Bold"]], "Text",
Evaluatable->False,
AspectRatioFixed->True],
Cell[TextData[{
StyleBox[
"\nSome 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 ",
Evaluatable->False,
AspectRatioFixed->True],
StyleBox["Mathematica",
Evaluatable->False,
AspectRatioFixed->True,
FontSlant->"Italic"],
StyleBox[" lighting.",
Evaluatable->False,
AspectRatioFixed->True]
}], "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]],
Cell[CellGroupData[{
Cell["Credits", "Section",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
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\
\>", "Text",
Evaluatable->False,
AspectRatioFixed->True]
}, Closed]]
}, Open ]]
},
FrontEndVersion->"Macintosh 3.0",
ScreenRectangle->{{0, 800}, {0, 580}},
AutoGeneratedPackage->None,
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{640, 563},
WindowMargins->{{44, Automatic}, {Automatic, 1}},
PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}},
ShowCellLabel->True,
ShowCellTags->False,
RenderingOptions->{"ObjectDithering"->True,
"RasterDithering"->False},
Magnification->1,
MacintoshSystemPageSetup->"\<\
00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001
0000I00000400`<300000BL?00400@00000000000000060001T1T00000000000
00000000000000000000000000000000\>"
]
(***********************************************************************
Cached data follows. If you edit this Notebook file directly, not using
Mathematica, you must remove the line containing CacheID at the top of
the file. The cache data will then be recreated when you save this file
from within Mathematica.
***********************************************************************)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[1731, 51, 110, 3, 136, "Title",
Evaluatable->False],
Cell[1844, 56, 110, 6, 89, "Subsubtitle",
Evaluatable->False],
Cell[CellGroupData[{
Cell[1979, 66, 87, 2, 50, "Section",
Evaluatable->False],
Cell[2069, 70, 8692, 180, 70, "Text",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[10798, 255, 94, 2, 30, "Section",
Evaluatable->False],
Cell[CellGroupData[{
Cell[10917, 261, 87, 2, 70, "Subsection",
Evaluatable->False],
Cell[CellGroupData[{
Cell[11029, 267, 375, 11, 70, "Input"],
Cell[11407, 280, 133, 7, 70, "Output",
Evaluatable->False]
}, Closed]],
Cell[11555, 290, 132, 3, 70, "Input"],
Cell[11690, 295, 557, 9, 70, "Input"],
Cell[12250, 306, 460, 8, 70, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[12747, 319, 94, 2, 70, "Subsection",
Evaluatable->False],
Cell[12844, 323, 3289, 48, 70, "Input",
InitializationCell->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[16182, 377, 100, 2, 30, "Section",
Evaluatable->False],
Cell[16285, 381, 706, 19, 70, "Text",
Evaluatable->False],
Cell[CellGroupData[{
Cell[17016, 404, 87, 2, 70, "Subsection",
Evaluatable->False],
Cell[17106, 408, 150, 3, 70, "Input"],
Cell[CellGroupData[{
Cell[17281, 415, 373, 11, 70, "Input"],
Cell[17657, 428, 133, 7, 70, "Output",
Evaluatable->False]
}, Closed]],
Cell[17805, 438, 155, 3, 70, "Input"],
Cell[17963, 443, 573, 9, 70, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[18573, 457, 94, 2, 70, "Subsection",
Evaluatable->False],
Cell[18670, 461, 4436, 66, 70, "Input",
InitializationCell->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[23155, 533, 95, 2, 30, "Section",
Evaluatable->False],
Cell[23253, 537, 1050, 15, 70, "Text",
Evaluatable->False],
Cell[24306, 554, 399, 7, 70, "Input"],
Cell[24708, 563, 465, 8, 70, "Input"],
Cell[25176, 573, 569, 9, 70, "Input"],
Cell[25748, 584, 782, 12, 70, "Input"],
Cell[26533, 598, 501, 8, 70, "Input"],
Cell[27037, 608, 296, 5, 70, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[27370, 618, 108, 2, 30, "Section",
Evaluatable->False],
Cell[27481, 622, 213, 4, 70, "Text",
Evaluatable->False],
Cell[CellGroupData[{
Cell[27719, 630, 347, 11, 70, "Input"],
Cell[28069, 643, 124, 6, 70, "Output",
Evaluatable->False]
}, Closed]],
Cell[28208, 652, 466, 7, 70, "Input"],
Cell[28677, 661, 150, 3, 70, "Input"],
Cell[CellGroupData[{
Cell[28852, 668, 260, 4, 70, "Input"],
Cell[29115, 674, 122, 6, 70, "Output",
Evaluatable->False]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[29286, 686, 106, 2, 30, "Section",
Evaluatable->False],
Cell[29395, 690, 1630, 23, 70, "Text",
Evaluatable->False],
Cell[31028, 715, 146, 3, 70, "Input"],
Cell[CellGroupData[{
Cell[31199, 722, 316, 11, 70, "Input"],
Cell[31518, 735, 151, 7, 70, "Output",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[31706, 747, 757, 12, 70, "Input"],
Cell[32466, 761, 153, 7, 70, "Output",
Evaluatable->False],
Cell[32622, 770, 318, 6, 70, "Print",
Evaluatable->False]
}, Closed]],
Cell[32955, 779, 218, 4, 70, "Input"],
Cell[CellGroupData[{
Cell[33198, 787, 810, 12, 70, "Input"],
Cell[34011, 801, 112, 2, 70, "Print",
Evaluatable->False],
Cell[34126, 805, 120, 6, 70, "Output",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[34283, 816, 221, 4, 70, "Input"],
Cell[34507, 822, 134, 3, 70, "Print",
Evaluatable->False]
}, Closed]],
Cell[34656, 828, 538, 9, 70, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[35231, 842, 105, 2, 30, "Section",
Evaluatable->False],
Cell[CellGroupData[{
Cell[35361, 848, 90, 2, 46, "Subsection",
Evaluatable->False],
Cell[35454, 852, 671, 16, 206, "Text",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[36162, 873, 95, 2, 46, "Subsection",
Evaluatable->False],
Cell[36260, 877, 285, 10, 110, "Text",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[36582, 892, 99, 2, 46, "Subsection",
Evaluatable->False],
Cell[36684, 896, 279, 5, 78, "Text",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[37000, 906, 94, 2, 46, "Subsection",
Evaluatable->False],
Cell[37097, 910, 249, 8, 94, "Text",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[37383, 923, 120, 2, 46, "Subsection",
Evaluatable->False],
Cell[37506, 927, 372, 15, 174, "Text",
Evaluatable->False]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[37927, 948, 108, 2, 30, "Section",
Evaluatable->False],
Cell[38038, 952, 4792, 104, 602, "Text",
Evaluatable->False],
Cell[42833, 1058, 444, 8, 192, "Input"],
Cell[CellGroupData[{
Cell[43302, 1070, 410, 7, 117, "Input"],
Cell[43715, 1079, 150, 4, 36, "Print",
Evaluatable->False]
}, Closed]],
Cell[43880, 1086, 2552, 53, 356, "Text",
Evaluatable->False],
Cell[46435, 1141, 1005, 16, 446, "Text",
Evaluatable->False],
Cell[47443, 1159, 1028, 23, 128, "Text",
Evaluatable->False]
}, Closed]],
Cell[CellGroupData[{
Cell[48508, 1187, 74, 2, 30, "Section",
Evaluatable->False],
Cell[48585, 1191, 670, 22, 286, "Text",
Evaluatable->False]
}, Closed]]
}, Open ]]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)