Image linear regression script

Image linear regression script

posted 06-11-2001 03:09 AM
This script takes two input images/data cubes — e.g., a cube derived from DMSP consisting of a series of detected fire counts from different DN thresholds, and another cube derived from Landsat consisting of a series of images resulting from GENERATE_BURNIMAGES (processed burn maps, producing count, total area, & avg area burned).
Anyway, it takes each band of each image, extracts the DN values from a user-created ROI that covers the data of both images (i.e., the extent of the original Landsat image), and then regresses the DNs against each other.

Further, it takes the DMSP count image (1st image input), rounds the values to whole numbers, counts up the number of observations of each burn count, averages all of the corresponding points on the landsat ETM-derived bands, and regresses those. This is done because the high per-pixel variability in the ETM-derived images is difficult to predict, while earlier work indicated that the mean and median could be predicted if DMSP count values (e.g., < 6 fires detected in the pixel) were high enough to get a decent number of observations (e.g., >1% of the data — often there is only 1 pixel with e.g. 14 detected fires).

At the moment I’ve got 8 bands of DMSP threshold running vs. 33 bands of Landsat- derived products, times 3 ways of cooking the above data; this is 792 regressions performed with minimal effort.

Further editions may run more than linear models, may apply filters to the DMSP dataset, etc. Pretty cool.

Nick

code:

; PRO GENERATE_BURNIMAGES_LINREGS3.PRO
;
; by Nick Matzke, May 2001
;
; PURPOSE: To regress the 99 images from generate_burnimages4.pro:
; – 3 resolutions: 1 km, 3 km, 5 km
; – 3 methods: total area, total count, average area
; – 11 burnscar minimum size thresholds
; = 99 images
;
; Against the DMSP observed fire count image (& that image resampled to 3 km & 5km)
;
; POSSIBLE FUTURE ADDITIONS:
;
; INPUTS:
;
; THE GENERAL PLAN FOR THIS SCRIPT:
; 1)
;
; HISTORY
; This script started out as a modification of GENERATE_BURNIMAGES4.PRO, with some
; code swiped from Kerry Halligan’s REGRESS_SAR_BIO.PRO.
pro generate_burnimage_linregs3
; Program title (must match filename)
; Started May, 2001.
; Written by Nick Matzke

; 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

; OPEN INPUT DMSP IMAGE
; 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 observed fire count image’, fid=fid1, $
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, fid1, data_type=data_type, nb=nb1, interleave=interleave, sname=sname1, $
fname=fname1;, h_map=h_map, handle_value, h_map, map_info
map_info = envi_get_map_info(fid=fid1)

;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

;
; OPEN INPUT PROCESSED BUNCH OF LANDSAT IMAGES
; 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 Landsat-derived data cube’, fid=fid2, $
pos=pos2, dims=dims2, /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, fid2, data_type=data_type2, nb=nb2, interleave=interleave2, sname=sname2, $
fname=fname2;, h_map=h_map, handle_value, h_map, map_info
;handle_value, h_map, map_info
map_info2 = envi_get_map_info(fid=fid2)
;calculate the number of lines and samples of the input image (this deals with subsetting)
ns2 = dims2(2)-dims2(1)+1
nl2 = dims2(4)-dims2(3)+1

;SELECT NAME OF OUTPUT FILE
;use pickfile to select the output image name (starting path will be same as input ascii rulesfile)
outimagename=dialog_pickfile(title=’Output File Name (DO NOT add *.img)’, get_path=gp, path=filename)
; don’t put apostrophes inside IDL quotes or they will truncate the quote!!!)
nbands=nb1

; CHECK TO MAKE SURE IMAGES HAVE SAME DIMENSIONS

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

; Run test to see if subsetted:
print,’Running test to see if subsetting needed…’
if ( (ns ne ns2) OR (nl ne nl2) ) then begin ; ‘ne’ = not equal
print,’Your DMSP image and your Landsat-derived image dont have the same’
print,’dimensions. You ninny!!’
print,’Actually, do this subsetting manually to save time (ENVI Tools –> Resize data).’
print,’The DMSP image dimensions are:’
print,’ DMSP # of samples = ‘,ns
print,’ DMSP # of lines = ‘,nl
print,’ Landsat-derived # of samples = ‘,ns2
print,’ Landsat-derived # of lines = ‘,nl2
print,’Ending program; re-run when ready.’

