GENERATE_BURNIMAGES_LINREGS5.PRO (backup)

GENERATE_BURNIMAGES_LINREGS5.PRO (backup)

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

; PRO GENERATE_BURNIMAGES_LINREGS5.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-
;
; INPUTS:
; Via ENVI_SELECT and PICKFILE, the user must select:
; 1- The first image/image cube to be compared (an image or meta file- — generally DMSP-derived
; 2- The second image/image cube to be compared — generally Landsat-derived.
; * All of the images being compared must have the same dimensions, and must be using the same
; ROI.
; 3- Pick the name of an output file. No extension (*.txt, etc.- should be added.
; 4 & 5- The user must then select the desired ROI for each image. The ROI must be created
; beforehand, and ought to cover just the part of the images containing data (e.g., Landsat
; will have nodata areas in the corners-.
;
;OUTPUTS:
; 1- A metadata file, containing information on the input images, time required for the run,
; image dimensions, etc.
; 2- A text file containing tab-delimited regression results, with headings at the top.
; 3- The above will also be output to the screen.
;
; THE GENERAL PLAN FOR THIS SCRIPT:
; 1- Get inputs as above.
; 2- Open files for output
; 3- Calculate metadata info
; 4- Extract data from the image cubes for the ROIs
; 5- Begin regressions (using the REGRESS command and perhaps others)
; Image-to-image regression:
; – Plain regression, no specified weights/measure_errors variable
; – Plain regression, gaussian weights/measure_errors variable
; – Plain regression, Poisson weights/measure_errors variable
; – SQRT both images, then do REGRESS with gaussian distribution
; – SQRT both images, then do REGRESS with poisson distribution
; Count-to-mean regression:
; – Here, all the DMSP pixels are binned into integers, and the means of the
; corresponding Landsat ETM-derived pixels are regressed against ’em.
; – Repeat each of the above for REGRESSing the means
;
; Count-to-mean regression, outliers removed…
; – Repeat each of the above for REGRESSing the means, minus outliers (often those w/ low n counts for x)
;
;
; POSSIBLE FUTURE ADDITIONS:
; Regressions:
; – median
; – other distributions
; – multiple regressions with multiple predictors, e.g.:
; – land cover
; – observed clouds
; – predicted precip. for the month
; – elevation/aspect/slope
; – variance?
;
; 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_linregs5
; 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
openw,u2,outimagename+’-results.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

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

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

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,u1,’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
outputtype=strarr(1)
;Create array for string names
outputstringarray=sindgen(4)
outputstringarray[0]=’Filename1 ‘ ;sname = short filname?
outputstringarray[1]=’ Bandname1 ‘
outputstringarray[2]=’ Filename2 ‘
outputstringarray[3]=’Bandname2 ‘
;Create array for output numbers
;const,r,rsquare,ftest,rmul,sigma,status
outputstr2array=sindgen(ncols)
outputstr2array[0]=’ File1band# ‘
outputstr2array[1]=’File2band# ‘
outputstr2array[2]=’result[0] ‘
outputstr2array[3]=’const ‘
outputstr2array[4]=’r ‘
outputstr2array[5]=’rsquare ‘
outputstr2array[6]=’pslope ‘
outputstr2array[7]=’ftest ‘
outputstr2array[8]=’rmul ‘
outputstr2array[9]=’sigma ‘
outputstr2array[10]=’status ‘
outputstr2array[11]=’maxsubs’
outputtype=’model_type’
;outputnumarray[9]=p ; put probability calc here…

;put in ncols in for 11, somehow…

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

,outputstr2array
,”chisq “,outputtype’
r = execute(fs)
print,”
fs = ‘printf,u2,FORMAT=”($,A,” “,4A45,” “,’ + string(nnames) + ‘(12A15, :, ” “),””,A15,A)”,” count”,outputstringarray

,outputstr2array
,”chisq “,outputtype’
r = execute(fs)
printf,u2,”
;printf,u1,FORMAT='($,(4A39))’,’sname1′,’bandnamearray1′,’sname2′,’bandnamearray2′
;printf,u1,FORMAT='($,(10A18))’,’Prod. type ‘,’threshhold ‘,’slope ‘,$
;’ y-intercept ‘,’r ‘,’ r^2 ‘,’ftest ‘,’rmul ‘,’sigma ‘,’status ‘;’p-value’
;printf,u1,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
regcount=0
outputtype=strarr(1)

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

;RUN ORIGINAL, UN-NEATIFIED REGRESSION…
;for k=0,points2-1 do begin
yb[0:points2-1,0]=float(final_array2[j,0:points2-1])
;endfor
y=yb

xmax=MAX(x)
ymax=MAX(y)

;Inversion error due to singular array trap
if ((xmax ne 0) AND (ymax ne 0)) then begin

errorarray=REPLICATE(1.0,N_ELEMENTS(y)); not used in original regress
maxsubscript=-10000
outputtype=’orig_image_run’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,x,y,error1,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;Alternative method:
; p = linfit(x,y, sigma=dp) ;; DP are parameter uncertainties;
; p returns slope & intercept…
; r2 = correlate(x,y) …actually gives you r, not r2

; print,’r2 =’,r2
; print,’p = ‘,p

; stop

;ORIGINAL REGRESSION, GAUSSIAN DISTRIBUTION
errorarray=REPLICATE(0.5,N_ELEMENTS(y)); for Gaussian (=normal) distribution
maxsubscript=-10000
outputtype=’gauss_image_run’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,x,y,errorarray,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, POISSON DISTRIBUTION
errorarray1=findgen(N_ELEMENTS(y))
errorarray1[0:points2-1]=sqrt(abs(y[0:points2-1,0]))
errorarray1=errorarray1+0.0000001 ; To avoid ‘singular array error’ with zeros
maxsubscript=-10000
outputtype=’poisson_image_run’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,x,y,errorarray1,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;commented out — previous regression style
;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
;regcount=regcount+1
;singular array error = array of all zeros; remove from data cube…

;Do it again for square root transform
xtemp=x
ytemp=y

x=(sqrt(x))
y=(sqrt(y))

errorarray=REPLICATE(1.0,N_ELEMENTS(y)); not used in original regress
maxsubscript=-10001
outputtype=’orig_image_run_sqrt_xy’
blah=call_function(‘RUNREGRESS’,x,y,error1,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, GAUSSIAN DISTRIBUTION
errorarray=REPLICATE(0.5,N_ELEMENTS(y)); for Gaussian (=normal) distribution
maxsubscript=-9001
outputtype=’gauss_image_run_sqrt_xy’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,x,y,errorarray,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, POISSON DISTRIBUTION
errorarray1=findgen(N_ELEMENTS(y))
errorarray1[0:points2-1]=sqrt(abs(y[0:points2-1,0]))
errorarray1=errorarray1+0.0000001 ; To avoid ‘singular array error’ with zeros
maxsubscript=-8001
outputtype=’poisson_image_run_sqrt_xy’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,x,y,errorarray1,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;Do it again for log transform of just x
;Scratched: do it again for log transform

;x=(xtemp+1)
;y=(ytemp+1)

x=(xtemp)
y=(sqrt(ytemp))

errorarray=REPLICATE(1.0,N_ELEMENTS(y)); not used in original regress
maxsubscript=-10001
outputtype=’orig_image_run_sqrt_y_only’
blah=call_function(‘RUNREGRESS’,x,y,error1,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, GAUSSIAN DISTRIBUTION
errorarray=REPLICATE(0.5,N_ELEMENTS(y)); for Gaussian (=normal) distribution
maxsubscript=-9001
outputtype=’gauss_image_run_sqrt_y_only’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,x,y,errorarray,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, POISSON DISTRIBUTION
errorarray1=findgen(N_ELEMENTS(y))
errorarray1[0:points2-1]=sqrt(abs(y[0:points2-1,0]))
errorarray1=errorarray1+0.0000001 ; To avoid ‘singular array error’ with zeros
maxsubscript=-8001
outputtype=’poisson_image_run_sqrt_y_only’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,x,y,errorarray1,maxsubscript,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

endif else begin
print,’Error — array of all zeros’
endelse
endfor
endfor

;GENERATE THE MEAN STATS…

;help,final_array1
;help,final_array2

;for taking only the points with enough samples:
n_threshold = float(0.01*total_pnts2)
if n_threshold LT 5 then begin
n_threshold = 5
endif else begin
zo=1
endelse

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
yb[0:points2-1,0]=float(final_array2[j,0:points2-1])
y=yb
x=round(x) ;bins x data (counts) into integers

xmax=MAX(x)
ymax=MAX(y)

;Inversion error due to singular array trap
if ((xmax ne 0) AND (ymax ne 0)) then begin

uniqxvalues = x[UNIQ(x, SORT(x))]

;create two arrays, for uniq x values & mean corresponding y values
blah=WHERE(uniqxvalues, uniqxvaluescount) ; get count
if uniqxvaluescount EQ 0 then begin
uniqxvaluescount=1
endif else begin
zo=0
endelse

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

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

;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))
if maxsubscript2 GT 3 then begin ;ERROR TRAP
errorarray=REPLICATE(1.0,N_ELEMENTS(y)); not used in original regress
outputtype=’mean_orig_regress-all’
blah=call_function(‘RUNREGRESS’,xarray,yarray,errors2,maxsubscript2,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, GAUSSIAN DISTRIBUTION
errorarray=REPLICATE(0.5,N_ELEMENTS(y)); for Gaussian (=normal) distribution
outputtype=’mean_gauss_regress-all’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,xarray,yarray,errorarray,maxsubscript2,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, POISSON DISTRIBUTION
errorarray1=findgen(N_ELEMENTS(y))
errorarray1[0:points2-1]=sqrt(abs(y[0:points2-1,0]))
errorarray1=errorarray1+0.0000001 ; To avoid ‘singular array error’ with zeros
outputtype=’mean_poisson_regress-all’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,xarray,yarray,errorarray1,maxsubscript2,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

endif else begin
print,’Too few points to regress (<3); setting Status to -1′
status2=-1
const2=0
r2=0
rsquare2=0
ftest2=0
rmul2=0
sigma2=0
result2=0
chisq2=0
endelse

;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

if ((maxsubscript3 GT 3) AND (TOTAL(ysubarray) NE 0) ) then begin ;ERROR TRAP
errorarray=REPLICATE(1.0,N_ELEMENTS(y)); not used in original regress
outputtype=’orig_mean_regress-subset’
blah=call_function(‘RUNREGRESS’,xsubarray,ysubarray,errors2,maxsubscript3,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, GAUSSIAN DISTRIBUTION
errorarray=REPLICATE(0.5,N_ELEMENTS(y)); for Gaussian (=normal) distribution
outputtype=’gauss_mean_regress-subset’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,xsubarray,ysubarray,errorarray,maxsubscript3,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)

;ORIGINAL REGRESSION, POISSON DISTRIBUTION
errorarray1=findgen(N_ELEMENTS(y))
errorarray1[0:points2-1]=sqrt(abs(y[0:points2-1,0]))
errorarray1=errorarray1+0.0000001 ; To avoid ‘singular array error’ with zeros
outputtype=’poisson_mean_regress-subset’
;blah=RUNREGRESS(x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount)
blah=call_function(‘RUNREGRESS’,xsubarray,ysubarray,errorarray1,maxsubscript3,outputtype,u2,i,j,$
sname1,bandnamearray1,sname2,bandnamearray2,regcount)
endif else begin
print,’Too few points to regress (<3); setting Status to -1′
status3=-1
const3=0
r3=0
rsquare3=0
ftest3=0
rmul3=0
sigma3=0
result3=0
chisq3=0
endelse
endif else begin
print,’Error — array of all zeros’
endelse

endfor
endfor

count1=count
print,’Number of single-band regressions : ‘,count1
printf,u1,’Number of single-band regressions : ‘,count1
printf,u1,’ending’
print,’ending’
;stop
envi_file_mng, id=fid3, /remove, /delete
envi_file_mng, id=fid4, /remove, /delete

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,u2
FREE_LUN,u1,u2

end

;START ‘RUN REGRESSION’ FUNCTION
function RUNREGRESS,x,y,errorarray,maxsubscript,outputtype,resultsfile,i,j,sname1,bandnamearray1,sname2,bandnamearray2,regcount

;print,’x’
;help,x
;print,’y’
;help,y
;weights=REPLICATE(1.0,N_ELEMENTS(y))
;error1=sqrt(abs(y))

result1=REGRESS(x,y,MEASURE_ERRORS=errorarray,YFIT=yfit1,CONST=const1,SIGMA=sigma1,FTEST=ftest1,$
CORRELATION=r1,MCORRELATION=rmul1,CHISQ=chisq1,STATUS=status1,/DOUBLE);/RELATIVE_WEIGHT)
rsquare1=r1*r1
pslope=yfit1/sigma1

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

regcount=regcount+1
; maxsubscript=-9999

ncols=12

;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]=pslope ;probability of slope being significant
outputnumarray[7]=ftest1
outputnumarray[8]=rmul1
outputnumarray[9]=sigma1
outputnumarray[10]=status1
outputnumarray[11]=maxsubscript+1
;outputtype=’image_x_vs_y’

;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)
fs1 = ‘print,FORMAT=”($,I,” “,4A45,” “,’ + string(ncols) + ‘(F, :, ” “),””,E20.6,” “,A)”,regcount,outputstringarray

,outputnumarray
,chisq1,outputtype’
r = execute(fs1)
fs1 = ‘printf,resultsfile,FORMAT=”($,I,” “,4A45,” “,’ + string(ncols) + ‘(F, :, ” “),””,E20.6,” “,A)”,regcount,outputstringarray
,outputnumarray
,chisq1,outputtype’
r = execute(fs1)
printf,resultsfile,” ; execute the print string
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)
return,1
end

_______________

matzke
Geo Member

Leave a Reply