;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; calc_healthtrendV2.pro ; ; mary.martin@unh.edu ; (603) 862-4508 ; 20091203 ; ; USAGE: at envi or idl prompt, type '.run calchealthtrend.pro' ; this will report that the function "healthtrend" has been compiled ; Then type: a=healthtrend("inputfilename","outputfilename","ancillaryfilename",numrows,numcols,numbands) ; be sure to enclose the file path/name with quotes. ; ; Here is a sample of the command - remember, on a windows machine, the slashes will be backslashes ; if the path names are long, and the commands wrap the line, you can use the "$" to signify ; that the command continues on the next line ; ; a=healthtrend("/net/nfs/gromit/d1/proj/nsrclandsat/testsubset/1989-2009stack", $ ; "/net/nfs/gromit/d1/proj/nsrclandsat/testsubset/healthtrend", $ ; "/net/nfs/gromit/d1/proj/nsrclandsat/testsubset/ancillarydata", 1021,916,21) ; ; INPUT FILE REQUIREMENTS: an input file containing health data for a series of years ; this file should be datatype=float and storage order=BSQ ; the different year-layers will be the 'bands' in the file. Be sure to buffer the ; missing years with a 'band' of zeros when you layerstack the years together. You can ; create this buffer layer by using band math and multiplying one of your stacked layers ; by zero. This is better than 'generate_test_data', since it preserves all of the ; map info. ; ; OUTPUT FILE: the healthtrend output file will hold a single band which is the slope of the fit ; line between year and health. ; The ancillary output file holds number_of_data_years, largest_gap_in_timeseries, and currently has ; a placeholder for chisq and prob - We need to better understand stdev in the linfit.pro routine ; to get the latter two to work - Jen, we need your stats expertise here, and a little testing with a small ; sample array. ; The authors make no representations about the ; suitability of this software for any purpose. It is ; provided "as is" without express or implied warranty. ; Users of this program are strongly encouraged ; to make an independent assessment of its suitability before ; use (for instance, do calculations by hand for a few pixels - ; its a pain, but worth the effort). ; ; IF YOU NEED HELP: every action in this program is preceeded by a comment. Hopefully, ; this is enough information to allow you to follow the process, and maybe even ; modify or fix things if necessary. Feel free to contact me if you run into ; problems! Mary ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION healthtrend,infile,healthtrendfile,ancillaryfile,numrows, numcols,numbands ; read in the data cube and store in a 3-D array called 'image' print, "reading the bsq health cube" openr, un, infile, /get_lun ; initialize an array of floats to hold the input data image=fltarr(numcols, numrows,numbands) readu, un, image free_lun, un ; open the output file for writing openw, un2, healthtrendfile, /get_lun openw, un3, ancillaryfile, /get_lun ; create an array with the year numbers (the first year of data will be year zero) years=indgen(numbands) ; create an array to hold the output (slope of fit line) ; and a second array to fill with other attributes (like number of years of data, and ; the largest data gap in the time series, and chisq) slopearr=fltarr(numcols,numrows) ancillarydata=fltarr(numcols,numrows,4) ; loop through the rows and columns, extract the time series of heath data, and ; make a linear fit between health and year. write out the slope of the fit line to ; an image file ; these first two commented out for loops are useful for debugging - you might have ; to set them to do another range of pixels if these are masked out, but combining ; this with some other debugging print statements below will help to verify that ; any changes you make are doing what you want them to. ;for y=21, 22 do begin ; for x=21, 22 do begin for y=0, numrows-1 do begin for x=0, numcols-1 do begin ; pull out all years of data for the current pixel ; transform will take the 3D array and pull out the 1D time series for the pixel ts=reform(image(x,y,*)) ;print ;print, years,ts ; determine what years have health data, so that we can filter these out when ; we perform the fit - idl 'where' reports the indices in the array where the condition ; is tru w_ok=where(ts(*) gt 0) ;print, w_ok ;print, n_elements(w_ok) ; this little loop reads through the ts data for the pixel and determine what the ; max string of zeros is ; set two variables, one to store the current number of zeros in a row, and then ; a variable to keep track of which string of zeros is the largest zerocount=0 maxzerostring=0 for val=0, numbands-1 do begin ; if the ts value is zero then begin keeping count of sequential zeros if ts(val) eq 0 then begin zerocount=zerocount+1 ; only overwrite maxzerostring, if current count is greater than previous ; this is important, since zerocount is reset to zero when the loop encounters ; valid data if zerocount ge maxzerostring then begin maxzerostring=zerocount endif ;print, ts(val),zerocount,maxzerostring endif else begin ; reset the counter when valid data is present zerocount=0 ;print, ts(val),zerocount,maxzerostring endelse end ; store the maxzerostring in the second element of the ancillary array ancillarydata(x,y,1)=maxzerostring ;print, ts, zerocount, maxzerostring ;condition added so that it skips instances where there is no valid data ;for only the years with data, perform the linear fit if n_elements(w_ok) gt 1 then begin ;fit the year and health data for the current pixel ; ***NOTE***, this contains a placeholder for chisq and prob are here, but I think it needs ; stdev info - see linfit.pro to see if you can figure this out. I'm stuck on this. fit=LINFIT(years(w_ok),ts(w_ok),chisq=chisq, prob=prob) ; every 100 lines/cols, print out some data so to indicate progress ; this just lets you know that the process is running, could be commented out ; needs to be commented out if you are debugging by running just a few pixels if x MOD 100 eq 0 and y MOD 100 eq 0 then print,x,y,fit(1),chisq, prob ;print the slope of the model ; store the slope of the fit line for this pixel in slopearr slopearr(x,y)=fit(1) ; store the number of elements used to develop the fit in the first element of the ; ancillary array ancillarydata(x,y,0)=n_elements(w_ok) ancillarydata(x,y,2)=chisq ancillarydata(x,y,3)=FLOAT(prob) endif else begin ; see comments above slopearr(x,y)=999 ancillarydata(x,y,0)=n_elements(w_ok) ancillarydata(x,y,2)=999 ancillarydata(x,y,3)=999 endelse end end ; print out the slopearr when all calculations have been finished writeu, un2, slopearr writeu, un3, ancillarydata ; close the output file close,un2 free_lun, un2 end