stop
;ending program early (saves time if subsetting done manually)
endif else begin
print,’Test passed.’
print,’Entering DMSP data into ENVI…’
inputimage = envi_get_data(fid=fid1, dims=[-1,0,ns-1,0,nl-1], pos=pos)
envi_file_query, fid1, bnames=bandname1
ENVI_ENTER_DATA, inputimage, descrip=’Max coarsenable subset of 30-m image’, $
file_type=0, map_info=map_info, xstart=0, ystart=0, bnames=bandname1, r_fid=fid3
print,’ENVI_ENTER_DATA for DMSP done.’

print,’Entering Landsat-derived data cube into ENVI…’
inputimage = envi_get_data(fid=fid2, dims=[-1,0,ns2-1,0,nl2-1], pos=pos2)

envi_file_query, fid2, bnames=bandname2
ENVI_ENTER_DATA, inputimage, descrip=’Max coarsenable subset of 30-m image’, $
file_type=0, map_info=map_info2, xstart=0, ystart=0, bnames=bandname2, r_fid=fid4
print,’ENVI_ENTER_DATA for Landsat-derived data cube done.’

print, ‘bandname1 = ‘
print, bandname1
print, ‘bandname2 = ‘
print, bandname2
endelse
;36 chars or so…

;For-loop: strip the ‘)’ off of the string in the first field and convert to float

length1=findgen(nb1)
length2=findgen(nb2)
bandnamearray1=sindgen(nb1)
bandnamearray2=sindgen(nb2)

;print,’length1 =’,length1
;help,length1
;print,’length2 =’,length2
;help,length2
;print,’nb1 =’,nb1
;print,’nb2 =’,nb2

;print,’help ban’

;Read in end of band name string for output purposes
i=0
for i=0,nb1-1 do begin
length1[i] = strlen(bandname1[i])
bandnamearray1[i] = strmid(bandname1[i],length1[i]-36,length1[i]-2)
endfor

;Read in end of band name string for output purposes
i=0
for i=0,nb2-1 do begin
length2[i] = strlen(bandname2[i])
bandnamearray2[i] = strmid(bandname2[i],length2[i]-36,length2[i]-2)
endfor
i=0
;print, ‘bandnamearray1 = ‘
;print, bandnamearray1
;print, ‘bandnamearray2 = ‘
;print, bandnamearray2

;FROM HELP FILE
;envi_open_file, ‘/data/test.image’, r_fid=fid
;envi_file_query, fid, ns=ns, nl=nl

;FROM HALLIGAN’S RANGECALC_part2.pro
;select the images
;envi_file_query, sr1_fid, nl=nl, ns=ns, h_map=h_map, fname=fname1, sname=sname1
;map_info = envi_get_map_info(fid=first_fid)

;ESTABLISH THRESHOLDS IN CODE
;(still need ’em for purposes of identifying results)
;
; Create an array of thresholds, to produce a group of output images
numberthresholds=byte(11)
burnareathresh = findgen(numberthresholds)
burnareathresh[0] = 1000000
burnareathresh[1] = 500000
burnareathresh[2] = 250000
burnareathresh[3] = 100000
burnareathresh[4] = 50000
burnareathresh[5] = 25000
burnareathresh[6] = 10000
burnareathresh[7] = 5000
burnareathresh[8] = 2500
burnareathresh[9] = 1000
burnareathresh[10] = 50 ;Using 0 screws up count images

; 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+’-regression_metadata.txt’, /GET_LUN

;BE CAREFUL WITH openw NOT TO BLAST A GOOD FILE

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

;print to file (= printf)
printf,u1, ‘Metadata file for GENERATE_BURNIMAGES_LINREGS1.PRO run’
printf,u1, ‘Time program run performed: ‘,SYSTIME()
printf,u1, ‘ ‘
printf,u1, ‘Name of inputted DMSP fire count image: ‘,sname1
printf,u1, ‘ fname = ‘,fname1
printf,u1, ‘ ‘
printf,u1, ‘Name of inputted Landsat-derived data cube: ‘,sname2
printf,u1, ‘ fname = ‘,fname2
printf,u1, ‘ ‘
printf,u1, ‘Number of lines : ‘, nl
printf,u1, ‘Number of samples : ‘, ns
printf,u1, ‘Number of bands : ‘, nb2
printf,u1, ‘Data type : ‘, data_type
printf,u1, ‘ ‘
printf,u1, ‘OUTPUT DATA’
printf,u1, ‘Thresholds:’
for a=0,numberthresholds-1 do begin
printf,u1, ‘Area threshold Output image names’
printf,u1, burnareathresh[a],’ ‘,outimagename,’-gt’,burnareathresh[a],’.img’
endfor

