Re: reading dem

[Sylvain Carette <sylvainc(at)total.net>]

Aaaaahh... Thank you so much
Whithout knowing it you answered a couple of unnasked questions
(I think I have one for each page of the documentation)

Craig Markwardt wrote:

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
 

Sorry if my requests sound confused; they are... I'm trying not to
trow all the questions floating in my head at the same time but probably
it will get better as time goes. Sorry also if a bit of frustration did
show up through the lines; I was getting anxious to get something done.
Now already feel better...

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.

I did realize that the data was in ascii but my problem was to make
idl realize that...
I tryed readf but could get only the value of the ascii code instead
of the value itself; of course I tryed to use "format=" but couldnt guess
the right syntax - I think I was trying somthing like "format=(6(i))" so...
- again, the documentation could be lot more helpful by providing a couple
of (meaningful) examples.
So I end up thinking it was not the right function especially that
all the dem doc say it is a binary file and I though that formatted input
mean a real ascii file (one you can read). So I can do lot with readf as
long I know how to express the "format" syntax. Now, where I can find a
good table & rule describing formatting (I have absolutely no background
in fortran - completely alien to me. I have the impression I'll need to
refer to that all the time.

 

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.

Sure I will! I knew already that what I was not on the right track.
Now I will be able to "start" the right way so its a big improvement. Just
this small snippet tell me lot. It show me how to do things right. Thanks
also for the gtopo answer; now I have two goods 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
--------------------------------------------------------------------------

Sylvain Carette
VRML designer-composer