imagetoGSLIBpts.pro

imagetoGSLIBpts.pro

posted 01-18-2002 01:31 PM
Hi all. Here’s the script we’ll be discussing. It takes an ENVI image and outputs it in a specific ASCII format (one pixel per line, with coordinates, and with the first pixel at the *lower* left-hand corner). This format is necessary for the geostatistical program GSLIB.
Note that I’m posting my script below inside “code” tags. Like this:

<code>
text, blah blah blah
</code>

…except that you replace the < > with [ ].

Nick

PS: Also click “disable smilies” when you post, or your code might have little happy faces in it.

code:

; PRO IMAGETOGSLIBPTS
;
; Purpose: To convert an ENVI image into a GSLIB points file.
;
; Started Oct. 15, 2001.
; written by Nick Matzke, Oct 15, 2001.
pro imagetoGSLIBpts

;SETTINGS SETTINGS SETTINGS SETTINGS

;get time that classification begins
starttime = systime(1)

;FOR FULL DMSP SCENE

; use the envi_select wwidget to select the image that contains the data,
; fid will be file id of the open file
; pos will be the number of bands
envi_select, title=’Input Filename of Image’, fid=fid, pos=pos, dims=dims ;, /NO_SPEC, /NO_DIMS

pos1=pos+1 ; add 1 to each position

twa=where(pos1,count1) ; count bands to be used
print,’nb count = ‘,count1

;use envi_file_query to find out number of rows, samples, bands,
;byte order, and note interleave: 0= bsq, 1=bil, 2=bip
;note data type: 1=byt, 2=integer, 4=floating point
; (see progguid.pdf page 330 for more)
envi_file_query, fid, data_type=data_type, nb=nb, interleave=interleave, sname=imagename, fname=fullpath;, h_map=h_map
;handle_value, h_map, map_info
map_info = envi_get_map_info(fid=fid)
proj = envi_get_projection(fid=fid)

help,map_info
print,map_info
help,proj
print,proj
mc=map_info.mc ; This gets the upper left corner info (reference corner, actually)
ps=map_info.ps

help,ps
print,’ps=’
print,ps(0) ;This is the size of the pixel in x
print,ps(1) ;This is the size of the pixel in y
help,mc
print,’mc =’
print,mc(0) ; UL corner pixel #, x (usually 0)
print,mc(1) ; UL corner pixel #, y (usually 0)
print,mc(2) ; UL map coordinate, x
print,mc(3) ; UL map coordinate, y

if ( (mc(0) ne 0) OR (mc(1) ne 0) ) then begin
print, ‘Upper left pixel with coordinate is not pixel 0,0’
print, ‘May God have mercy on your soul.’
stop
endif else begin
twa=2
endelse

;calculate the number of lines and samples of the input image (this deals with subsetting)
ns = dims(2)-dims(1)+1
nl = dims(4)-dims(3)+1

;use pickfile to select the output image (starting path will be same as input ascii rulesfile)
outimage = dialog_pickfile(title=’Enter Filename of Output File (NO *.dat!!)’, get_path=gp, path=filename)

nb=nb
ncb=count1 ; ncb=number of counted bands

;open a metadata file that has the same name as the output image
;with a .meta at the end and start filli ng it
openw,u1,outimage+’.dat’, /GET_LUN ; output ASCII data file
openw,u2,outimage+’-meta.txt’, /GET_LUN ; metadata file
openw,u3,outimage+’.img’, /GET_LUN ; metadata file
;BE CAREFUL WITH openw NOT TO BLAST A GOOD FILE

date_time = systime()
;print to file (= printf)
printf,u2,’Metadata file for IMAGE TO GSLIB (ASCII) run, ‘,date_time
printf,u2, ‘ ‘
printf,u2, ‘Input image : ‘,imagename
printf,u2, ‘Input image pathname: ‘,fullpath
printf,u2, ‘Number of lines : ‘, nl
printf,u2, ‘Number of samples : ‘, ns
printf,u2, ‘Number of bands : ‘, nb
printf,u2, ‘Data type : ‘, data_type
printf,u2, ‘ ‘
printf,u2, ‘Output classified image : ‘,outimage
printf,u2, ‘Number of lines : ‘, nl
printf,u2, ‘Number of samples : ‘, ns
printf,u2, ‘Number of bands : ‘, nb
printf,u2, ‘Number of counted bands : ‘, ncb
printf,u2, ‘Data type : byte’

;print to screen
print,’Metadata file for IMAGE TO GSLIB (ASCII) run, ‘,date_time
print, ‘ ‘
print, ‘Input image : ‘,imagename
print, ‘Input image pathname: ‘,fullpath
print, ‘Number of lines : ‘, nl
print, ‘Number of samples : ‘, ns
print, ‘Number of bands : ‘, nb
print, ‘Number of counted bands : ‘, ncb
print, ‘Data type : ‘, data_type

print, ‘ ‘
print, ‘Output classified image : ‘,outimage
print, ‘Number of lines : ‘, nl
print, ‘Number of samples : ‘, ns
print, ‘Number of bands : ‘, 1
print, ‘Data type : byte’

printf,u2, ‘ ‘
printf,u2, ‘ENVI Image bands used:
;printf,u2, FORMAT='(I2)’, envivars
printf,u2, ‘ ‘

