In article <3klgmi$noh@newsbf02.news.aol.com> arthurb733@aol.com (ArthurB733) writes:
> I'm new to IDL and to this group. Please forgive me if this subject has
> come up before. I was wondering if anyone knows of existing IDL
> procedures for reading USGS DEMs, or does someone have pointers to offer
> on the topic?
Well, I may have something for you to chew on as well. I wrote these
awhile back to work with a CDROM put out by MicroMap. They have DEM
files on the disk, but I have no idea if it a the "Standard Format".
Anyway here are two programs. CVTDEM gets the data off of the CDROM
and writes it to a binary files (the originals were in ASCII). DEM
then plots the data up.
Look it over, (watch out of the .sig at the end of this message)
pro cvtdem,file
; Open the database
openr,1,'/mnt/USA/DEM/W109W130.DEM'
openw,2,'/work/W109W130.DEM.bin'
; Record structures
in = assoc(1,bytarr(6012))
b = {DEM,lat:0.,lon:0.,z:intarr(120,10)}
out = assoc(2,b)
i = 0
; Loop until EOF
while not eof(1) do begin
; Grab a record
a = in(i)
; Extract the lat and lon fields and convert to floating point
b.lat = float(string(a(0:3)))+float(string(a(4:7)))/60.
b.lon = float(string(a(8:11)))
j = 7
for k=0,9 do begin
for l=0,119 do begin
; Update offset
j = j+5
; Convert altitude field to integer
b.z(119-l,k) = fix(string(a(j:j+4)))
endfor
endfor
out(i) = b
i = i+1
endwhile
close,1
close,2
end
pro dem,print=print,color=color,logscale=logscale
; Process output options
if keyword_set(print) then begin
set_graph,/cps,/land
cos37 = cos(37.5/!radeg)
ydel = (1.-cos37)/2.
pos = [0,ydel,1.,1.-ydel]
image = bytarr(1200,600)
thk = 5
endif else begin
; Make a window to display the image in
set_graph,/x
window,xsiz=1200,ysiz=600,colors=256
pos = [0,0,1.,1.]
thk = 3
endelse
if keyword_set(color) then loadct,color else loadct,0
tvlct,[255,0],[255,0],[255,0],254
; Open the database
openr,1,'/usr0/map-data/usa/0.5-min/W060W086.DEM.bin'
; Record structure
b = {DEM,lat:0.,lon:0.,z:intarr(120,10)}
in = assoc(1,b)
i = 1284
while not eof(1) do begin
; Grab a record
b = in(i)
; Convert lat-lon to x-y
xo = fix(1080-(b.lon-75.)*120) ; one pixel = 0.5 minutes
yo = fix((b.lat-35.)*120+.5) ; one pixel = 0.5 minutes
if((xo gt 1080) or (xo lt 0)) then goto,jump
; Convert altitude to color via a log-scale
if keyword_set(logscale) then begin
sub = byte(alog10(b.z*100.+1.)*35) < 253B
endif else begin
sub = (byte((b.z < 5000)/20)+3B)
endelse
if keyword_set(print) then begin
image(xo,yo) = sub
endif else begin
tv,sub,xo,yo
endelse
JUMP: i = i+1
if(i gt 2200) then goto,next
; Loop until EOF
endwhile
NEXT: close,1
if keyword_set(print) then tv,image,0,pos(1),xsiz=1,ysiz=cos37,/norm
!x.margin = [0,0]
!y.margin = [0,0]
maps,0,0,/cyl,ran=[-85,-75,35,40],/noeras,/hires,grid=1,pos= pos $
,mcolor=254,mthick=thk,gcolor=254,gres=.25,gthick=thk
maps,0,0,/cyl,ran=[-85,-75,35,40],/noeras,/hires,pos=pos $
,mcolor=255,gcolor=255,gres=.25
end
--
Joseph M. Zawodny (KO4LW) NASA Langley Research Center
Internet: j.m.zawodny@larc.nasa.gov MS-475, Hampton VA, 23681-0001
TCP/IP: ko4lw@ko4lw.ampr.org Packet: ko4lw@n4hog.va.usa.na
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
