DMSP_INDICATOR1.PRO (backup)

DMSP_INDICATOR1.PRO (backup)

posted 07-24-2001 01:40 PM
code:

; PRO DMSP_INDICATOR1.PRO
;
; Basic idea:
;
; Read in a DMSP “Cumulative Fire Count” image, and turn each band into
; an indicator-variable image suitable for ANOVA analysis. Output headers
; for said images as well.
;
;
; well as ENVI_SELECT, and DIALOG_PICKFILE to allow the user to select input/output
; files
;
;
; Started July 20th, 2001.
; written by Nick Matzke, July, 2001.
pro dmsp_indicator1

; use envi_check_save to load the save files with the registration and band math routines
ENVI_CHECK_SAVE, /registration
ENVI_CHECK_SAVE, /envi_math
ENVI_CHECK_SAVE, /utility
ENVI_CHECK_SAVE, /spectral

; INPUT FILES & GET HEADER INFO
; 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 DMSP FIRECOUNT IMAGE’, fid=fid1, pos=pos1, dims=dims1, /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, fid1, data_type=data_type1, nb=nb1, interleave=interleave1, sname=sname1, fname=fname1, $
xstart=xstart1, ystart=ystart1, pixel_size=pixel_size1;, h_map=h_map
map_info1 = envi_get_map_info(fid=fid1)
proj1 = envi_get_projection(fid=fid1)

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

;CORNER COORDINATES; CHANGE FOR EACH RUN ON A DIFFERENT IMAGE
;Establish upper left corner coords for map projection (should be done automatically, but
;DO BY HAND FOR NOW!

leftpixel=double(0)
toppixel=double(0)

p158r73fullcoord_left=double(1241703.972)
p158r73fullcoord_top=double(3000163.484)
p158r73temp2_left=double(1291953.972)
p158r73temp2_top=double(2921863.484)

;CHANGE FOR EACH RUN ON A DIFFERENT IMAGE
leftcoord=p158r73fullcoord_left
topcoord=p158r73fullcoord_top

;use pickfile to select the output image name (starting path will be same as input ascii rulesfile)
outimagename=dialog_pickfile(title=’Enter Path/row/version (e.g. p158r73v2) of Output Image (DO NOT add *.img)’, get_path=gp, path=fullpath)
; don’t put apostrophes inside IDL quotes or they will truncate the quote!!!)

;figure out max fine dimensions that can be fit into ouput coarse image

rescoarse=1000 ; low-resolution, in meters

;resfine=30 ; high-resolution, in meters
;finefract=(float(resfine)/rescoarse)
;ns_coarse_max=(floor(ns*finefract))
;ns_fine_max=CEIL(ns_coarse_max/(float(finefract)))
;nl_coarse_max=(floor(nl*resfine/rescoarse))
;nl_fine_max=CEIL(nl_coarse_max/(float(finefract)))

;print,’ns_coarse_max = ‘,ns_coarse_max
;print,’nl_coarse_max = ‘,nl_coarse_max
;print,’ ‘
;print,’ns_fine_max = ‘,ns_fine_max
;print,’nl_fine_max = ‘,nl_fine_max
;print,’ ‘

;SET UP INITIAL 1-KM HEADER
;MAKE HEADERS FOR 1-KM IMAGES
pixel_size1=[rescoarse,rescoarse] ; This pixel size is for the header; the other (‘ps’)
;is for the projection. Intelligent system, no? NO.

;
ps1=dblarr(2)
ps1(0)=rescoarse
ps1(1)=rescoarse
mc1=dblarr(4)
mc1(0)=leftpixel
mc1(1)=toppixel
mc1(2)=leftcoord
mc1(3)=topcoord

; THESE ARE THE ORIGINAL THRESHOLDS USED IN THE dmspedit2.pro PROGRAM
; THAT PRODUCED THE dmsp count IMAGE
thresharraycount=nbands1 ;number of thresholds to be used (8 in this case)
thresharray=findgen(thresharraycount) ;holding the desired DN thresholds
thresharray[0]=1
thresharray[1]=5
thresharray[2]=10
thresharray[3]=20
thresharray[4]=35
thresharray[5]=50
thresharray[6]=60
thresharray[7]=63