;print to screen
print,’ ‘
print,’===========================================================================’
print, ‘Metadata file for GENERATE_BURNIMAGES_LINREGS1.PRO run’
print, ‘Time program run performed: ‘,SYSTIME()
print, ‘ ‘
print, ‘Name of inputted DMSP fire count image: ‘,sname1
print, ‘ fname = ‘,fname1
print, ‘ ‘
print, ‘Name of inputted Landsat-derived data cube: ‘,sname2
print, ‘ fname = ‘,fname2
print, ‘ ‘
print, ‘Number of lines : ‘, nl
print, ‘Number of samples : ‘, ns
print, ‘Number of bands : ‘, nb2
print, ‘Data type : ‘, data_type
print, ‘ ‘
print, ‘OUTPUT DATA’
print, ‘Thresholds:’
for a=0,numberthresholds-1 do begin
print, ‘Area threshold Output image names’
print, burnareathresh[a],’ ‘,outimagename,’-gt’,burnareathresh[a],’.img’
endfor
print,’ ‘
print,’===========================================================================’
print,’ ‘

;This shouldn’t be necessary if above test was passed
;ns=ns_fine_max
;nl=nl_fine_max

;SELECT ROIs TO COMPARE
;use roi_ids to retrieve all of the defined roi’s for the image
;for the DMSP image
roi_ids1 = envi_get_roi_ids(fid=fid1, roi_names=roi_names1)
if (roi_ids1[0] eq -1) then begin
print, ‘No regions associated with the selected file’
return
endif

;build a compound widget for selecting the desired ROIs
base1 = widget_auto_base(title=’ROI Selection’)
wm1 = widget_multi(base1, list=roi_names1, uvalue=’list1′, /auto)
result1 = auto_wid_mng(base1)

;count is an output of the ROI selection widget that is the number of ROIs selected
;ptr is equal to count for all all valid lists (lists with value of 1 rather than -1)
ptr1 = where(result1.list1 eq 1, count)

;determine the total number of ROI points and make an array that is nb x number of points
tot_pnts1=intarr(count)
for i=0,count-1 do begin
roi_addr1 = envi_get_roi(roi_ids1[i])
tot_pnts1(i) = n_elements(roi_addr1)
endfor
total_pnts1 = total(tot_pnts1)
final_array1 = make_array(nb1,total_pnts1,type=data_type)

;begin the loops that retrieve all of the data by band for each ROI
;for each RIO
for i=0,count-1 do begin
;query the ROI for the number of points
roi_addr1 = envi_get_roi(roi_ids1[i])
points1 = n_elements(roi_addr1)
;create a temp. array to store the band data for all that is nb by points
temp_array1 = make_array(nb1,points1,type=data_type)
;for each band of the image
for j=0, nb1-1 do begin
;retrieve the image data for all points in the image
data1 = envi_get_roi_data(roi_ids1[ptr1[i]], fid=fid1, pos=pos[j])
;fill columns with data for each band
temp_array1(j,0 :points1-1) = data1(0 :points1-1)
endfor
;copy values from temp_array to the appropriate cells of final_array
if (i ne 0) then begin
ystart1 = total(tot_pnts1(0:i-1))
endif else begin
ystart1 = 0
endelse
yend1 = (total(tot_pnts1(0:i))-1)
final_array1(0:nb1-1,ystart1:yend1) = temp_array1
endfor

print,’help,final_array1=’
help,final_array1
print,’print,final_array1=’
;print,final_array1
print,’ending’
;stop
print,’nb1=’,nb1

print,”
print,’Getting second ROI dataset’
print,”
;SELECT ROIs TO COMPARE
;use roi_ids to retrieve all of the defined roi’s for the image
;for the DMSP image
roi_ids2 = envi_get_roi_ids(fid=fid2, roi_names=roi_names2)
if (roi_ids2[0] eq -1) then begin
print, ‘No regions associated with the selected file’
return
endif

