(******************************************************************* This file was generated automatically by the Mathematica front end. It contains Initialization cells from a Notebook file, which typically will have the same name as this file except ending in ".nb" instead of ".m". This file is intended to be loaded into the Mathematica kernel using the package loading commands Get or Needs. Doing so is equivalent to using the Evaluate Initialization Cells menu command in the front end. DO NOT EDIT THIS FILE. This entire file is regenerated automatically each time the parent Notebook file is saved in the Mathematica front end. Any changes you make to this file will be overwritten. ***********************************************************************) 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];\.1b 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::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." ; \.1b 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->{"-"}]; \.1b 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"]] ]; \.1b 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 ] ]