print, ‘ ‘
print, ‘ENVI Image bands used:
;print, FORMAT='(I2)’, envivars
print, ‘ ‘

;WON’T NEED THIS
;set up and write the ASCII header for output image
bandname=strarr(1)
bandname[0]=’binary_class (‘+imagename+’)’
envi_setup_head, data_type=4, $
descrip=’Output image from imagetoGSLIB.PRO using bands of image’, $
file_type=0, fname=outimage+’.hdr’, interleave=0, nl=nl, ns=ns, nb=ncb, map_info=map_info, $
xstart=dims(1), ystart=dims(3), bnames=bandname, /WRITE
;nb = nb. Original nb=1

;make an array with the same dimensions as original image
inputimage = findgen(ns,nl,ncb) ;floating point array
S = SIZE(arr)
filename=strarr(1)

;fs = ‘printf,3,FORMAT=”(1A,’ + string(count1) + ‘F15.6)”,roi_names(i),final_array(0:count1-1,ystart+j)’
;fs = ‘
;printf,3,FORMAT='(1A, F15.6)’,roi_names(i),final_array(0:count1-1,ystart+j)

z=’Input image, ‘
a=’ ns=’
aa=string(ns, FORMAT='(I5)’)
b=’ nl=’
bb=string(nl, FORMAT='(I5)’)
c=’ nb=’
cc=string(ncb, FORMAT='(I5)’)
y=imagename

;printf,u1, ‘Input image : ‘,imagename, ‘ ns=’,ns,’ nl=’,nl,’ nb=’,ncb ;1st line GSLIB
fs = ‘printf,u1,FORMAT=”(8A)”,z+y+a+aa+b+bb+c+cc’
r=execute(fs)

printf,u1, FORMAT='(1I2.0)’,ncb+3 ;2nd line GSLIB (+2 cols for coords, +1 for index ID)
printf,u1, ‘index ID’ ;3nd line GSLIB
printf,u1, ‘Easting (m)’ ;4nd line GSLIB
printf,u1, ‘Northing (m)’ ;5th line GSLIB
for band=0,ncb-1 do begin
;fs = ‘printf,u1,FORMAT=”(8A)”,’
;r=execute(fs)
dd = string(pos1[band], FORMAT='(I3.3)’) ;Band headers
printf,u1,’b’,dd ;band # used
endfor

help,inputimage
print,’Putting data in output matrix’
;retrieve image data for only those bands used in the tree
;this reduces the data volume read into memory to the bands that will actually be used
for band=0,ncb-1 do begin
;****double check that this is getting the right bands****
inputimage[*,*,band] = ENVI_GET_DATA(dims=dims,fid=fid,pos=band) ;idlvars[i])
print,’band output = ‘,band
endfor
help,inputimage
;stop

;write the array to the output image
print,’writing output image file’

;Create array with UTM coords
;map corner and pixelsize
;
;oa is output array
;oax is UTM coord x in middle of each pixel
;oay is UTM coord y in middle of each pixel

;oa=inputimage

oa=REVERSE(inputimage,2, /OVERWRITE) ;OVERWRITES old array

oax=findgen(ns) ; x coordinates — 1 row
oay=findgen(nl) ; y coordinates — 1 column

;findgen = floating array with index values
;findarr = floating array with zeros or random values

;Fill in X coords
for samp=0,ns-1 do begin
oax[samp]=(samp*ps(0))
endfor
oax=oax+mc(2)+0.5*ps(0) ; Add UL coordinate plus 1/2 of 1 pixel (to get middle of pixel)

;Fill in Y coords
for line=0,nl-1 do begin
oay[line]=(line*ps(1))
endfor
oay=oay+mc(3)+0.5*ps(1) ; Add UL coordinate plus 1/2 of 1 pixel (to get middle of pixel)

i=float(1)
;for samp, line, allbands
for line=0,nl-1 do begin
for samp=0,ns-1 do begin

;fs = format statments
;fs = ‘printf,u1,FORMAT=”(1I4.0,” “,1I4.0,” “,’ + string(ncb) + ‘(D10.4, :, ” “))”,su,lu,oa[samp,line,*]’
fs = ‘printf,u1,FORMAT=”(1I7.0,” “,D10.2,” “,D10.2,” “,’ + string(ncb) + ‘(D13.4, :, ” “))”,i,oax[samp],oay[line],oa[samp,line,*]’
;(Tried the 3I4.0 thing.)
r = execute(fs)

i=i+1 ; add 1 to index

endfor
endfor

;Figure out corners

print
print,’Xmin/Xmax = ‘, mc(2), oax[ns-1]+ps(0)
print,’Ymin/Ymax = ‘, mc(3), oay[nl-1]+ps(1)
print

writeu,u3,oa ;oa=outarray
;calculate time that classification took
endtime = systime(1)
time = endtime – starttime
minutes = time/60
printf,u2,’Total time for classification: ‘,minutes, ‘ minutes’
print,’done, total time for classification: ‘,minutes, ‘ minutes’

close,u2,u1,u3
FREE_LUN,u2,u1,u3

print
print,’Xmin/Xmax = ‘, mc(2), oax[ns-1]+ps(0)
print,’Ymin/Ymax = ‘, mc(3), oay[nl-1]+ps(1)
print

end

_______________________

matzke
Geo Member

Leave a Reply