;build a compound widget for selecting the desired ROIs
base2 = widget_auto_base(title=’ROI Selection’)
wm2 = widget_multi(base2, list=roi_names2, uvalue=’list2′, /auto)
result2 = auto_wid_mng(base2)

;count is an output of the ROI selection widget that is the number of ROIs selected
;ptr is equal to count for all all valid lists (lists with value of 1 rather than -1)
ptr2 = where(result2.list2 eq 1, count)

;determine the total number of ROI points and make an array that is nb x number of points
tot_pnts2=intarr(count)
for i=0,count-1 do begin
roi_addr2 = envi_get_roi(roi_ids2[i])
tot_pnts2(i) = n_elements(roi_addr2)
endfor
total_pnts2 = total(tot_pnts2)
final_array2 = make_array(nb2,total_pnts2,type=data_type)

;begin the loops that retrieve all of the data by band for each ROI
;for each ROI
for i=0,count-1 do begin
;query the ROI for the number of points
roi_addr2 = envi_get_roi(roi_ids2[i])
points2 = n_elements(roi_addr2)
;create a temp. array to store the band data for all that is nb by points
temp_array2 = make_array(nb2,points2,type=data_type)
;for each band of the image
for j=0, nb2-1 do begin
;retrieve the image data for all points in the image
data2 = envi_get_roi_data(roi_ids2[ptr2[i]], fid=fid2, pos=pos2[j])
;fill columns with data for each band
temp_array2(j,0 :points2-1) = data2(0 :points2-1)
endfor
;copy values from temp_array to the appropriate cells of final_array
if (i ne 0) then begin
ystart2 = total(tot_pnts2(0:i-1))
endif else begin
ystart2 = 0
endelse
yend2 = (total(tot_pnts2(0:i))-1)
final_array2(0:nb2-1,ystart2:yend2) = temp_array2
endfor

print,’help,final_array2=’
help,final_array2
;print,’print,final_array2=’
;print,final_array2

print,’points1=’,points1
print,’points2=’,points2
print,’nb2=’,nb2

Thanks, Nick

[This message has been edited by matzke (edited 06-11-2001).]

________________

matzke
Geo Member

Image linear regression script

Permalink Submitted by admin on Wed, 01/19/2011 – 12:31

posted 06-11-2001 03:09 AM
…here’s the rest:
code:
;NOW, REGRESS EACH BAND OF 33-BAND ROI AGAINST DMSP ROI
;PREDICTOR/INDEPENDENT VARIABLE = x = DMSP
;DEPENDENT VARIABLE = y = Landsat-derived

;from Halligan’s program
;;open the output file and print the header row
openw,u2,outimagename+’-regression_metadata.txt’, /get_lun

;PRINT STATEMENT TITLES
;set number of columns in code
ncols=11

dasharray=strarr(ncols)
for i=0,ncols-1 do begin
dasharray[i]=(“——————————“)
endfor
dash=string(“————————————————“)
dashes = ‘print,FORMAT=”(1A,””,’ + string(ncols) + ‘A10,””,A)”,dash,dasharray

