[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: reading dem




Sylvain, 

Part of the problem is that very few people have any direct experience
with what you are trying.  Some, including me, offered some general
advice, hoping that would help, but reading byte-level DEM's is not
exactly common knowledge.  Also realize that a topic so focussed and
specific as yours is not likely to pique anybody's interst.  In the
future you would do better to describe what you need to do at a higher
level, and help us by providing documentation, like I do now.

The documentation is here:
http://rockyweb.cr.usgs.gov/nmpstds/demstds.html

The DEM's are digital elevation maps, and are stored in 1024 byte
blocks.  The first block contains a type "A" record.  The next blocks
contain a type "B" record which have the actual elevation data.  The
actual data start at offset 144 of the record (where offset 0 is the
beginning); there are 1201x1201 elements stored, 146 in the first
1024-block, and then 170 in the successive 1024-blocks.  Elevation
values are stored as ASCII format=(I6).

The key thing to realize that all the data is in ASCII, separated by
blanks.  Therefore, while we could read each 1024-block in turn, it's
better to just do a formatted READF.  There is no need to work at the
byte level, except at the beginning to get to the right file offset.
Unfortunately IDL can only read 32k elements at a time, so I read a
row at a time as a compromise.

pro readdem250, file, im
  m = 1201L & n = 1201L ;; Really get this from Type A, element 16
  im = lonarr(m, n)     ;; Formally defined as INTEGER*4

  openr, lun, file, /get_lun
  ;; Skip past type A record and 144 bytes of Type B record
  point_lun, lun, 1024L + 144L 

  ;; Read data, one row at a time
  row = im(*,0)
  for i = 0, 1200 do begin
    readf, lun, row, format='(1201(I6))'
    im(*,i) = row
  endfor
  free_lun, lun
end

Of course I haven't tested this, but you can use this as a starting
point.

Craig

P.S.  If you are still interested in the block structure, then you
would read the blocks like this:

bb = bytarr(1024)*nblocks
readu, lun, bb

This avoids some extraneous STRING calls.

-- 
--------------------------------------------------------------------------
Craig B. Markwardt, Ph.D.         EMAIL:    craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
--------------------------------------------------------------------------