;; just a nice program to plot the output of ndvi vs century in a couple of different ways

pro plotnvc, datafile, outfile, byteit=byteit, nclasses=nclasses, nb=nb, calndvi=calndvi, plotboth=plotboth, nodetail=nodetail
  
  old=!D.name
  set_plot, 'ps'
  device,file= outfile ,xsize=7.5, ysize=10, xoff=.5, yoff=.5, /inches, $
         set_font='Times', /tt_font, font_size=14
  
  outfile=string(datafile, '3.dat')
  
  openr, un, /get, datafile
  openw, oun, /get, outfile
  
  if not keyword_set(plotboth) then plotboth=0
  if not keyword_set(nclasses) then nclasses=12
  if not keyword_set(nb) then nb=8
  count=0l
  
  save=!p.multi
  if plotboth eq 1 then begin
     !p.multi=[0,4,round(nclasses/2.)]
  endif else begin
;     !p.multi=[0,3,round(nclasses/3. + 0.4)]
     !p.multi=[0,3,4]
  endelse
  
  for class=0, nclasses-1 do begin
     if not eof(un) then begin
        
        readu, un, count

        IF count GT 10000000l THEN BEGIN 
           print, 'ERROR, number fo points is too large!'
           print, '   Count = ', strcompress(count)
           device, /close
           !p.multi=save
           set_plot, old
           
           close, un, oun
           free_lun, un, oun
           return
        ENDIF
 

        ndvi=intarr(count, nb)
        if keyword_set(byteit) then begin
           cent=bytarr(count, nb)
        endif else $
          cent=fltarr(count, nb)
        
        
        readu, un, ndvi
        readu, un, cent

;; set up arrays for averaging into bins
        Cdata=indgen(150)
        Ndata=intarr(150)
        
        if plotboth eq 1 then begin
           for curval=0, 149 do begin
              thiscount=0
              for year=0, nb-1 do begin
                 index=where(fix(cent(*,year)) eq curval)
                 
                 if index(0) ne -1 and n_elements(index) gt 200 then begin
                    thiscount=thiscount+1
                    Ndata(curval)=Ndata(curval)+fix(mean(ndvi(index)))
                 endif
              endfor
;	if curval mod 15 eq 0 then print, fix(float(curval)/1.49), '% of class', class+1
              if thiscount ne 0 then $
                Ndata(curval) = Ndata(curval)/thiscount
           endfor
           
           index=where(Ndata ne 0)
           
           if index(0) ne -1 then begin
              tmp=Ndata(index);/255.
              plot, Cdata(index), tmp, psym=1, $
                    yrange=[min(tmp)-10, max(tmp)+10], /ystyle, $
                    xrange=[min(Cdata(index))-5, max(Cdata(index))+5], /xstyle, $
                    xtitle='Century Aglivc+0.5Stded', ytitle='NDVI'
              
;;Compute the best regression on the means for each century bin, and plot
              if n_elements(index) gt 2 then begin
                 result=regress(transpose(Cdata(index)), tmp, replicate(1.0, n_elements(tmp)),$
                                yfit, const, sigma, ftst, rval)
                 print, 'COMPUTED BY CENTURY BIN FOR CLASS : ', class+1
                 print, 'fit',result
                 print, 'const', const
;		print, 'sigma', sigma
                 print, 'F test', ftst
                 print, 'p = ', F_PDF(ftst, 1, n_elements(tmp)-1)
                 print, 'R sq', (rval^2)
                 
                 writeu, oun, 150, Cdata, Ndata
                 y=[const, result*200 + const]
                 oplot, [0,200], y
              endif else print, 'Could not regress by bin for class:',class+1
              
           endif $ ;;no elements to plot!
           else begin
              dummy=[[0,1],[1,0]]
              plot, dummy(0,*),dummy(1,*)
              print, 'ERROR, no valid points in class, ', class+1
              print, '   Computing valid points, there are a total of ', count, ' points in this class'
              for curval=0, 149 do begin
                 thiscount=0
                 for year=0, nb-1 do begin
                    index=where(fix(cent(*,year)) eq curval)
                    
                    if index(0) ne -1 and n_elements(index) le 200 then begin
                       print, n_elements(index)
                    endif
                 endfor
              endfor
           endelse
        endif ;;plotboth
        
        Cdata=fltarr(nb)
        Ndata=fltarr(nb)
        
        for year=0, nb-1 do begin 
           Cdata(year)=mean(cent(*,year))
           Ndata(year)=mean(ndvi(*,year))
        endfor
        
        
        title=string('Class'+strcompress(class+1))
        index = where(Ndata ne 0)
        if index(0) ne -1 then begin
;		tmp=float(Ndata(index))/255.
           tmp = Ndata(index)
           xrng=[min(Cdata(index))-5, max(Cdata(index))+5]
           yrng=[min(tmp)-100, max(tmp)+100]
           xsz=xrng[1]-xrng[0]
           ysz=yrng[1]-yrng[0]
;offset the top so we have room to print r^2 and pvals
           yrng[1]=yrng[1]+ysz/10  

           plot, Cdata(index), tmp, psym=1, $
                 yrange=yrng, /ystyle, $
                 xrange=xrng, /xstyle, $
                 xtitle='Century Aglivc+0.5Stded', ytitle='NDVI', title=title, $
                 thick=3        ;, xtick=1, ytick=3
           
;    i=linfit(Cdata(index), tmp)
;	y=[i(0), i(1)*200 + i(0)]
;	oplot, [0,200], y
;    print, i
           if n_elements(index) gt 2 then begin
              result=regress(transpose(Cdata(index)), tmp, replicate(1.0, n_elements(tmp)),$
                             yfit, const, sigma, ftst, rval)
              print, 'COMPUTED BY YEARLY MEAN FOR CLASS : ', class+1
              print, 'Number of Points= ', count
              print, 'fit               ',result
              print, 'const             ', const
;			print, 'sigma', sigma
;              print, 'F test            ', ftst
              pval=F_PDF(ftst, 1, n_elements(tmp)-1)
              print, 'p =               ', pval
              print, 'R sq =            ', (rval^2)

              print, ''
              print, ''
              print, ''
;              print, Ndata(index)
              writeu, oun, n_elements(index), Cdata(index), Ndata(index)
              y=[const, result*200 + const]
              oplot, [0,200], y

              xyouts, xrng[0]+xsz/20., yrng[1]-ysz/20., $
                      string("R^2=", strcompress(rval^2), $
                             '                      npixels=', strcompress(count)), $
                      charsize=0.5
              xyouts, xrng[0]+xsz/20., yrng[1]-ysz/20.-ysz/20, $
                      string("p  =", strcompress(pval)), charsize=0.5

              IF NOT keyword_set(nodetail) THEN BEGIN
                 plot, Cdata(index), /ystyle
                 plot, Ndata(index), /ystyle
              ENDIF

           endif else print, 'Could not regress by yearly mean for class:',class+1
        endif
        
     endif                      ;check for null classes
     
  endfor
  
  device, /close
  !p.multi=save
  set_plot, old
  
  close, un, oun
  free_lun, un, oun
  
  
end