,dash’
r = execute(dashes)
print
print
print,’REGRESSION OUTPUT RESULTS’
printf,u2,’REGRESSION OUTPUT RESULTS’
r = execute(dashes)
;ptf=’f,u2′ ;ptf=print to file (use this code if necessary
;titlestring=’bnum1 sname1 band1 bnum2 sname2 band2 slope y-int r r^2 ftest rmul sigma status’

;fs = ‘print,FORMAT=”($,1A,” “,’ + string(ncols) + ‘(A10, :, ” “),””,A)”,titlestring’
;r = execute(fs)
;r = execute(dashes)

nnames=4
;Create array for string names
outputstringarray=sindgen(4)
outputstringarray[0]=’File name 1 ‘ ;sname = short filname?
outputstringarray[1]=’ Band name 1 ‘
outputstringarray[2]=’ File name 2 ‘
outputstringarray[3]=’Band name 2 ‘
;Create array for output numbers
;const,r,rsquare,ftest,rmul,sigma,status
outputstr2array=sindgen(ncols)
outputstr2array[0]=’ File 1 band # ‘
outputstr2array[1]=’File 2 band # ‘
outputstr2array[2]=’result[0] ‘
outputstr2array[3]=’const ‘
outputstr2array[4]=’r ‘
outputstr2array[5]=’rsquare ‘
outputstr2array[6]=’ftest ‘
outputstr2array[7]=’rmul ‘
outputstr2array[8]=’sigma ‘
outputstr2array[9]=’status ‘
outputstr2array[10]=’maxsubs’
;outputnumarray[9]=p ; put probability calc here…

;put in ncols in for 11, somehow…

fs = ‘print,FORMAT=”($,4A45,” “,’ + string(nnames) + ‘(11A15, :, ” “),””,A15)”,outputstringarray

,outputstr2array
,”chisq1′
r = execute(fs)
print,”
;printf,u2,FORMAT='($,(4A39))’,’sname1′,’bandnamearray1′,’sname2′,’bandnamearray2′
;printf,u2,FORMAT='($,(10A18))’,’Prod. type ‘,’threshhold ‘,’slope ‘,$
;’ y-intercept ‘,’r ‘,’ r^2 ‘,’ftest ‘,’rmul ‘,’sigma ‘,’status ‘;’p-value’
;printf,u2,FORMAT='(A)’,’chisq’
;print,FORMAT='($,(4A39))’,’sname1′,’bandnamearray1′,’sname2′,’bandnamearray2′
;print,FORMAT='($,(10A18))’,’Prod. type’,’threshhold’,’slope’,’y-intercept’,’r’,’r^2′,’ftest’,’rmul’,’sigma’,’status’;’p-value’
;print,FORMAT='(A18)’,’chisq’

;loop through the bands and run regressions for each solo band with each bio param
count=0
yb=findgen(points2,1)
type=’avg_area’
threshnum=0
threshnumfract=float(0)
for i=0,nb1-1 do begin ;for 1 band
x=float(final_array1(i,0 :points1-1))
for j=0,nb2-1 do begin ;for 33 bands
count=count+1
;for k=0,points2-1 do begin
yb[0 :points2-1,0]=float(final_array2[j,0 :points2-1])
;endfor
y=yb
;print,’x’
;help,x
;print,’y’
;help,y
weights=REPLICATE(1.0,N_ELEMENTS(y))
result1=REGRESS(x,y,MEASURE_ERRORS=error1,YFIT=yfit1,CONST=const1,SIGMA=sigma1,FTEST=ftest1,$
CORRELATION=r1,MCORRELATION=rmul1,CHISQ=chisq1,STATUS=status1,/DOUBLE);/RELATIVE_WEIGHT)
rsquare1=r1*r1
p1=0

;singular array error = array of all zeros; remove from data cube…

;if (j eq 11) then begin
; type=’ttl_area’
;endif else begin
; z=0
;endelse
;if (j eq 22) then begin
; type=’ttl_count’
;endif else begin
; z=0
;endelse

maxsubscript=-9999

;PRINT STATEMENTS
;Create array for string names
outputstringarray=sindgen(4)
outputstringarray[0]=sname1 ;sname = short filname?
outputstringarray[1]=bandnamearray1[i]
outputstringarray[2]=sname2
outputstringarray[3]=bandnamearray2[j]

;Create array for output numbers
;const,r,rsquare,ftest,rmul,sigma,status
outputnumarray=findgen(ncols)
outputnumarray[0]=i+1
outputnumarray[1]=j+1
outputnumarray[2]=result1[0]
outputnumarray[3]=const1
outputnumarray[4]=r1
outputnumarray[5]=rsquare1
outputnumarray[6]=ftest1
outputnumarray[7]=rmul1
outputnumarray[8]=sigma1
outputnumarray[9]=status1
outputnumarray[10]=maxsubscript+1

;for row=0,ncols-1 do begin
;build a print format statement using the ncols variable
;for the number of classes to print pixels for, first variable is
;a string (trainpixelsarray)
fs = ‘print,FORMAT=”($,4A45,” “,’ + string(ncols) + ‘(F, :, ” “),””,E20.10)”,outputstringarray

,outputnumarray
,chisq1’
r = execute(fs)
print,” ; execute the print string
;endfor
;r = execute(dashes)
;fs = ‘print,FORMAT=”(1A,” “,’ + string(ncols) + ‘(I10, :, ” “),””,I)”,”BEST”,TOTAL(sumrows)’
;r = execute(fs)
;r = execute(dashes)

