posted 02-27-2001 04:01 PM
Hi all,
Here’s my DMSP totals program — adds up the number of “hits” over a time series of DMSP fire data, given a user-defined (in the code) threshold.

Messy but functional.

Thanks, Nick


; Initial hopes: to read in a large time-series of DMSP data and
; output a mean & variance image, and perhaps further useful products.
; well as ENVI_SELECT, and DIALOG_PICKFILE to allow the user to select input/output
; files
; Started Feb 22, 2001.
; written by Nick Matzke, Feb 22, 2001.
pro dmspedit

; use the envi_select widget 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

;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

;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 Classified Image’, get_path=gp, path=filename)


;set up and write the envi header for output image
bandname[0]=’binary_class (‘+imagename+’)’
envi_setup_head, data_type=1, $
descrip=’Output image from DMSP-EDIT.PRO using DMSP time series input’, $
file_type=0, fname=outimage+’.hdr’, interleave=0, nl=nl, ns=ns, nb=nbands, map_info=map_info, $
xstart=dims(1), ystart=dims(3), bnames=bandname, /WRITE
;nb = nb. Original nb=1

;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+’-mean.img’, /GET_LUN
openw,u2,outimage+’-mean.txt’, /GET_LUN

date_time = systime()
;print to file (= printf)
printf,u2,’Metadata file for BINARY_CLASS 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 : ‘, 1
printf,u2, ‘Data type : byte’

;print to screen
print,’Metadata file for BINARY_CLASS 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, ‘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, ‘ ‘

;make an array with the same dimensions as original image
inputimage = findgen(ns,nl,nbands)
S = SIZE(arr)

;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,nb-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,’importing envi image band number: ‘,band+1 ;envivars[i]

;create an array to hold the classified image
threshximage = findgen(ns,nl,nb)
;thresh20image = make_array(ns,nl,nb,/FLOAT)
;thresh30image = make_array(ns,nl,nb,/FLOAT)

print,’all parameters set, beginning classification’

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

threshx=float(1) ;setting threshhold
;for band=0,nb-1 do begin
threshximage[*,*,0:nb-1]=(inputimage[*,*,0:nb-1]/threshx) ;)[*,*,*] [(0:ns-1),(0:nl-1),band]

fireimagesx = bindgen(ns,nl,nbands)
;fireimage20 = make_array(ns,nl,nb,/byte)
;fireimage30 = make_array(ns,nl,nb,/byte)

;for band=0,nb-1 do begin
fireimagesx[*,*,0:nb-1] = bytscl(threshximage[*,*,0:nb-1],max=1.0000,min=0.9999,/NAN,TOP=1)

totalfiresx = bindgen(ns,nl)
totalfiresx = TOTAL(fireimagesx,3,/NAN)


;****double check that this is getting the right bands****
;image[*,*,i] = ENVI_GET_DATA(dims=dims,fid=fid,pos=idlvars[i])
;print,’importing envi image band number: ‘,envivars[i]

;create an array to hold the output mean image
;fltmeanarray = make_array(ns,nl, /float)
;fltmeanarray = AVERAGE(image)

print,’all parameters set, beginning DMSP time series calculations’

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

;begin a loop that processes the ENVI image pixel by pixel
;for i=0,nl-1 do begin

; show percent completed to screen
; percent = float(float(i)/float(nl)*100)
; print,FORMAT='($(1A))’,’percent completed: ‘
; print,FORMAT='(1F4.0)’,percent

;write the array to the output image
print,’writing output image file’
;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’


; AVERAGE function sucked from help files & modified.

[This message has been edited by matzke (edited 03-16-2001).]

Geo Member

One comment on “

  1. posted 03-21-2001 04:15 PM
    Minor bug discovered:
    Sometimes when you run the S+ decision tree algorithm, one of the classes may be completely classified as another class. If this happens, sruleconv will bomb because there will be (e.g.) 9 columns in the decision tree input text file fractions, but only 8 names in the tree.

    (so actually, the problem is with S+ & your training data, not my script. Ha!)

    Geo Member

Leave a Reply