; SET THE MINIMUM COUNT THRESHOLD FOR PRODUCING INDICATOR IMAGES
; Create an array of thresholds, to produce a group of output images
countthresh_count=byte(4)
countthresh = findgen(countthresh_count)
countthresh[0] = 0
countthresh[1] = 1
countthresh[2] = 5
countthresh[3] = 10

;countthresh_max = float(10^37) ; set in program code; not currently used
;maximum in floating point (single precision) = +/- 10^38

;Here, we are going to put all of the requisite names and numbers for headers into arrays

si1=strarr(9) ;stringinfo-1km
fi1=fltarr(7) ;floatinfo-1km

;Get header info for 1-km images.
map_info_1km = envi_map_info_create(proj=ENVI_GET_PROJECTION(fid=fid1), ps=ps1, mc=mc1)
pixel_size1=pixel_size1
nb=1

fi1[0] = ns
fi1[1] = nl
fi1[2] = 1
fi1[3] = data_type1
fi1[4] = 0 ;interleave by band (methinks)
fi1[5] = xstart1
fi1[6] = ystart1

si1[0] = ‘-‘ ; prefix; was “_ind-”
si1[1] = outimagename ; filename user entered at prompt
si1[2] = ‘DN_gt’ ; DN value threshold in original DMSP count image
si1[3] = ‘threshstring’ ; string of min. area threshold
si1[4] = ‘-1km’ ; resolution suffix
si1[5] = ‘.img’ ; .img suffix
si1[6] = ‘.hdr’ ; .hdr suffix
si1[7] = ‘-ind_gt’ ; bandnames
si1[8] = ‘Generate burnimages output’ ; description of file/bandnames

;bandname=strarr(1)
;bandname[0]=’binary_class (‘+imagename+’)’
;envi_setup_head, data_type=1, $
;descrip=’Output image from DMSP_INDICATOR1.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

; WRITING METADATA FILES
; open a metadata file that has the same name as the output image
; with a *-metadata.txt at the end and start filling it
; One metadata file for the whole shebang.

openw,u1,outimagename+’-metadata.txt’, /GET_LUN

;get time/date that program run begins
starttime = systime(1)

;print to file (= printf)
printf,u1, ‘Metadata file for DMSP_INDICATOR1.PRO run’
printf,u1, ‘Time program run performed: ‘,SYSTIME()
printf,u1, ‘ ‘
printf,u1, ‘Input DMSP burncount image: ‘,sname1
printf,u1, ‘ fname = ‘,fname1
printf,u1, ‘Input image pathname: ‘,gp
printf,u1, ‘ ‘
printf,u1, ‘Number of lines : ‘, nl
printf,u1, ‘Number of samples : ‘, ns
printf,u1, ‘Number of bands : ‘, nb
printf,u1, ‘Data type : ‘, data_type1
printf,u1, ‘ ‘
printf,u1, ‘OUTPUT IMAGES’
printf,u1, ‘Thresholds for binary-indicator images, DNs used =’
print, thresharray
for a=0,thresharraycount-1 do begin
printf,u1, ‘Area threshold Output image names’
printf,u1, a,’ ‘,outimagename,’-gt’,thresharray[a],’.img’
endfor
printf,u1, ‘Output mask image name : ‘,outimagename,’-tempmask.img’
printf,u1, ‘Output masked image name : ‘,outimagename,’-tempmasked.img’
printf,u1, ‘Number of lines : ‘, nl
printf,u1, ‘Number of samples : ‘, ns
printf,u1, ‘Number of bands : ‘, 1
printf,u1, ‘Data type : byte’