;printf,u2,FORMAT='($,(4A39))’,sname1,bandnamearray1[i],sname2,bandnamearray2[j]
;printf,u2,FORMAT='($,(10A18))’,type,burnareathresh[threshnum],result[0],const,r,rsquare,ftest,rmul,sigma,status
;printf,u2,FORMAT='(A)’,chisq
;print,FORMAT='($,(4A39))’,sname1,bandnamearray1[i],sname2,bandnamearray2[j]
;print,FORMAT='($,(10A18))’,type,burnareathresh[threshnum],result[0],const,r,rsquare,ftest,rmul,sigma,status
;print,FORMAT='(A18)’,chisq
;threshnumbfract=1/3+threshnumfract
;threshnum=FLOOR(threshnumfract)
;junk
;,bandname1(i),bandname2(j)

endfor
endfor

;Generate mean stats.

;help,final_array1
;help,final_array2

;for taking only the points with enough samples:
n_threshold = float(0.01*total_pnts2)

count=0
yb=findgen(points2,1)
type=’avg_area’
threshnum=0
threshnumfract=float(0)
for i=0,nb1-1 do begin ;for 1 band
x=float(final_array1(i,0 :points1-1))

for j=0,nb2-1 do begin ;for 33 bands
count=count+1
;for k=0,points2-1 do begin
yb[0 :points2-1,0]=float(final_array2[j,0 :points2-1])
;endfor
y=yb
;print,’x’
;help,x
;print,’y’
;help,y

x=round(x) ;bins x data (counts) into integers

uniqxvalues = x[UNIQ(x, SORT(x))]
;help,’help uniqxvalues’
;help,uniqxvalues

;create two arrays, for uniq x values & mean corresponding y values
blah=WHERE(uniqxvalues, uniqxvaluescount) ; get count

maxsubscript2=uniqxvaluescount

;print,’Number of unique x-values in file 1, band ‘,i+1,’ = ‘,uniqxvaluescount

xarraytemp=uniqxvalues ;findgen(uniqxvaluescount)
yarraytemp=findgen(uniqxvaluescount)
;xvalcount=findgen(uniqxvaluescount)

xvalcount=float(0)

for k=0,uniqxvaluescount-1 do begin
xnum=uniqxvalues[k] ; Get the current uniq value
;print,’xnum = ‘, xnum

yarraytemp[k]=MEAN(y[WHERE(x eq xnum,xvalcount)])
;xvalcount[k]=xvalcounttemp

;help,xvalcount
;help,n_threshold

if (xvalcount ge n_threshold) then begin
maxsubscript3=k
endif else begin
zo=0
endelse

; print,’xarray[‘,k,’] = ‘,xarray[k],’ yarray[‘,k,’] = ‘,yarray[k],’ count(xarray[‘,k,’]) = ‘,xvalcount[k]
;get y total for that xnum
endfor

;switching rows & cols to make regression work…
xarray=findgen(1,uniqxvaluescount)
yarray=findgen(uniqxvaluescount,1)
xarray[0,0:uniqxvaluescount-1]=uniqxvalues[0:uniqxvaluescount-1]

yarray[0:uniqxvaluescount-1,0]=float(yarraytemp[0:uniqxvaluescount-1])
;endfor
y=yb

weights=REPLICATE(1.0,N_ELEMENTS(y))
result2=REGRESS(xarray,yarray,MEASURE_ERRORS=errors2,YFIT=yfit2,CONST=const2,SIGMA=sigma2,$
FTEST=ftest2,CORRELATION=r2,MCORRELATION=rmul2,CHISQ=chisq2,STATUS=status2,/DOUBLE);/RELATIVE_WEIGHT)
rsquare2=r2*r2
p2=0

;number of points threshold
;
;underthresh=WHERE((xvalcount lt n_threshold), uniqxvaluescount)

xsubarray=xarray[0,0:maxsubscript3]
ysubarray=yarray[0:maxsubscript3,0]

;print,’k=’,k,’ maxsubscript= ‘,maxsubscript

result3=REGRESS(xsubarray,ysubarray,MEASURE_ERRORS=errors3,YFIT=yfit3,CONST=const3,SIGMA=sigma3,$
FTEST=ftest3,CORRELATION=r3,MCORRELATION=rmul3,CHISQ=chisq3,STATUS=status3,/DOUBLE);/RELATIVE_WEIGHT)
rsquare3=r3*r3
p3=0

