Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

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

Search the Archive

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"]


  • Prev by Date: Re: RandomChoice does not accept null weights
  • Next by Date: Re: Plot to Plot3D problem
  • Previous by thread: trying to Import[] USGS-resident DEM data
  • Next by thread: Re: trying to Import[] USGS-resident DEM data