Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
arthurb733@aol.com (ArthurB733) wrote: > > ... I was wondering if anyone knows of existing IDL > procedures for reading USGS DEMs, ..... > > Arthur Brandt > art@pima.gov The ENVI layered product, also available from RSI, will read in DEM's and display them. VERY nice package, but real pricey...I took a hack at it once, and found the format a little obscure - the USGS site has the documentation.. rc olsen naval postgraduate school no US gov't endorsement implied.....
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? I wrote the attached short bit of code a while ago to do this. It's not very fancy, but it should do the job for you. Jeff Moersch Astronomy and Space Sciences Cornell University moersch@astrosun.tn.cornell.edu ----cut here---- pro demread, filename, demarr ; Reads in a standard 1201x1201 USGS Digital Elevation Map. hdr='' linehdr='' alt1=intarr(146) alt2=intarr(170) alt3=intarr(170) alt4=intarr(170) alt5=intarr(170) alt6=intarr(170) alt7=intarr(170) alt8=intarr(35) openr, 1, filename readf,1,hdr for i=0,1200 do begin readf,1,linehdr,alt1,format='(a147,146i6)' readf,1,alt2,format='(170i6)' readf,1,alt3,format='(170i6)' readf,1,alt4,format='(170i6)' readf,1,alt5,format='(170i6)' readf,1,alt6,format='(170i6)' readf,1,alt7,format='(170i6)' readf,1,alt8,format='(35i6)' demarr(i,0:145)=alt1 demarr(i,146:315)=alt2 demarr(i,316:485)=alt3 demarr(i,486:655)=alt4 demarr(i,656:825)=alt5 demarr(i,826:995)=alt6 demarr(i,996:1165)=alt7 demarr(i,1166:1200)=alt8 endfor close,1 end ----cut here---- --
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