;print to screen
print,’ ‘
print,’===========================================================================’
;print to file (= printf)
print, ‘Metadata file for DMSP_INDICATOR1.PRO run’
print, ‘Time program run performed: ‘,SYSTIME()
print, ‘ ‘
print, ‘Input DMSP burncount image: ‘,sname1
print, ‘ fname = ‘,fname1
print, ‘Input image pathname: ‘,gp
print, ‘ ‘
print, ‘Number of lines : ‘, nl
print, ‘Number of samples : ‘, ns
print, ‘Number of bands : ‘, nb
print, ‘Data type : ‘, data_type1
print, ‘ ‘
print, ‘OUTPUT IMAGES’
print, ‘Thresholds for binary-indicator images, DNs used =’
print, thresharray
for a=0,thresharraycount-1 do begin
print, ‘Area threshold Output image names’
print, a,’ ‘,outimagename,’-gt-‘,thresharray[a],’.img’
endfor
print, ‘Output mask image name : ‘,outimagename,’-tempmask.img’
print, ‘Output masked image name : ‘,outimagename,’-tempmasked.img’
print, ‘Number of lines : ‘, nl
print, ‘Number of samples : ‘, ns
print, ‘Number of bands : ‘, 1
print, ‘Data type : byte’
print,’===========================================================================’

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

print,’Entering ENVI data…’

;inputimage = envi_get_data(fid=fid1, dims=[-1,0,ns-1,0,nl-1], pos=pos1[bigband])
;ENVI_ENTER_DATA, inputimage, descrip=’Max coarsenable subset of 30-m image’, $
;file_type=0, map_info=map_info, xstart=0, ystart=0, bnames=bandname, r_fid=fid2;
;print,’ENVI_ENTER_DATA done.’

;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)

;can’t use t=0 ’cause you can’t divide by 0…
for t=0,thresharraycount-1 do begin
thresharraystring=string(thresharray[t],FORMAT='(I3.3)’)
print,’importing envi image band number: ‘,t+1 ;envivars[i]

inputimage = ENVI_GET_DATA(dims=dims1,fid=fid1,pos=t) ;idlvars[i])

for i=0,countthresh_count-1 do begin
print,’export output band number: ‘,i+1 ;envivars[i]

zo = inputimage
;for some reason inputimage is wiped when stuck in CREATE_INDICATOR_MAPS
RUNFUNCTION=call_function(‘CREATE_INDICATOR_MAPS’, zo, countthresh[i], $
thresharraystring, outputimage_indicator, fi1, si1, map_info_1km, pixel_size1)
endfor

; openw,u2,outimagename+’-fire_hits-DN’+thresharraystring+’.img’, /GET_LUN
; threshx=float(thresharray[t]) ;setting threshhold (count this number or greater)
; help,threshximage
; help,inputimage
;for band=0,nb-1 do begin
; threshximage[*,*,0:nb-1]=(inputimage[*,*,0:nb-1]/threshx) ;)[*,*,*] [(0:ns-1),(0:nl-1),band]
;endfor
; help,threshximage

; fireimagesx = bindgen(ns,nl,nbands)

;Get fires gt or equal to threshold
;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)
;endfor

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

;bytscl

;****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]
;endfor

;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’
; print,’DN threshold = ‘,thresharray[t]

;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
;endfor

;write the array to the output image
print,’writing output image file’
;for
; writeu,u1,totalfiresx
;calculate time that classification took
; close,u1
; FREE_LUN,u1
endfor

endtime = systime(1)
time = endtime – starttime
minutes = time/60
printf,u1,’Total time for classification: ‘,minutes, ‘ minutes’
print,’done, total time for classification: ‘,minutes, ‘ minutes’

close,u1
FREE_LUN,u1
end

; AVERAGE function sucked from help files & modified.
FUNCTION AVERAGE, arr
RETURN, TOTAL(arr)/N_ELEMENTS(arr)
END

;ACCESSORY FUNCTIONS

;START CREATE_INDICATOR_MAPS FUNCTION
function CREATE_INDICATOR_MAPS,inputindicatorimage, countthresh, thresharraystring, $
outputindicatorimage, fi, si, map_info, pixel_size

