[Date Index] [Thread Index] [Author Index]

Re: trying to Import[] USGS-resident DEM data

• To: mathgroup at smc.vnet.net
• Subject: [mg83415] Re: trying to Import[] USGS-resident DEM data
• From: mcmcclur at unca.edu
• Date: Tue, 20 Nov 2007 03:44:44 -0500 (EST)
• References: <fhp2fk$1ov$1@smc.vnet.net>

On Nov 18, 5:00 am, congruentialumina... at yahoo.com wrote:
> I am trying to Import[] and ReliefPlot[] the DEM data resident on the
> USGS site using the following statment:
>
> Import["http://bard.wr.usgs.gov/bard/dem/dems24k/sfbay_10meter/
> gzip_dems/sf_north.dem.gz", "ColorFunction" -> "DarkTerrain",
> "DownsamplingFactor" -> 4]

Most DEMs representing data in the contiguous US decompose
simply into a rectangular array, but the DEM of interest
does not.  I strongly suspect that the basic code for
importation of DEMs is not written quite well enough to
handle this DEM.

DEM is an ascii format, so it's not too bad to write your
own code to read DEMs. The basic file format of a DEM is
explained on Wikipedia and on the USGS site here:
http://edc.usgs.gov/guides/usgs_dem_supplement.html

Based on that, I pieced together the code below.  It's
certainly not perfect; you might notice a bit of strangeness
on the left side of the image, but it should get you close.
Just download and gunzip the file to "sf_north.dem" first.


strm = OpenRead["sf_north.dem"];
contents = Read[strm, String];
Close[strm];

dataString = StringDrop[StringDrop[contents, 1024], -1204];
dataPlus = StringSplit[dataString, Whitespace];
step1 = Split[dataPlus, ! (StringMatchQ[#1, __ ~~ "D" ~~ __]
  || StringMatchQ[#2, __ ~~ "D" ~~ __]) &];
step2 = Map[ToExpression, Last /@ Partition[Rest[step1], 6], {2}];
step3 = MapAt[dd @@ # &,
   MapAt[{dd @@ Drop[#, -4], hh @@ Take[#, -4]} &, step2,
    List /@ Range[Length[step2] - 1]], -1];
step4 = Partition[Flatten[Prepend[step3, hh @@ Take[dataPlus, 4]]],
   2];
maxLength = Max[Length[Last[#]] & /@ step4];
row[{hh[_, _, len_, _], dd[data__]}] := PadRight[{data}, maxLength];
data = Transpose[row /@ step4];

ReliefPlot[data, "ColorFunction" -> "DarkTerrain"]