pro make_images

; ToDo: 
; add JJA, DJF line for +- 23.5 degrees (swap in south). Color blue + red

; file = 'v2.temperature.inv.txt'
; ;               012  7   10      42   48   56
; fmt = '          A3, A5, A3,     A32, F7.2,F8.2,I5, I4,  A1,I5, A2,A2,A2,I2,A1,I2,A16,A1'
; readfmt,file,fmt,icc,wmo,wmo_mod,name,lat, lon, ele,tele,p, pop,tp,v, lo,co,a, ds,veg,bi,skipline=38
; save, icc,wmo,wmo_mod,name,lat,lon,ele,tele,p,pop,tp,v,lo,co,a,ds,veg,bi
restore
name = STRCOMPRESS(name)
restore, 'data_struct.sav'

cred = 253
cblue = 10

;for i=0,n_elements(icc)-1 do begin
for i=5000,5010 do begin
    idx =  where( s.icc eq icc[i] AND $
                  s.wmo eq wmo[i] AND $
                  s.wmo_mod eq wmo_mod[i], count )
    if count eq 0 then continue
    station = icc[i] + wmo[i] + wmo_mod[i]
    
    set_plot,'ps'
    device, /color, bits=8
    
    ; get the min/max of all the temperatures, so we can pre-set the Y range
    t_all = [[s[idx].jja,s[idx].djf,s[idx].ann]]
    gd = where( t_all ne 999.90 AND t_all NE 0, cnt )
    if cnt eq 0 then continue 
    trange = minmax(t_all[gd])

    ; plot annual data
    plotg, s[idx].year[gd], s[idx].ann[gd], $
           /yno, yrange=trange, $
           xrange=[1880,2010], /xst,psym=2, xminor=2, $
           xtitle='Time (Years)', ytitle='Temperature (Degrees C)', $
           chars=1.2, charth=3, title=name[i]
    
    ; 1st order fit to ANN
    gd = where( s[idx].ann ne 999.90 AND s[idx].year NE 0 )
    line = indgen(n_elements(s[idx].year))
    r1 = poly_fit( line[gd], s[idx].ann[gd], 1 )
    fit1 = r1[0] + r1[1]*line
    oplot, s[idx].year[gd], fit1[gd], color=0
    
    ; if non-tropical, then 1st order fit to JJA and DJF
     if abs( lat[i] ) GE 23.26 then begin
        ; 1st order fit to JJA
        gd = where( s[idx].JJA ne 999.90 AND s[idx].year NE 0, count )
        if count ne 0 then begin
            r1 = poly_fit( line[gd], s[idx].jja[gd], 1 )
            fit1 = r1[0] + r1[1]*line
            color=(lat[i] ge 0)? cred: cblue
            oplot, s[idx].year[gd], s[idx].JJA[gd], color=color, psym=2, th=2
            oplot, s[idx].year[gd], fit1[gd], color=color, thick=2
        endif
        ; 1st order fit to DJF
        gd = where( s[idx].DJF ne 999.90 AND s[idx].year NE 0, count )
        if count ne 0 then begin
            r1 = poly_fit( line[gd], s[idx].djf[gd], 1 )
            fit1 = r1[0] + r1[1]*line
            color=(lat[i] lt 0)? cred: cblue
            oplot, s[idx].year[gd], s[idx].DJF[gd], color=color, psym=2, th=2
            oplot, s[idx].year[gd], fit1[gd], color=color, thick=2
       endif
    endif

;     r2 = poly_fit( line[gd], ann[gd], 2 )
;     fit2 = r2[0]+r2[1]*line+r2[2]*line^2
;     oplot, year, fit2, color=254
    
    device,/close
    set_plot,'x'
    spawn, 'convert idl.ps  ./images/' + station + '.png'
;    spawn, 'open ./images/' + station + '.png'
endfor

end

make_images
end