;print, ‘Creating indicator (binary 1/0) maps…’
;si[1] has the PATHNAME, darn it!

ns=fi[0]
nl=fi[1]

thresh=float(countthresh)
countthreshstring = string(thresh,FORMAT='(I3.3)’)

; Since this function is for indicator variables, the output image
; should have this “indicator variable file” marker
si[0] = ‘-‘ ; Then put this non-blank in the header name
test1=1 ; This will pass the info to the header-creator function

;Stick filename in “name” variable
;name=(si[1]+si[0]+si[2]+thresharraystring+’-‘+countthreshstring+si[4]+si[5])
name=string(si[1]+si[0]+si[2]+thresharraystring+si[7]+countthreshstring+si[4]+si[5])
openw,u100,(name),/GET_LUN

; enter data into ENVI
ENVI_ENTER_DATA, inputindicatorimage, descrip=’Input data for thresholding to indicator var’, $
file_type=0, map_info=map_info, xstart=0, ystart=0, bnames=bandname, r_fid=fid30;

; Replace all nonzero pixels with ‘1’

maskexp='(b1 gt ‘+countthreshstring+’)’
math_doit, fid=fid30, exp=maskexp, out_name=’temp-imask.img’, $
dims=[-1,0,ns-1,0,nl-1], pos=0, r_fid=fid31

;Create array for fine-res subsetted, thresholded input image
outputindicatorimage = envi_get_data(fid=fid31, dims=[-1,0,ns-1,0,nl-1], pos=0)

writeu,u100,outputindicatorimage

; WRITE HEADER FOR INDICATOR IMAGE
;CALLING APPLY_HEADER (creating 1/0 binary indicator variables…) FUNCTION
RUNFUNCTION=call_function(‘APPLY_HEADER’, fi, si, countthreshstring, map_info, pixel_size, thresharraystring, test1)

;return these to default state
si[0] = ‘-‘
test1=0

close,u100
FREE_LUN,u100

envi_file_mng, id=fid30, /remove, /delete
envi_file_mng, id=fid31, /remove, /delete

return,1
end

;inputimage, countthresh[i], $
; thresharraystring, outputimage_indicator, fi1, si1, map_info_1km, pixel_size11

;START APPLY_HEADER
function APPLY_HEADER,fi, si, countthreshstring, map_info, pixel_size, thresharraystring, test

if (test eq 1) then begin ; If user specifies “indicator variable file”
si[0] = ‘-‘ ; Then put this non-blank in the header name
fi[3] = 1 ; 4 = data_type = byte
endif else begin
si[0] = ‘_raw’
fi[3] = 4 ; 4 = data_type = Floating-point
endelse

;APPLY HEADER TO 1-KM IMAGES

name=string(si[1]+si[0]+si[2]+thresharraystring+si[7]+countthreshstring+si[4]+si[6])

envi_setup_head, fname=name, $
ns=fi[0],$
nl=fi[1], $
nb=fi[2], $
data_type=fi[3], $
interleave=fi[4], $
bnames=[si[0]+si[2]+thresharraystring+si[7]+countthreshstring+si[4]], $
descrip=si[8], $
xstart=fi[5], ystart=fi[6], map_info=map_info, pixel_size=pixel_size, /write
;DEF_STRETCH={dstretch_struct, type:0L, vals:fltarr(4)}

; print,’fi’
; print,fi
; print,”

; print,’si’
; print,si
; print,”

;CALLING CREATE_INDICATOR_MAPS (creating 1/0 binary indicator variables…) FUNCTION
;RUNFUNCTION=call_function(‘APPLY_HEADER’, fi, si, map_info, pixel_size)

;, , $
;outimagenametemp1 + ‘_ind_’ + outimagenametemp2a + thresharraystring, $
;’-1km.img’, outcoarseimage_indicator,ns_coarse_max,nl_coarse_max)

; Return to defaults
si[0] = ‘-‘
fi[3] = 4

return,1
end

____________________

matzke
Geo Member

Leave a Reply