;print,’xarray = , yarray =’
;print,xarray,’ ‘,yarray

;top

;PRINT STATEMENTS
;Create array for string names
outputstringarray=sindgen(4)
outputstringarray[0]=sname1 ;sname = short filname?
outputstringarray[1]=bandnamearray1[i]
outputstringarray[2]=sname2
outputstringarray[3]=bandnamearray2[j]

;Create array for output numbers
;const,r,rsquare,ftest,rmul,sigma,status
outputnumarray=findgen(ncols)
outputnumarray[0]=i+1
outputnumarray[1]=j+1
outputnumarray[2]=result2[0]
outputnumarray[3]=const2
outputnumarray[4]=r2
outputnumarray[5]=rsquare2
outputnumarray[6]=ftest2
outputnumarray[7]=rmul2
outputnumarray[8]=sigma2
outputnumarray[9]=status2
outputnumarray[10]=maxsubscript2+1; put probability calc here…

;for row=0,ncols-1 do begin
;build a print format statement using the ncols variable
;for the number of classes to print pixels for, first variable is
;a string (trainpixelsarray)
fs = ‘print,FORMAT=”($,4A45,” “,’ + string(ncols) + ‘(F, :, ” “),””,E20.10)”,outputstringarray

,outputnumarray
,chisq2’
r = execute(fs)
print,” ; execute the print string
;endfor
;r = execute(dashes)
;fs = ‘print,FORMAT=”(1A,” “,’ + string(ncols) + ‘(I10, :, ” “),””,I)”,”BEST”,TOTAL(sumrows)’
;r = execute(fs)
;r = execute(dashes)

;PRINT STATEMENTS
;Create array for string names
outputstringarray=sindgen(4)
outputstringarray[0]=sname1 ;sname = short filname?
outputstringarray[1]=bandnamearray1[i]
outputstringarray[2]=sname2
outputstringarray[3]=bandnamearray2[j]

;Create array for output numbers
;const,r,rsquare,ftest,rmul,sigma,status
outputnumarray=findgen(ncols)
outputnumarray[0]=i+1
outputnumarray[1]=j+1
outputnumarray[2]=result3[0]
outputnumarray[3]=const3
outputnumarray[4]=r3
outputnumarray[5]=rsquare3
outputnumarray[6]=ftest3
outputnumarray[7]=rmul3
outputnumarray[8]=sigma3
outputnumarray[9]=status3
outputnumarray[10]=maxsubscript3+1

;for row=0,ncols-1 do begin
;build a print format statement using the ncols variable
;for the number of classes to print pixels for, first variable is
;a string (trainpixelsarray)
fs = ‘print,FORMAT=”($,4A45,” “,’ + string(ncols) + ‘(F, :, ” “),””,E20.10)”,outputstringarray

,outputnumarray
,chisq3′
r = execute(fs)
print,” ; execute the print string
;printf,u2,FORMAT='($,(4A39))’,sname1,bandnamearray1[i],sname2,bandnamearray2[j]
;printf,u2,FORMAT='($,(10A18))’,type,burnareathresh[threshnum],result[0],const,r,rsquare,ftest,rmul,sigma,status
;printf,u2,FORMAT='(A)’,chisq
;print,FORMAT='($,(4A39))’,sname1,bandnamearray1[i],sname2,bandnamearray2[j]
;print,FORMAT='($,(10A18))’,type,burnareathresh[threshnum],result[0],const,r,rsquare,ftest,rmul,sigma,status
;print,FORMAT='(A18)’,chisq
;threshnumbfract=1/3+threshnumfract
;threshnum=FLOOR(threshnumfract)
;junk
;,bandname1(i),bandname2(j)

endfor
endfor

count1=count
print,’Number of single-band regressions : ‘,count1

print,’ending’
;stop
envi_file_mng, id=fid3, /remove, /delete
envi_file_mng, id=fid4, /remove, /delete
close,u2,fid1,fid2,fid3,fid4
FREE_LUN,u2,fid1,fid2,fid3,fid4

print,’all calcs done set, beginning file writing’
;calculate time that classification took
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

[This message has been edited by matzke (edited 06-11-2001).]

________________

matzke
Geo Member

Leave a Reply