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